Bug Summary

File:libraries/AMPTOOLS_AMPS/TwoPiWt_sigma.cc
Location:line 142, column 11
Description:Value stored to 't' 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
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
22complex<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
83TwoPiWt_sigma::TwoPiWt_sigma( const vector< string >& args ) :
84UserAmplitude< 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
100complex< GDouble >
101TwoPiWt_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);
Value stored to 't' during its initialization is never read
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
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
184void
185TwoPiWt_sigma::updatePar( const AmpParameter& par ){
186
187 // could do expensive calculations here on parameter updates
188
189}
190
191