Bug Summary

File:programs/Simulation/genpi/genpi+pi-.cc
Location:line 149, column 4
Description:Dereference of null pointer

Annotated Source Code

1
2
3#include <iostream>
4#include <fstream>
5#include <string>
6#include <cmath>
7#include <stdlib.h>
8#include <time.h>
9#include <vector>
10using namespace std;
11
12#include <TVector3.h>
13#include <TLorentzVector.h>
14
15double PI_CHARGED_MASS = 0.139568;
16double PROTON_MASS = 0.938;
17double MESON_MASS = 0.770;
18unsigned int MAX_EVENTS=10000;
19int NUM_TO_GEN=2;
20double E_BEAM_MIN=4.0*PI_CHARGED_MASS;
21double E_BEAM_MAX=1.0;
22int RUN_NUMBER=100;
23string OUTPUT_FILENAME="genpiX.ascii";
24double gINTEGRAL_FRACTION;
25
26#define GAMMA_TYPE1 1
27#define PI_PLUS_TYPE8 8
28#define PI_MINUS_TYPE9 9
29
30#define _DBG_cout<<"genpi+pi-.cc"<<":"<<30<<" " cout<<__FILE__"genpi+pi-.cc"<<":"<<__LINE__30<<" "
31#define _DBG__cout<<"genpi+pi-.cc"<<":"<<31<<" "<<
endl _DBG_cout<<"genpi+pi-.cc"<<":"<<31<<" "<<endl
32
33class piX{
34 public:
35 double px,py,pz,E; // piX
36};
37
38
39extern "C"{
40double rtsafe(void (*funcd)(double, double *, double *), double x1, double x2, double xacc);
41double rtnewt(void (*funcd)(double, double *, double *), double x1, double x2, double xacc);
42double zbrent(double (*func)(double), double x1, double x2, double tol);
43void funcd(double, double *f, double *df);
44double func(double);
45}
46
47double SampleSin2Theta(void);
48void ParseCommandLineArguments(int narg, char* argv[]);
49void Usage(void);
50
51
52//----------------------------
53// main
54//----------------------------
55int main(int narg, char* argv[])
56{
57
58 // Parse the command line
59 ParseCommandLineArguments(narg, argv);
60
61 // Open file for output
62 ofstream of(OUTPUT_FILENAME.c_str());
63 if(!of.is_open()){
1
Taking false branch
→
64 cout<<"Unable to open \""<<OUTPUT_FILENAME<<"\" for writing!"<<endl;
65 exit(-1);
66 }
67
68 // Seed random number generator
69 srandom(time(NULL__null));
70
71
72 // Loop over events
73 unsigned int nevents;
74 for(nevents=1; nevents<=MAX_EVENTS; nevents++){
2
←
Loop condition is true. Entering loop body
→
6
←
Loop condition is true. Entering loop body
→
10
←
Loop condition is true. Entering loop body
→
14
←
Loop condition is true. Entering loop body
→
75
76 // Determine incident photon energy. (This is just
77 // evenly sampled from the specified range.)
78 double Eg = (double)random()/(double)RAND_MAX2147483647*(E_BEAM_MAX-E_BEAM_MIN) + E_BEAM_MIN;
79
80
81#if 0
82 // tmin occurs when the target proton goes out perpendicular to
83 // the beam. This corresponds to the smallest theta angle the
84 // rho can go to. Calculate th rho angle and momentum here.
85 double rho_mom = 2.0*pow(Eg, 2.0);
86 rho_mom -= pow(PROTON_MASS,2.0) - pow(MESON_MASS, 2.0);
87 rho_mom /= 2.0*Eg;
88 rho_mom = sqrt(pow(rho_mom, 2.0) - pow(MESON_MASS, 2.0));
89#endif
90
91 // Calculate rho theta and phi angles
92 double rho_theta = 0.0;
93 double rho_phi = 2.0*M_PI3.14159265358979323846*((double)random()/(double)RAND_MAX2147483647);
94#if 1
95 // Calculate rho momentum based on angle
96 double M_N = 938.0;
97 double L = Eg/(Eg + M_N);
98 double K = (Eg*M_N + pow(MESON_MASS,2.0)/2.0)/(Eg + M_N);
99 double A = pow(L*cos(rho_theta),2.0) - 1.0;
100 double B = 2.0*K*L*cos(rho_theta);
101 double C = pow(K,2.0) - pow(MESON_MASS,2.0);
102 double rho_mom = ((-B) - sqrt(B*B - 4.0*A*C))/(2.0*A);
103#endif
104 TVector3 rho_p;
105 rho_p.SetMagThetaPhi(rho_mom, rho_theta, rho_phi);
106
107 // Find decay angle of pions in rest frame of rho
108 double pi1_theta = SampleSin2Theta();
109 double pi1_phi = 2.0*M_PI3.14159265358979323846*((double)random()/(double)RAND_MAX2147483647);
110 double pi2_theta = M_PI3.14159265358979323846 - pi1_theta;
111 double pi2_phi = pi1_phi + M_PI3.14159265358979323846;
112 if(pi2_phi>2.0*M_PI3.14159265358979323846)pi2_phi-=2.0*M_PI3.14159265358979323846;
3
←
Taking false branch
→
7
←
Taking false branch
→
11
←
Taking false branch
→
15
←
Taking false branch
→
113 double pi_E = MESON_MASS/2.0;
114 double pi_mom = sqrt(pow(pi_E,2.0) - pow(PI_CHARGED_MASS,2.0));
115 TLorentzVector pi1( pi_mom*sin(pi1_theta)*cos(pi1_phi)
116 , pi_mom*sin(pi1_theta)*sin(pi1_phi)
117 , pi_mom*cos(pi1_theta)
118 , pi_E);
119 TLorentzVector pi2( pi_mom*sin(pi2_theta)*cos(pi2_phi)
120 , pi_mom*sin(pi2_theta)*sin(pi2_phi)
121 , pi_mom*cos(pi2_theta)
122 , pi_E);
123
124 // Boost the pions into the lab frame
125 TVector3 beta = (1.0/sqrt(rho_p.Mag2() + pow(MESON_MASS,2.0)))*rho_p;
126 pi1.Boost(beta);
127 pi2.Boost(beta);
128
129 // Generate piXs
130 vector<piX> piXs;
131
132 piX p;
133 p.E = pi1.E();
134 p.px = pi1.Px();
135 p.py = pi1.Py();
136 p.pz = pi1.Pz();
137 piXs.push_back(p);
138
139 p.E = pi2.E();
140 p.px = pi2.Px();
141 p.py = pi2.Py();
142 p.pz = pi2.Pz();
143 piXs.push_back(p);
144
145 // Write event to file
146 unsigned int type = PI_PLUS_TYPE8;
147 of<<RUN_NUMBER<<" "<<nevents<<" "<<piXs.size()<<endl;
148 for(unsigned int j=0; j<piXs.size(); j++){
4
←
Loop condition is false. Execution continues on line 160
→
8
←
Loop condition is false. Execution continues on line 160
→
12
←
Loop condition is false. Execution continues on line 160
→
16
←
Loop condition is true. Entering loop body
→
149 piX &p = piXs[j];
17
←
Dereference of null pointer
150 unsigned int index = j+1;
151
152 of<<index<<" "<<type<<" "<<0<<endl;
153 of<<" "<<0<<" "<<p.px<<" "<<p.py<<" "<<p.pz<<" "<<p.E<<endl;
154
155 // alternate bewtween pi+ and pi-
156 type = type==PI_PLUS_TYPE8 ? PI_MINUS_TYPE9:PI_PLUS_TYPE8;
157 }
158
159 // Update screen
160 if(nevents%1000 == 0){
5
←
Taking false branch
→
9
←
Taking false branch
→
13
←
Taking false branch
→
161 cout<<" "<<nevents<<" events generated"<<endl;
162 }
163
164 }
165
166 of.close();
167
168 cout<<endl<<""<<nevents-1<<" events written to "<<OUTPUT_FILENAME<<endl;
169
170 return 0;
171}
172
173
174//----------------------------
175// SampleSin2Theta
176//----------------------------
177double SampleSin2Theta(void)
178{
179 // Randomly sample a value of theta between 0 and pi
180 // such than the distributions goes like sin^2.
181
182 // Get random number evenly sampled from 0 to 1
183 gINTEGRAL_FRACTION = (double)random()/(double)RAND_MAX2147483647;
184
185 // The integral fraction of sin^2 as a function
186 // of theta over the range zero pi is given by
187 //
188 // f = (theta - sin(2*theta))/pi
189 //
190 // Unfortunately, this is a trancendental equation
191 // so we must find the root numerically.
192 // We use the numerical recipes routine rtsafe
193 // and define the function and derivative in
194 // funcd.
195 return zbrent(func, 0.0, M_PI3.14159265358979323846, 1.0E-5);
196}
197
198//----------------------------
199// funcd
200//----------------------------
201void funcd(double theta, double *f, double *df)
202{
203 *f = gINTEGRAL_FRACTION - (theta - sin(2.0*theta))/M_PI3.14159265358979323846;
204 *df = -(1.0 - 2.0*cos(2.0*theta))/M_PI3.14159265358979323846;
205}
206
207//----------------------------
208// func
209//----------------------------
210double func(double theta)
211{
212 return gINTEGRAL_FRACTION - (theta - 0.5*sin(2.0*theta))/M_PI3.14159265358979323846;
213}
214
215//----------------------------
216// ParseCommandLineArguments
217//----------------------------
218void ParseCommandLineArguments(int narg, char* argv[])
219{
220 if(narg==1)Usage();
221
222 for(int i=1; i<narg; i++){
223 string arg = argv[i];
224 if(arg=="-h" || arg=="--help"){
225 Usage();
226 }else if(arg=="-M"){
227 if(i==narg-1){cout<<"-M requires an argument!"<<endl; Usage();}
228 MAX_EVENTS = atoi(argv[++i]);
229 }else if(arg=="-N"){
230 if(i==narg-1){cout<<"-N requires an argument!"<<endl; Usage();}
231 NUM_TO_GEN = atoi(argv[++i]);
232 }else if(arg=="-Emin"){
233 if(i==narg-1){cout<<"-Emin requires an argument!"<<endl; Usage();}
234 E_BEAM_MIN = atof(argv[++i]);
235 }else if(arg=="-Emax"){
236 if(i==narg-1){cout<<"-Emax requires an argument!"<<endl; Usage();}
237 E_BEAM_MAX = atof(argv[++i]);
238 }else if(arg=="-o"){
239 if(i==narg-1){cout<<"-o requires an argument!"<<endl; Usage();}
240 OUTPUT_FILENAME = argv[++i];
241 }else if(arg=="-m"){
242 if(i==narg-1){cout<<"-m requires an argument!"<<endl; Usage();}
243 PI_CHARGED_MASS = atof(argv[++i]);
244 }
245 }
246
247 cout<<"---- genpiX will use the following settings: ----"<<endl;
248 cout<<"MAX_EVENTS = "<<MAX_EVENTS<<endl;
249 cout<<"NUM_TO_GEN = "<<NUM_TO_GEN<<endl;
250 cout<<"E_BEAM_MIN = "<<E_BEAM_MIN<<endl;
251 cout<<"E_BEAM_MAX = "<<E_BEAM_MAX<<endl;
252 cout<<"PI_CHARGED_MASS = "<<PI_CHARGED_MASS<<endl;
253 cout<<"OUTPUT_FILENAME = \""<<OUTPUT_FILENAME<<"\""<<endl;
254 cout<<"-------------------------------------------------"<<endl;
255
256}
257
258//----------------------------
259// Usage
260//----------------------------
261void Usage(void)
262{
263 cout<<endl;
264 cout<<"Usage:"<<endl;
265 cout<<" genpiX -M numEvents [options]"<<endl;
266 cout<<endl;
267 cout<<" options:"<<endl;
268 cout<<" -h print this help message"<<endl;
269 cout<<" -M numEvents set the number of events to produced (required)"<<endl;
270 cout<<" -N numpiXs set the number of piXs to generate per event."<<endl;
271 cout<<" -Emin E set the lower beam energy limit in GeV"<<endl;
272 cout<<" -Emax E set the upper beam energy limit in GeV"<<endl;
273 cout<<" -o filename set the output filename"<<endl;
274 cout<<" -m mass set the rest mass of the pi(GeV) (e.g. make it an eta!)"<<endl;
275 cout<<endl;
276 cout<<"This program will produce events with one or more piXs and write"<<endl;
277 cout<<"out the resulting decay photons in an ASCII file of the same"<<endl;
278 cout<<"format as produced by genr8. There is pretty much no real physics"<<endl;
279 cout<<"here other than the piX decays isotropically in its rest frame."<<endl;
280 cout<<"The total beam energy is evenly sampled from the given range and"<<endl;
281 cout<<"The piXs divide that up in a semi-random way so that each should"<<endl;
282 cout<<"get a reasonable amount of kinetic energy. The total momentum is"<<endl;
283 cout<<"NOT conserved however as each pion is given a random angle"<<endl;
284 cout<<"distributed isotropically."<<endl;
285 cout<<endl;
286 cout<<"This is intended for testing the calorimeter (BCAL and FCAL)"<<endl;
287 cout<<"reconstruction code."<<endl;
288
289 cout<<endl;
290 exit(0);
291}
292
293