File: | libraries/AMPTOOLS_AMPS/TwoPiAngles_primakoff.cc |
Location: | line 86, column 24 |
Description: | Value stored to 'psi' is never read |
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 | |
17 | TwoPiAngles_primakoff::TwoPiAngles_primakoff( const vector< string >& args ) : |
18 | UserAmplitude< 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 | |
39 | complex< GDouble > |
40 | TwoPiAngles_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(cos(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 = sin(angles.Theta())*sin(angles.Theta()); |
80 | // GDouble sin2Theta = 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 = sqrt(1-polFrac)*(-sin(Phi)* Y( 0, Mrho, CosTheta, phi) ); |
101 | break; |
102 | case 1: |
103 | Mrho = m_rho; |
104 | Amp = sqrt(1+polFrac)*(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 = 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 = 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 |