Bug Summary

File:libraries/AMPTOOLS_AMPS/TwoPiWt_primakoff.cc
Location:line 315, column 11
Description:Value stored to 'Thpipi' during its initialization is never read

Annotated Source Code

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
21Double_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
54Double_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
84Double_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
127Double_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
207TwoPiWt_primakoff::TwoPiWt_primakoff( const vector< string >& args ) :
208UserAmplitude< 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
230complex< GDouble >
231TwoPiWt_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
333void
334TwoPiWt_primakoff::updatePar( const AmpParameter& par ){
335
336 // could do expensive calculations here on parameter updates
337
338}
339
340