Bug Summary

File:libraries/AMPTOOLS_AMPS/TwoPiAngles_primakoff.cc
Location:line 86, column 24
Description:Value stored to 'psi' is never read

Annotated Source Code

1
2#include <cassert>
3#include <iostream>
4#include <string>
5#include <sstream>
6#include <cstdlib>
7#include <complex.h>
8
9#include "TLorentzVector.h"
10#include "TLorentzRotation.h"
11
12#include "IUAmpTools/Kinematics.h"
13#include "AMPTOOLS_AMPS/TwoPiAngles_primakoff.h"
14#include "AMPTOOLS_AMPS/clebschGordan.h"
15#include "AMPTOOLS_AMPS/wignerD.h"
16
17TwoPiAngles_primakoff::TwoPiAngles_primakoff( const vector< string >& args ) :
18UserAmplitude< TwoPiAngles_primakoff >( args )
19{
20 assert( args.size() == 5 )((args.size() == 5) ? static_cast<void> (0) : __assert_fail
("args.size() == 5", "libraries/AMPTOOLS_AMPS/TwoPiAngles_primakoff.cc"
, 20, __PRETTY_FUNCTION__));
21
22 phipol = atof(args[0].c_str() )*3.14159/180.; // azimuthal angle of the photon polarization vector in the lab. Convert to radians.
23 polFrac = AmpParameter( args[1] ); // fraction of polarization (0-1)
24 m_rho = atoi( args[2].c_str() ); // Jz component of rho
25 PhaseFactor = AmpParameter( args[3] ); // prefix factor to amplitudes in computation
26 flat = atoi( args[4].c_str() ); // flat=1 uniform angles, flat=0 use YLMs
27
28 assert( ( phipol >= 0.) && (phipol <= 2*3.14159))((( phipol >= 0.) && (phipol <= 2*3.14159)) ? static_cast
<void> (0) : __assert_fail ("( phipol >= 0.) && (phipol <= 2*3.14159)"
, "libraries/AMPTOOLS_AMPS/TwoPiAngles_primakoff.cc", 28, __PRETTY_FUNCTION__
));
29 assert( ( polFrac >= 0 ) && ( polFrac <= 1 ) )((( polFrac >= 0 ) && ( polFrac <= 1 )) ? static_cast
<void> (0) : __assert_fail ("( polFrac >= 0 ) && ( polFrac <= 1 )"
, "libraries/AMPTOOLS_AMPS/TwoPiAngles_primakoff.cc", 29, __PRETTY_FUNCTION__
));
30 assert( ( m_rho == 1 ) || ( m_rho == 0 ) || ( m_rho == -1 ))((( m_rho == 1 ) || ( m_rho == 0 ) || ( m_rho == -1 )) ? static_cast
<void> (0) : __assert_fail ("( m_rho == 1 ) || ( m_rho == 0 ) || ( m_rho == -1 )"
, "libraries/AMPTOOLS_AMPS/TwoPiAngles_primakoff.cc", 30, __PRETTY_FUNCTION__
));
31 assert( ( PhaseFactor == 0 ) || ( PhaseFactor == 1 ) || ( PhaseFactor == 2 ) || ( PhaseFactor == 3 ))((( PhaseFactor == 0 ) || ( PhaseFactor == 1 ) || ( PhaseFactor
== 2 ) || ( PhaseFactor == 3 )) ? static_cast<void> (0
) : __assert_fail ("( PhaseFactor == 0 ) || ( PhaseFactor == 1 ) || ( PhaseFactor == 2 ) || ( PhaseFactor == 3 )"
, "libraries/AMPTOOLS_AMPS/TwoPiAngles_primakoff.cc", 31, __PRETTY_FUNCTION__
));
32 assert( (flat == 0) || (flat == 1) )(((flat == 0) || (flat == 1)) ? static_cast<void> (0) :
__assert_fail ("(flat == 0) || (flat == 1)", "libraries/AMPTOOLS_AMPS/TwoPiAngles_primakoff.cc"
, 32, __PRETTY_FUNCTION__));
33
34 // need to register any free parameters so the framework knows about them
35 registerParameter( polFrac );
36}
37
38
39complex< GDouble >
40TwoPiAngles_primakoff::calcAmplitude( GDouble** pKin ) const {
41
42 complex< GDouble > i( 0, 1 );
43 complex< GDouble > factor( 0, 0 );
44 complex< GDouble > Amp( 0, 0 );
45 Int_t Mrho=0;
46
47 if (flat == 1) { // no computations needed
48 Amp = 1;
49 return Amp;
50 }
51
52
53 // for Primakoff, all calculations are in the lab frame. Keep recoil but remember that it cannot be measured by detector.
54
55 TLorentzVector beam ( pKin[0][1], pKin[0][2], pKin[0][3], pKin[0][0] );
56 TLorentzVector p1 ( pKin[1][1], pKin[1][2], pKin[1][3], pKin[1][0] );
57 TLorentzVector p2 ( pKin[2][1], pKin[2][2], pKin[2][3], pKin[2][0] );
58 TLorentzVector recoil ( pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0] );
59 TLorentzVector resonance = p1 + p2;
60
61 TVector3 eps(G_COS(phipol)cos(phipol), G_SIN(phipol)sin(phipol), 0.0); // beam polarization vector in lab
62
63 TLorentzRotation resonanceBoost( -resonance.BoostVector() );
64
65 TLorentzVector beam_res = resonanceBoost * beam;
66 TLorentzVector recoil_res = resonanceBoost * recoil;
67 TLorentzVector p1_res = resonanceBoost * p1;
68
69 // choose helicity frame: z-axis opposite recoil target in rho rest frame. Note that for Primakoff recoil is defined as missing P4
70 TVector3 y = (beam.Vect().Unit().Cross(-recoil.Vect().Unit())).Unit();
71 TVector3 z = -1. * recoil_res.Vect().Unit();
72 TVector3 x = y.Cross(z).Unit();
73 TVector3 angles( (p1_res.Vect()).Dot(x),
74 (p1_res.Vect()).Dot(y),
75 (p1_res.Vect()).Dot(z) );
76
77 GDouble CosTheta = angles.CosTheta();
78 GDouble phi = angles.Phi();
79 // GDouble sinSqTheta = G_SIN(angles.Theta())*G_SIN(angles.Theta());
80 // GDouble sin2Theta = G_SIN(2.*angles.Theta());
81
82 GDouble Phi = atan2(y.Dot(eps), beam.Vect().Unit().Dot(eps.Cross(y)));
83
84 GDouble psi = Phi - phi; // define angle difference
85 if(psi < -1*PI3.141592653589793) psi += 2*PI3.141592653589793;
86 if (psi > PI3.141592653589793) psi -= 2*PI3.141592653589793;
Value stored to 'psi' is never read
87
88 /*cout << " recoil_res Angles="; recoil_res.Vect().Print();
89 cout << " p1_res Angles="; p1_res.Vect().Print();
90 cout << "Phi_pip= " << Phi_pip << endl;
91 cout << "Phi= " << Phi << endl;
92 cout << "Phi_prod= " << Phi_prod << endl;
93 cout << "phi= " << phi << endl;
94 cout << " psi=" << psi << endl;*/
95
96
97 switch (PhaseFactor) {
98 case 0:
99 Mrho = m_rho;
100 Amp = G_SQRT(1-polFrac)sqrt(1-polFrac)*(-G_SIN(Phi)sin(Phi)* Y( 0, Mrho, CosTheta, phi) );
101 break;
102 case 1:
103 Mrho = m_rho;
104 Amp = G_SQRT(1+polFrac)sqrt(1+polFrac)*(G_COS(Phi)cos(Phi)* Y( 0, Mrho, CosTheta, phi) );
105 break;
106 case 2:
107 Mrho = m_rho;
108 factor = exp(-i*Phi)* Y( 1, Mrho, CosTheta, phi);
109 Amp = G_SQRT(1-polFrac)sqrt(1-polFrac)* imag(factor);
110 break;
111 case 3:
112 Mrho = m_rho;
113 factor = exp(-i*Phi)* Y( 1, Mrho, CosTheta, phi);
114 Amp = G_SQRT(1+polFrac)sqrt(1+polFrac)* real(factor);
115 break;
116 }
117
118
119 if (abs(Amp) <= 0) cout << " m_rho=" << m_rho << " CosTheta=" << CosTheta << " phi=" << phi << " factor=" << factor << " Amp=" << Amp << endl;
120
121 return Amp;
122}
123