File: | libraries/AMPTOOLS_AMPS/TwoPiWt_primakoff.cc |
Location: | line 315, column 11 |
Description: | Value stored to 'Thpipi' 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 | #include "particleType.h" |
10 | |
11 | #include "barrierFactor.h" |
12 | #include "breakupMomentum.h" |
13 | |
14 | #include "IUAmpTools/Kinematics.h" |
15 | #include "AMPTOOLS_AMPS/TwoPiWt_primakoff.h" |
16 | |
17 | // Class modeled after BreitWigner amplitude function provided for examples with AmpTools. |
18 | // Dependence of swave 2pi cross section on W (mass of 2pi system) |
19 | // Elton 4/17/2017 |
20 | |
21 | Double_t sigma_ggpi0pi0_func (Double_t *x, Double_t *par){ |
22 | |
23 | // Parameterization for data from Masiske Phys Rev D 41 (1990) 3324. |
24 | // Returns cross section in units of nb |
25 | |
26 | // constants |
27 | // Double_t const PI = 3.14159; |
28 | // parin[0] = A; // parameter 1: normalization |
29 | // parin[x] = mu; // parameter 2: mean of Fermi Function. simplify to 2 parameters |
30 | // parin[1] = kT // parameter 3: transition parameter |
31 | |
32 | Double_t MPI = ParticleMass(Pi0); |
33 | // Double_t PI=3.14159; |
34 | |
35 | Double_t A = par[0]; |
36 | Double_t mu = 2*MPI; |
37 | Double_t kT = par[1]; |
38 | Double_t Wpipi = x[0] ; |
39 | Double_t f; |
40 | |
41 | if (Wpipi < 2*MPI) { |
42 | f = 0.5; |
43 | } |
44 | else { |
45 | f = 1 / (exp((Wpipi-mu)/kT) + 1); // Use Fermi function. Data approx flat out to about 0.8 GeV, ignore beyond that. Data for costhe<0.8 |
46 | } |
47 | |
48 | |
49 | // cout << " Wpipi=" << Wpipi << " A=" << A << " mu=" << mu << " kT=" << kT << " f=" << f << endl; |
50 | |
51 | return 2*A*(0.5-f)/0.8; // Convert to nb, assuming isotropic angular distribution of pions. Correct for 80% coverage |
52 | } |
53 | |
54 | Double_t sigma_ggpipi_func (Double_t *x, Double_t *par){ |
55 | |
56 | // Parameterization from Rory for cross section for gamma gamma -> pi+ pi- |
57 | // Returns cross section in units of nb/sr |
58 | |
59 | // constants |
60 | // Double_t const PI = 3.14159; |
61 | // Double_t MPI = ParticleMass(Pi0); // pi0 |
62 | Double_t MPI = ParticleMass(PiPlus); // pi+ |
63 | Double_t W0 = 0.3; |
64 | |
65 | Double_t expon = par[0]; |
66 | // Double_t par2 = par[1]; |
67 | Double_t Wpipi = x[0] ; |
68 | Double_t f; |
69 | |
70 | if (Wpipi < 2*MPI) { |
71 | f = 0; |
72 | } |
73 | else if (Wpipi < W0) { |
74 | f = 300./(0.6*4.*PI3.141592653589793)*(Wpipi-2.*MPI)/(W0-2.*MPI); // linear rise, isotropic, CB data only 60% coverage |
75 | } |
76 | else { |
77 | f = 300./(0.6*4.*PI3.141592653589793)*pow(W0/Wpipi,expon); // power fall off, isotropic |
78 | } |
79 | |
80 | return f; |
81 | } |
82 | |
83 | |
84 | Double_t ff_func (Double_t *x, Double_t *par){ |
85 | |
86 | // return the nuclear form factor accourding to 2 parameter Fermi distribution |
87 | // See Journall of Research of the National Bureau of Standards - B Mathenatics and Mathematical Physics |
88 | // Vol. 70B, No. 1, Jan-Mar 1966. "The Form Factor of the Fermi Model Spatial Distribution," by Maximon and Schrack |
89 | // |
90 | // Function is a function of q, the three-momentum transfer to the nucleus. |
91 | // Input argument is t |
92 | // Note that q is the 3-vector momentum, but for low -t, q ~ sqrt(-t). |
93 | |
94 | // constants |
95 | // Double_t alpha = 1/137.; |
96 | // Double_t pi = 3.14159; |
97 | Double_t hbarc = 0.19733; // GeV*fm |
98 | Double_t q = sqrt(x[0])/hbarc; // take q to be in fm^-1. Input variable is positive (-t) |
99 | |
100 | Double_t R0 = par[0]; // half-density radius |
101 | Double_t a0 = par[1]; // skin or diffuseness parameter |
102 | |
103 | Double_t rho0; |
104 | Double_t sum=0; |
105 | Int_t jmax=4; |
106 | for (Int_t j=1;j<jmax;j++) { // truncate after 3 terms, first term dominates. |
107 | Double_t sum1 =pow(-1.,j-1)*exp(-j*R0/a0)/(j*j*j); |
108 | sum += sum1; |
109 | // cout << "jmax=" << jmax << " j=" << j << " R0=" << R0 << " a0=" << a0 << " sum1=" << sum1 << " sum=" << sum << endl; |
110 | } |
111 | |
112 | rho0 = (4*PI3.141592653589793*R0/3)*( PI3.141592653589793*a0*PI3.141592653589793*a0 + R0*R0 + 8*PI3.141592653589793*a0*a0*a0*sum); |
113 | rho0 = 1/rho0; |
114 | |
115 | Double_t f = 0; |
116 | |
117 | f = (4*PI3.141592653589793*PI3.141592653589793*rho0*a0*a0*a0)/(q*q*a0*a0*sinh(PI3.141592653589793*q*a0)*sinh(PI3.141592653589793*q*a0)) |
118 | * (PI3.141592653589793*q*a0 * cosh(PI3.141592653589793*q*a0) * sin(q*R0) - q*R0 *cos(q*R0) * sinh(PI3.141592653589793*q*a0) ) |
119 | + 8*PI3.141592653589793*rho0*a0*a0*a0*sum; |
120 | |
121 | // cout << " q=" << q << " f=" << f << endl; |
122 | return f; |
123 | |
124 | } |
125 | |
126 | |
127 | Double_t sigmat_func (Double_t *x, Double_t *par){ |
128 | |
129 | // return the cross section for Primakoff production of pi+pi-. CPP proposal PR12-13-008 Eq 8. |
130 | // independent variable is momentum transfer -t, in GeV2; |
131 | |
132 | // constants |
133 | Double_t alpha = 1/137.; |
134 | // Double_t pi = 3.14159; |
135 | Double_t betapipi = 0.999; // beta=0.999 (W=0.3), beta=0.997 (W=0.4), beta= 0.986 (W=1.0 GeV) |
136 | // Double_t sigmagg = 1; // take sigma (gg -> pi+ pi-) = 1 for now |
137 | // Double_t Z = 50; // Z of Sn, target |
138 | Double_t Z = 82; // Z of Pb, target |
139 | Double_t coef = 4*alpha*Z*Z/(PI3.141592653589793); |
140 | |
141 | Double_t Wpipi = par[0]; |
142 | Double_t Eg = par[1]; |
143 | Double_t t = -x[0] ; // theta of pipi system in the lab. Input Variable is positive (i.e. -t) |
144 | |
145 | Double_t xin[2]; |
146 | Double_t parin[2]; |
147 | |
148 | xin[0] = -t; // input variable to ff is positive (-t) |
149 | parin[0] = par[2]; // half-density radius, fm |
150 | parin[1] = par[3]; // diffuseness paramter, fm |
151 | |
152 | Double_t FF = ff_func(xin,parin); |
153 | |
154 | // include other masses here |
155 | |
156 | Double_t m1 = 0; // mass of photon, incident beam |
157 | // Double_t m2 = 108.*0.931494; // mass of 116Sn, target |
158 | Double_t m2 = 208.*0.931494; // use Pb mass because it is in the particle list |
159 | Double_t m3 = Wpipi; // mass of 2pi system, scattered system |
160 | Double_t m4 = m2; // recoil target |
161 | |
162 | |
163 | Double_t f = 0; |
164 | |
165 | Double_t s = m2*m2 + 2*Eg*m2; |
166 | Double_t sqrts = s > 0? sqrt(s) : 0; |
167 | if (s < 0) { |
168 | cout << "*** sigma_func: s =" << s << " < 0!" << endl; |
169 | return f; |
170 | } |
171 | |
172 | Double_t E1cm = (s+m1*m1-m2*m2)/(2*sqrts); |
173 | // Double_t E2cm = (s+m2*m2-m1*m1)/(2*sqrts); |
174 | Double_t E3cm = (s+m3*m3-m4*m4)/(2*sqrts); |
175 | // Double_t E4cm = (s+m4*m4-m3*m3)/(2*sqrts); |
176 | |
177 | Double_t p1cm = E1cm*E1cm - m1*m1? sqrt(E1cm*E1cm - m1*m1) : 0; |
178 | // Double_t p2cm = E2cm*E2cm - m2*m2? sqrt(E2cm*E2cm - m2*m2) : 0; |
179 | Double_t p3cm = E3cm*E3cm - m3*m3? sqrt(E3cm*E3cm - m3*m3) : 0; |
180 | // Double_t p4cm = E4cm*E4cm - m4*m4? sqrt(E4cm*E4cm - m4*m4) : 0; |
181 | |
182 | Double_t arg = (m1*m1-m3*m3-m2*m2+m4*m4)/(2*sqrts); |
183 | Double_t t0 = arg*arg - (p1cm - p3cm)*(p1cm - p3cm); |
184 | // Double_t t1 = arg*arg - (p1cm + p3cm)*(p1cm + p3cm); |
185 | |
186 | Double_t betastar = Eg/(Eg + m2); |
187 | Double_t gammastar = (Eg + m2)/sqrts; |
188 | Double_t betapipicm = p3cm/E3cm; |
189 | |
190 | // Double_t thepipicm = t0 -t > 0? sqrt((t0 -t)/(p1cm*p3cm)) : 0; |
191 | |
192 | Double_t conv = 1./(gammastar*(1 + betastar/betapipicm)); |
193 | |
194 | if (-t > -t0) { |
195 | f = (coef/2)* Eg*Eg*Eg*Eg * (t0-t)* betapipi*betapipi * (FF*FF/(t*t))*conv*conv*conv*conv/(p1cm*p3cm*p1cm*p3cm); |
196 | } |
197 | else { |
198 | f = 0; |
199 | } |
200 | |
201 | // if (f <= 0) cout << " t=" << t << " t0=" << t0 << " betastar=" << betastar << " gammastar=" << gammastar << " betapipicm=" |
202 | // << betapipicm << " FF=" << FF << " f=" << f << endl; |
203 | return f; |
204 | } |
205 | |
206 | |
207 | TwoPiWt_primakoff::TwoPiWt_primakoff( const vector< string >& args ) : |
208 | UserAmplitude< TwoPiWt_primakoff >( args ) |
209 | { |
210 | |
211 | assert( args.size() == 6)((args.size() == 6) ? static_cast<void> (0) : __assert_fail ("args.size() == 6", "libraries/AMPTOOLS_AMPS/TwoPiWt_primakoff.cc" , 211, __PRETTY_FUNCTION__)); |
212 | m_par1 = AmpParameter( args[0] ); |
213 | m_par2 = AmpParameter( args[1] ); |
214 | Bgen = AmpParameter( args[2] ); |
215 | mtmax = AmpParameter( args[3] ); |
216 | m_daughters = pair< string, string >( args[4], args[5] ); |
217 | |
218 | // need to register any free parameters so the framework knows about them |
219 | registerParameter( m_par1 ); |
220 | registerParameter( m_par2 ); |
221 | registerParameter( Bgen ); |
222 | registerParameter( mtmax ); |
223 | |
224 | // make sure the input variables look reasonable |
225 | // assert( ( m_orbitL >= 0 ) && ( m_orbitL <= 4 ) ); |
226 | assert( Bgen >= 1)((Bgen >= 1) ? static_cast<void> (0) : __assert_fail ("Bgen >= 1", "libraries/AMPTOOLS_AMPS/TwoPiWt_primakoff.cc" , 226, __PRETTY_FUNCTION__)); |
227 | assert( mtmax > 0)((mtmax > 0) ? static_cast<void> (0) : __assert_fail ("mtmax > 0", "libraries/AMPTOOLS_AMPS/TwoPiWt_primakoff.cc" , 227, __PRETTY_FUNCTION__)); |
228 | } |
229 | |
230 | complex< GDouble > |
231 | TwoPiWt_primakoff::calcAmplitude( GDouble** pKin ) const |
232 | { |
233 | TLorentzVector P1, P2, Ptot, Ptemp, Precoil; |
234 | |
235 | for( unsigned int i = 0; i < m_daughters.first.size(); ++i ){ |
236 | |
237 | string num; num += m_daughters.first[i]; |
238 | int index = atoi(num.c_str()); |
239 | Ptemp.SetPxPyPzE( pKin[index][1], pKin[index][2], |
240 | pKin[index][3], pKin[index][0] ); // pi+ is index 1 |
241 | P1 += Ptemp; |
242 | Ptot += Ptemp; |
243 | |
244 | /* cout << " 1i=" << i << " num=" << num << " index=" << index << " P1.M=" << P1.M() << endl; |
245 | P1.Print(); |
246 | Ptot.Print();*/ |
247 | } |
248 | |
249 | for( unsigned int i = 0; i < m_daughters.second.size(); ++i ){ |
250 | |
251 | string num; num += m_daughters.second[i]; |
252 | int index = atoi(num.c_str()); |
253 | Ptemp.SetPxPyPzE( pKin[index][1], pKin[index][2], |
254 | pKin[index][3], pKin[index][0] ); // pi- is index 2 |
255 | P2 += Ptemp; |
256 | Ptot += Ptemp; |
257 | |
258 | /* cout << " 2i=" << i << " num=" << num << " index=" << index << " P2.M=" << P2.M() << endl; |
259 | P2.Print(); |
260 | Ptot.Print();*/ |
261 | } |
262 | |
263 | GDouble Wpipi = Ptot.M(); |
264 | GDouble mass1 = P1.M(); |
265 | // GDouble mass2 = P2.M(); |
266 | GDouble Ppipi = Ptot.E() > Wpipi? sqrt(Ptot.E()*Ptot.E() - Wpipi*Wpipi): 0; |
267 | GDouble Thetapipi = Ptot.Theta()*180./PI3.141592653589793; |
268 | |
269 | // get momentum transfer |
270 | Precoil.SetPxPyPzE (pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]); // Recoil is particle 3 |
271 | GDouble Et = Precoil.E(); |
272 | GDouble Mt = Precoil.M(); |
273 | GDouble t = -2*Precoil.M()*(Et - Mt); |
274 | |
275 | // cout << "Precoil.M()=" << Precoil.M() << " T=" << Precoil.E() - Precoil.M() << " t=" << t << endl; Precoil.Print();cout << endl << endl; |
276 | |
277 | // call sigma (gamma gamma -> pi pi) cross section |
278 | |
279 | |
280 | Int_t const npar = 4; |
281 | Double_t xin[1]; |
282 | xin[0] = Wpipi; // W, 2pi mass |
283 | Double_t Eg = pKin[0][0]; // incident photon energy |
284 | Double_t parin[npar]; |
285 | // parin[0] A= 9.6; // pi0 fit to data for costhe<0.8 Marsiske Phys Rev D 41 (1990) 3324 |
286 | // parin[x] mu= 2*MPI; // pi0 should go to zero at threshold. Move to function |
287 | // parin[1] kT= 0.028; // pi0 transition over about 100 MeV. |
288 | |
289 | // parin[0] = 1.29; // charged parameter 1: exponent |
290 | // parin[1] = 0.; // charged parameter 2: par2 (spare) |
291 | parin[0] = m_par1; // parameter 1: exponent |
292 | parin[1] = m_par2; // parameter 2: par2 (spare) |
293 | // Double_t Wmin=0.2 ; |
294 | // Double_t Wmax=0.8; |
295 | |
296 | |
297 | // cout << " TwoPiWt_primakoff: m_par1=" << m_par1 << " m_par2=" << m_par2 << " Bgen=" << Bgen << endl; |
298 | |
299 | GDouble sig_ggpipi = sigma_ggpi0pi0_func(xin,parin); |
300 | |
301 | parin[0] = Wpipi; |
302 | parin[1] = Eg; |
303 | // Double_t R0 = 6.62; // Pb half-density radius, fm |
304 | // Double_t a0 = 0.546; // Pb difuseness parameter, fm |
305 | // Double_t R0 = 5.358; // Sn half-density radius, fm |
306 | // Double_t a0 = 0.550; // Sn difuseness parameter, fm |
307 | parin[2] = 6.62; |
308 | parin[3] = 0.546; |
309 | xin[0] = -t; // input positive value of t |
310 | |
311 | GDouble xnorm = 0.001; |
312 | GDouble sigmat = sigmat_func (xin,parin) * xnorm; // normlize amplitude to about unity |
313 | |
314 | GDouble tpar = (mass1*mass1/(2*Eg)) * (mass1*mass1/(2*Eg)); |
315 | GDouble Thpipi = -t > tpar? (180/PI3.141592653589793)*sqrt( (-t-tpar)/(Eg*Ppipi) ): 0; |
Value stored to 'Thpipi' during its initialization is never read | |
316 | |
317 | double epsilon = 1e-7; |
318 | complex<GDouble> RealOne(1,0); |
319 | complex<GDouble> ImagOne(0,1); |
320 | complex<GDouble> Csig; |
321 | |
322 | double arg = Bgen*t > -100? Bgen*t : -100; // limit exponential |
323 | Csig = isfinite(sigmat*sig_ggpipi/Wpipi/exp(arg))? sqrt(sigmat*sig_ggpipi/Wpipi/exp(arg)) * RealOne : 0; // Return complex double, sqrt (cross section). Divide out exponential |
324 | if (-t > mtmax) Csig = 0; // eliminate events at high t with large weights |
325 | |
326 | |
327 | // cout << "calcAmplitude: 2pi mass=" << Wpipi << " Eg=" << Eg << " t=" << t << " BGen=" << Bgen << " exp(Bgen*t)=" << exp(Bgen*t) << " m_par1=" << m_par1 << " m_par2=" |
328 | //<< m_par2 << " sig_ggpipi=" << sig_ggpipi << " sigmat=" << sigmat << " Thpipi=" << Thpipi << " Thetapipi=" << Thetapipi << " Csig=" << Csig << endl; |
329 | |
330 | return( Csig + epsilon); // return non-zero value to protect from likelihood calculation |
331 | } |
332 | |
333 | void |
334 | TwoPiWt_primakoff::updatePar( const AmpParameter& par ){ |
335 | |
336 | // could do expensive calculations here on parameter updates |
337 | |
338 | } |
339 | |
340 |