File: | libraries/AMPTOOLS_AMPS/TwoPiWt_sigma.cc |
Location: | line 148, column 14 |
Description: | Value stored to 'Eg' during its initialization is never read |
1 | |
2 | #include <cassert> |
3 | #include <iostream> |
4 | #include <string> |
5 | #include <complex> |
6 | #include <cstdlib> |
7 | |
8 | #include "TLorentzVector.h" |
9 | |
10 | #include "barrierFactor.h" |
11 | #include "breakupMomentum.h" |
12 | |
13 | #include "particleType.h" |
14 | |
15 | #include "IUAmpTools/Kinematics.h" |
16 | #include "AMPTOOLS_AMPS/TwoPiWt_sigma.h" |
17 | |
18 | // Class modeled after BreitWigner amplitude function provided for examples with AmpTools. |
19 | // Dependence of swave 2pi cross section on W (mass of 2pi system) Elton 4/17/2017 |
20 | // Version for sigma f0(500) production Elton 9/9/2018 |
21 | |
22 | complex<double> Aw_func (Double_t *x, Double_t *par) { |
23 | // Returns the amplitude for the sigma photon production. GlueX-doc-XXXX DRAFT October 1, 2018 |
24 | // "Amplitude generation of the $\sigma$ meson f$_0$(500)" |
25 | |
26 | double alpha1Re=par[0]; |
27 | double alpha1Im=par[1]; |
28 | double alpha2Re=par[2]; |
29 | double alpha2Im=par[3]; |
30 | double alpha3Re=par[4]; |
31 | double alpha3Im=par[5]; |
32 | double alpha4Re=par[6]; |
33 | double alpha4Im=par[7]; |
34 | double W=x[0]; |
35 | |
36 | complex<double> Amp1, Amp2; |
37 | complex<double> alpha1(alpha1Re,alpha1Im); |
38 | complex<double> alpha2(alpha2Re,alpha2Im); |
39 | complex<double> alpha3(alpha3Re,alpha3Im); |
40 | complex<double> alpha4(alpha4Re,alpha4Im); |
41 | |
42 | |
43 | // S-wave I=0 phase shift as a function of center of mass energy |
44 | |
45 | Double_t Mpi = ParticleMass(Pi0); |
46 | Double_t Mk = ParticleMass(KPlus); |
47 | Double_t s = W*W; // center of mass energy ^2 |
48 | Double_t s0 = (2*Mk)*(2*Mk); |
49 | |
50 | if (s <= 4*Mpi*Mpi || s >= s0) return 0; |
51 | |
52 | Double_t k = s/4 - Mpi*Mpi >= 0? sqrt(s/4 - Mpi*Mpi) : 0; // center-of-mass momentum |
53 | |
54 | // Reference Ananthanaryan, Phys Reports 353 (2001) 207, Appendix D. tan(delta) |
55 | // S-wave I=0 phase shift as a function of center of mass energy |
56 | Double_t A00=0.225; |
57 | Double_t B00=12.651; |
58 | Double_t C00=-43.8454; |
59 | Double_t D00=-87.1632; |
60 | Double_t s00=0.715311; |
61 | |
62 | Double_t k2=k*k; |
63 | Double_t f = (2*k/sqrt(s))* (A00 + B00*k2 + C00*k2*k2 + D00*k2*k2*k2)*(4*Mpi*Mpi-s00)/(s-s00); |
64 | |
65 | Double_t delta0 = atan(f); |
66 | if (delta0 < 0) delta0 += 3.14159; |
67 | |
68 | // cout << " s=" << s << " delta0=" << delta0 << endl; |
69 | |
70 | Double_t sind = sin(delta0); |
71 | Double_t cosd = cos(delta0); |
72 | complex<double> expdelta(cos(delta0),sin(delta0)); |
73 | |
74 | Amp1 = (W/(2*k)) * sind * expdelta * (alpha1 + alpha2*s); |
75 | Amp2 = cosd * expdelta* (alpha3 + alpha4*s); |
76 | complex<double> Amp = Amp1 + Amp2; |
77 | |
78 | return Amp; |
79 | |
80 | } |
81 | |
82 | |
83 | TwoPiWt_sigma::TwoPiWt_sigma( const vector< string >& args ) : |
84 | UserAmplitude< TwoPiWt_sigma >( args ) |
85 | { |
86 | |
87 | assert( args.size() == 4 )((args.size() == 4) ? static_cast<void> (0) : __assert_fail ("args.size() == 4", "libraries/AMPTOOLS_AMPS/TwoPiWt_sigma.cc" , 87, __PRETTY_FUNCTION__)); |
88 | m_par1 = AmpParameter( args[0] ); |
89 | m_par2 = AmpParameter( args[1] ); |
90 | m_daughters = pair< string, string >( args[2], args[3] ); |
91 | |
92 | // need to register any free parameters so the framework knows about them |
93 | registerParameter( m_par1 ); |
94 | registerParameter( m_par2 ); |
95 | |
96 | // make sure the input variables look reasonable |
97 | // assert( ( m_orbitL >= 0 ) && ( m_orbitL <= 4 ) ); |
98 | } |
99 | |
100 | complex< GDouble > |
101 | TwoPiWt_sigma::calcAmplitude( GDouble** pKin ) const |
102 | { |
103 | TLorentzVector P1, P2, Ptot, Ptemp, Precoil; |
104 | |
105 | for( unsigned int i = 0; i < m_daughters.first.size(); ++i ){ |
106 | |
107 | string num; num += m_daughters.first[i]; |
108 | int index = atoi(num.c_str()); |
109 | Ptemp.SetPxPyPzE( pKin[index][1], pKin[index][2], |
110 | pKin[index][3], pKin[index][0] ); // pi+ is index 1 |
111 | P1 += Ptemp; |
112 | Ptot += Ptemp; |
113 | |
114 | /* cout << " 1i=" << i << " num=" << num << " index=" << index << " P1.M=" << P1.M() << endl; |
115 | P1.Print(); |
116 | Ptot.Print();*/ |
117 | } |
118 | |
119 | for( unsigned int i = 0; i < m_daughters.second.size(); ++i ){ |
120 | |
121 | string num; num += m_daughters.second[i]; |
122 | int index = atoi(num.c_str()); |
123 | Ptemp.SetPxPyPzE( pKin[index][1], pKin[index][2], |
124 | pKin[index][3], pKin[index][0] ); // pi- is index 2 |
125 | P2 += Ptemp; |
126 | Ptot += Ptemp; |
127 | |
128 | /* cout << " 2i=" << i << " num=" << num << " index=" << index << " P2.M=" << P2.M() << endl; |
129 | P2.Print(); |
130 | Ptot.Print();*/ |
131 | } |
132 | |
133 | GDouble Wpipi = Ptot.M(); |
134 | GDouble mass1 = P1.M(); |
135 | GDouble mass2 = P2.M(); |
136 | |
137 | // get momentum transfer |
138 | Precoil.SetPxPyPzE (pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]); // Recoil is particle 3 |
139 | // next three lines commented out, unused variables |
140 | GDouble Et = Precoil.E(); |
141 | GDouble Mt = Precoil.M(); |
142 | GDouble t = -2*Precoil.M()*(Et - Mt); |
143 | |
144 | |
145 | Int_t const npar = 8; |
146 | Double_t xin[1]; |
147 | xin[0] = Wpipi; // W, 2pi mass |
148 | Double_t Eg = pKin[0][0]; // incident photon energy |
Value stored to 'Eg' during its initialization is never read | |
149 | |
150 | Double_t alpha1Re=0.378129; |
151 | Double_t alpha1Im=0; |
152 | Double_t alpha2Re=0.751557; |
153 | Double_t alpha2Im=0; |
154 | Double_t alpha3Re=0.244899; |
155 | Double_t alpha3Im=0; |
156 | Double_t alpha4Re=0.179788; |
157 | Double_t alpha4Im=0; |
158 | |
159 | Double_t parin[npar]; |
160 | parin[0] = alpha1Re; |
161 | parin[1] = alpha1Im; |
162 | parin[2] = alpha2Re; |
163 | parin[3] = alpha2Im; |
164 | parin[4] = alpha3Re; |
165 | parin[5] = alpha3Im; |
166 | parin[6] = alpha4Re; |
167 | parin[7] = alpha4Im; |
168 | |
169 | complex<GDouble> Aw = Aw_func(xin,parin); |
170 | |
171 | /*Double_t Bslope=3.7; |
172 | Double_t Bgen=6.0; |
173 | |
174 | Aw = Aw * exp(Bslope*t)/exp(6.0*t); // Divide out generated exponential. Move to separate amplitude.*/ |
175 | Aw = isfinite(real(Aw)) && isfinite(imag(Aw))? Aw : 0; // protect against infitinites |
176 | |
177 | if (Wpipi < mass1+mass2) Aw = 0; |
178 | |
179 | // cout << "TwoPiWt_sigma: calcAmplitude: 2pi mass=" << Wpipi << " Eg=" << Eg << " t=" << t << " AwNorm=" << std::norm(Aw) << " AwPhase=" << std::arg(Aw) << endl; |
180 | |
181 | return( Aw ); |
182 | } |
183 | |
184 | void |
185 | TwoPiWt_sigma::updatePar( const AmpParameter& par ){ |
186 | |
187 | // could do expensive calculations here on parameter updates |
188 | |
189 | } |
190 | |
191 |