Bug Summary

File:programs/Simulation/gen_ee/code/gen_ee.cc
Location:line 463, column 7
Description:Function call argument is an uninitialized value

Annotated Source Code

1#include <string>
2#include <iostream>
3#include "TTree.h"
4#include "TFile.h"
5#include "TLorentzVector.h"
6#include "Riostream.h"
7#include "TH1.h"
8#include "TH2.h"
9#include "TGraph.h"
10#include "TRandom3.h"
11#include "qDevilLib.h"
12#include "devilTreePT.h"
13#include "HddmOut.h"
14#include "UTILITIES/BeamProperties.h"
15#include "UTILITIES/MyReadConfig.h"
16
17//STRUCTURE TO KEEP THE CONFIGURATION SETTINGS
18struct genSettings_t {
19 int beamType; //Type of beam (1 -> single energy; 2-> bremstrahlung spectrum)
20 int runNum; //run number
21 int polDir; //beam polarization direction (0-> unpolarized; 1-> pol in x; 2->pol in y)
22 double eGammaInit; //incident photon energy
23 double eLower; //incident photon energy min (only used for beamType = 1)
24 double eUpper; //incident photon energy max (only used for beamType = 1)
25 int corrYes; //do screening and radiative corrections if corrYes = 1
26 int nToGen; //number of events to generate
27 int prescale; //number of events between printing to terminal
28 int rSeed; //seed for random number generator
29 int tOut; //type of output file (1->ROOT; 2->HDDM)
30 int reaction; //reaction (2->pair; 3->triplet)
31 char outFile[80]; //name of output file
32 //char inFileBrem[80]; //name of input root file that contains spectra histogram cobrem_vs_E
33 string beamconfigfile = "beam.cfg"; //beam cfg file
34 TString genconfigfile = "gen.cfg"; //generator cfg file
35};
36
37//FUNCTION PROTOTYPES
38double ampSqPT(int type, int polDir, TLorentzVector target, TLorentzVector beam,
39 TLorentzVector recoil,TLorentzVector q1,TLorentzVector q2);
40void printUsage(genSettings_t genSettings, int goYes);
41
42int main(int argc, char **argv){
43
44 char *argptr;
45 //SET THE DEFAULT CONFIGURATION SETTINGS
46 genSettings_t genSettings;
47 genSettings.reaction = 2;
48 genSettings.runNum = 9001;
49 genSettings.tOut = 1;
50 genSettings.polDir = 2;
51 genSettings.beamType = 1;
52 genSettings.eGammaInit = 9.0;
53 genSettings.eLower = 8.0;
54 genSettings.eUpper = 9.0;
55 genSettings.corrYes = 1;
56 genSettings.nToGen = 100000;
57 genSettings.prescale = 1000;
58 genSettings.rSeed = 103;
59 sprintf(genSettings.outFile,"genOut.root");
60 //sprintf(genSettings.inFileBrem,"cobrems.root");
61
62 char rootFile[80];
63 sprintf(rootFile,"genOut.root");
64
65 char hddmFile[80];
66 sprintf(hddmFile,"genOut.hddm");
67
68 int outFileSet = 0;
69 //COMMAND LINE PARSING
70 for (int i=1; i<argc; i++) {
1
Assuming 'i' is >= 'argc'
→
2
←
Loop condition is false. Execution continues on line 134
→
71 argptr = argv[i];
72 if (*argptr == '-') {
73 argptr++;
74 switch (*argptr) {
75 case 'c':
76 //strcpy(genSettings.beamconfigfile,++argptr);
77 genSettings.beamconfigfile = ++argptr;
78 break;
79 case 'd':
80 //strcpy(genSettings.genconfigfile,++argptr);
81 genSettings.genconfigfile = ++argptr;
82 break;
83 case 'r':
84 genSettings.rSeed = atoi(++argptr);
85 break;
86 case 'z':
87 genSettings.runNum = atoi(++argptr);
88 break;
89 case 'R':
90 genSettings.reaction = atoi(++argptr);
91 break;
92 case 'n':
93 genSettings.nToGen = atoi(++argptr);
94 break;
95 case 'p':
96 genSettings.prescale = atoi(++argptr);
97 break;
98 case 'P':
99 genSettings.polDir = atoi(++argptr);
100 break;
101 case 'b':
102 genSettings.beamType = atoi(++argptr);
103 break;
104 case 't':
105 genSettings.tOut = atoi(++argptr);
106 break;
107 case 'e':
108 genSettings.eGammaInit = atof(++argptr);
109 break;
110 case 'l':
111 genSettings.eLower = atof(++argptr);
112 break;
113 case 'u':
114 genSettings.eUpper = atof(++argptr);
115 break;
116 case 'o':
117 outFileSet = 1;
118 strcpy(genSettings.outFile,++argptr);
119 break;
120 //case 's':
121 //strcpy(genSettings.inFileBrem,++argptr);
122 //break;
123 case 'h':
124 // Do nothing right now
125 break;
126 default:
127 fprintf(stderrstderr,"Unrecognized argument: [-%s]\n\n",argptr);
128 printUsage(genSettings,0);
129 break;
130 }
131 }
132 }
133
134 for (int i=1; i<argc; i++) {
3
←
Loop condition is false. Execution continues on line 146
→
135 argptr = argv[i];
136 if (*argptr == '-') {
137 argptr++;
138 switch (*argptr) {
139 case 'h':
140 printUsage(genSettings,0);
141 break;
142 }
143 }
144 }
145
146 if (genSettings.beamconfigfile == "") {
4
←
Taking false branch
→
147 cout << "No beam configuration file: run gen_ee -h for help " << endl;
148 exit(1);
149 }
150
151 if (genSettings.genconfigfile == "") {
5
←
Taking false branch
→
152 cout << "No generator configuration file: run gen_ee -h for help " << endl;
153 exit(1);
154 }
155
156 // Get beam properties from configuration file
157 TH1D * cobrem_vs_E = 0;
158 if (genSettings.beamconfigfile != "") {
6
←
Taking true branch
→
159 BeamProperties beamProp( genSettings.beamconfigfile );
160 cobrem_vs_E = (TH1D*)beamProp.GetFlux();
161 }
162
163 TGraph * grXS_pair = new TGraph();
164 TGraph * grXS_trip = new TGraph();
165 double Z = 0, A = 0, rho = 0, Ltarget = 0;
166 // Get generator config file
167 if (genSettings.genconfigfile != "") {
7
←
Taking false branch
→
168 MyReadConfig * ReadFile = new MyReadConfig();
169 ReadFile->ReadConfigFile( genSettings.genconfigfile );
170 Double_t * m_target = ReadFile->GetConfig4Par("target");
171 Double_t * m_polDir = ReadFile->GetConfig1Par("polDir");
172 Double_t * m_corrYes = ReadFile->GetConfig1Par("corrYes");
173 Double_t * m_reaction = ReadFile->GetConfig1Par("reaction");
174 Double_t * m_beamType = ReadFile->GetConfig1Par("beamType");
175 Double_t * m_tOut = ReadFile->GetConfig1Par("tOut");
176 TString m_XS_pair = ReadFile->GetConfigName("XS_pair");
177 TString m_XS_trip = ReadFile->GetConfigName("XS_trip");
178 Z = m_target[0];
179 A = m_target[1];
180 rho = m_target[2];
181 Ltarget = m_target[3];
182 grXS_pair = new TGraph(m_XS_pair);
183 grXS_trip = new TGraph(m_XS_trip);
184 genSettings.polDir = (int) m_polDir[0];
185 genSettings.corrYes = (int) m_corrYes[0];
186 genSettings.reaction = (int) m_reaction[0];
187 genSettings.beamType = (int) m_beamType[0];
188 genSettings.tOut = (int) m_tOut[0];
189 }
190
191 if (genSettings.tOut == 1){
8
←
Taking false branch
→
192 if (outFileSet == 1) sprintf(rootFile,genSettings.outFile);
193 sprintf(hddmFile,"/dev/null");
194 }
195 if (genSettings.tOut == 2){
9
←
Taking false branch
→
196 if (outFileSet == 1) sprintf(hddmFile,genSettings.outFile);
197 sprintf(rootFile,"/dev/null");
198 }
199
200 double eGamma = genSettings.eGammaInit;
201
202 //GET THE HISTOGRAM FOR COHERENT BREMSTRAHLUNG SPECTRUM
203 TH1D* hGvsE;
204 //TFile *inCoBrem=new TFile(genSettings.inFileBrem); //Using spectrum from Richard
205 //hGvsE=(TH1D*)inCoBrem->Get("cobrem_vs_E");
206 hGvsE=(TH1D*)cobrem_vs_E;
207 TH1D* hGvsEout = (TH1D*)hGvsE->Clone("hGvsEout");
208 //hGvsEout->Reset();
209 //hGvsEout->Rebin(20);
210 int eBinLow = hGvsE->GetXaxis()->FindBin(genSettings.eLower);
211 int eBinHigh = hGvsE->GetXaxis()->FindBin(genSettings.eUpper);
212 hGvsE->GetXaxis()->SetRange(eBinLow,eBinHigh);
213 double gMax = hGvsE->GetMaximum();
214
215 //GET THE TRIPLET TO PAIR FRACTION HISTOGRAM NEEDED FOR RADIATIVE CORRECTIONS
216 //TH1D* hcsFraction;
217 //TFile *inCSfrac=new TFile("csFraction.root");
218 //hcsFraction=(TH1D*)inCSfrac->Get("hcsFraction");
219
220 //DEFINE OUTPUT FILE
221
222 HddmOut hddmGo(hddmFile);
223 int evtNumber = 0;
224
225 TFile *fout = new TFile(rootFile,"RECREATE");
226
227 // DEFINE TREE TO STORE THE DATA (SEE qDevilLib.h)
228 TTree *t1 = new TTree("t1","genDevilPairs");
229 devilTreePT_t devilTree;
230 setBranchesT1(t1, &devilTree);
231
232 //PRINT OUT THE SETTINGS
233 printUsage(genSettings,1);
234
235 double PIval =2*atan2(1,0);
236 double alphaQED = 1.0/137.036;
237 double hbarcSqr = 389.37966;
238 double crossSection = 0.0;
239 double fullWeight;
240 double mElectron = 0.51099907e-3;
241 double mProton = 0.938;
242 //SETTING THE CUT PARAMETERS IS A
243 //BALANCING ACT BETWEEN SPEED OF CONVERGENCE
244 //AND MAKING SURE THAT THE PHASE SPACE IS COMPLETE
245 double Mcut= 5.0e-3;
246 if (genSettings.reaction == 3) Mcut= 20.0e-3;
10
←
Taking false branch
→
247 if (genSettings.reaction == 2) Mcut= 1.0;
11
←
Taking false branch
→
248 double qRcut = 1.0e-3;
249 if (genSettings.reaction == 2) qRcut= 2.0;
12
←
Taking false branch
→
250
251 //DEFINE SOME FOUR VECTORS
252 TLorentzVector beam(0.0,0.0,eGamma,eGamma);
253 TLorentzVector target;
254 if (genSettings.reaction == 2) target.SetPxPyPzE(0.0,0.0,0.0,mProton);
13
←
Taking false branch
→
255 if (genSettings.reaction == 3) target.SetPxPyPzE(0.0,0.0,0.0,mElectron);
14
←
Taking false branch
→
256 TLorentzVector wVec = beam + target;
257 TLorentzVector q1; //electron from pair
258 TLorentzVector q2; //positron from pair
259 TLorentzVector recoil;
260 TLorentzVector q12;
261 TLorentzVector q23;
262 TLorentzVector moTransfer;
263
264 if (genSettings.reaction == 2) recoil.SetPxPyPzE(0.0,0.0,0.0,mProton);
15
←
Taking false branch
→
265 if (genSettings.reaction == 3) recoil.SetPxPyPzE(0.0,0.0,0.0,mElectron);
16
←
Taking false branch
→
266
267 //DEFINE SOME HISTGRAMS
268 double wMax = sqrt(target.Mag2() + 2*eGamma*target.Mag());
269 if (genSettings.beamType == 2) {
17
←
Taking false branch
→
270 wMax = sqrt(target.Mag2() + 2*genSettings.eUpper*target.Mag());
271 }
272 double m12MinSq = pow(q1.Mag() + q2.Mag(),2);
273 double m23MinSq = pow(q2.Mag() + recoil.Mag(),2);
274 double m12MaxSq = pow(wMax - recoil.Mag(),2);
275 double m23MaxSq = pow(wMax - q1.Mag(),2);
276 int nBinXY = sqrt(genSettings.nToGen/100);
277 TH2D* hPhaseSpaceR = new TH2D("hPhaseSpaceR","",nBinXY,m12MinSq,m12MaxSq,nBinXY,m23MinSq,m23MaxSq);
278
279
280 TRandom3 *random = new TRandom3;
281 random->SetSeed(genSettings.rSeed);
282 double sigmaVal;
18
←
'sigmaVal' declared without an initial value
→
283 PhasePT phaseGenR; //phase space object (based off Richards calculations)
284 phaseGenR.SetRCut(qRcut);
285 phaseGenR.SetM12Cut(Mcut);
286 phaseGenR.SetBeam(beam);
287 phaseGenR.SetTarget(target);
288 double phaseSpaceWeight = 0.0;
289 double sum=0.0;
290 double sum2=0.0;
291
292 double sumTest1=0.0;
293 double sumTest2=0.0;
294 double sumTest3=0.0;
295 double sumTest4=0.0;
296 double sumTest5=0.0;
297
298 int genVal;
299 int nTest = 0;
300 double yMax,yVal,testValY;
301 int eBin;
302 int nGen = 0;
303 int nSkip = 0;
304 //nSkip = 79992;//ASDF
305 for (int nGenTmp = 1; nGenTmp <= genSettings.nToGen; nGenTmp++) {
19
←
Loop condition is true. Entering loop body
→
306 genVal = -1;
307 //MONTE CARLO THE COHERENT BREMSTRAHLUNG SPECTRUM
308 if (genSettings.beamType == 2){
20
←
Taking false branch
→
309 yMax = gMax*1.02;
310 yVal = 0.0;
311 testValY = yMax + 10.0;
312 while(testValY>yVal){//Monte Carlo the event to get brem spectrum
313 eGamma = random->Uniform(genSettings.eLower,genSettings.eUpper); //Grab a photon energy
314 testValY = random->Uniform(0.0,yMax); //Grab a test value
315 eBin = hGvsE->GetXaxis()->FindBin(eGamma);
316 yVal = hGvsE->GetBinContent(eBin);
317 }
318 hGvsEout->Fill(eGamma);
319 beam.SetPxPyPzE(0,0,eGamma,eGamma);
320 wVec = beam + target;
321 phaseGenR.SetBeam(beam);
322 }
323 while (genVal < 0) {
21
←
Loop condition is true. Entering loop body
→
28
←
Assuming 'genVal' is >= 0
→
29
←
Loop condition is false. Execution continues on line 334
→
324 if (nGenTmp >= nSkip) nTest++;
22
←
Taking true branch
→
325 //GENERATE THE PHASE SPACE EVENT
326 genVal = phaseGenR.Gen(random);
327 if (genVal == -1) sumTest1 += 1.0;
23
←
Taking false branch
→
328 if (genVal == -2) sumTest2 += 1.0;
24
←
Taking false branch
→
329 if (genVal == -3) sumTest3 += 1.0;
25
←
Taking false branch
→
330 if (genVal == -4) sumTest4 += 1.0;
26
←
Taking false branch
→
331 if (genVal == -5) sumTest5 += 1.0;
27
←
Taking false branch
→
332 //NOTE: if (genVal < 0) then phase space not physical
333 }
334 if (nGenTmp < nSkip) continue;
30
←
Taking false branch
→
335 nGen++;
336
337 //GET PHASE SPACE WEIGHT
338 phaseSpaceWeight = phaseGenR.GetWeight();
339
340 //GET THE FOUR-VECTORS
341 beam = phaseGenR.GetBeam();
342 target = phaseGenR.GetTarget();
343 q1 = phaseGenR.GetQ1();
344 q2 = phaseGenR.GetQ2();
345 q12 = phaseGenR.GetQ12();
346 q23 = phaseGenR.GetQ23();
347 recoil = phaseGenR.GetRecoil();
348 moTransfer = recoil - target;
349
350 //FILL PHASE SPACE HISTOGRAM
351 hPhaseSpaceR->Fill(q12.Mag2(),q23.Mag2(),phaseSpaceWeight);
352
353 //NOTE: Richard's calculation for the phase space gives small recoil momentum
354 //very high probability. To account for this, the phase space weight is large
355 //when the recoil is large. The problem is that these rare events have a large
356 //weight. The rare events should not be a problem, except that one needs to
357 //generate a huge number of events to get the proper value. Instead, in the case
358 //of triplet production, we can notice that every diagram has a corresponding
359 //diagram that has the recoil electon momentum switched with the momentum of the
360 //produced electron. This means that the cross section for the phase space where
361 // q1.P > recoil.P will give identical results to the cross section for the
362 //phase space where q1.P < recoil.P. Because of this symmetry between
363 //q1.P and recoil.P, we can calculate the cross sections for the case
364 //q1.P > recoil.P and just make sure that we account for the "lost" phase space.
365 //Since we already have to account for the unphysical parts of the phase space
366 //generation, accounting for the case where a1.P < recoil.P does not require
367 //any additional work at this point in the code. For pair production, luckily
368 //the rare events (high momentum protons) are so rare that I have not seen
369 //any instability in the cross section results due to the rare large-weight
370 //events.
371 if (genSettings.reaction == 3 && recoil.P() > q1.P()) continue;
372 //If you comment out the line above and generate triplets, you should
373 //notice that everything is fine with the cross sections except that
374 //at some point (usualy after a very large number of events) there is
375 //a spike in the cross section.
376
377 //ADDITIONAL FACTORS TO GET CROSS SECTION (micro barns)
378 double fluxFactor = 4*beam.E()*(target.P() + target.E());
379 double rhoFactor = 1.0/(8*recoil.E()*q12.P());
380 double piFactor = pow(2*PIval,4-9)*pow(4*PIval,3);
381
382 //SCREENING AND RADIATIVE CORRECTIONS
383 double bohrRadius = (1.0/alphaQED)*(1.0/mElectron);
384 double fH = 1.0/pow(1 + pow(bohrRadius*moTransfer.P()/2.0,2),2);
385
386 double sHfactor = 1.0;
387 double sHfactorPair = pow(1.0 - fH,2); //Screening for pair production
388 double sHfactorTrip = 1.0 - pow(fH,2); //Screening for triplet production
389
390 //RADIATIVE CORRECTIONS FOR PAIRS IS SIMPLY A COMMON FACTOR :)
391 double radFactor = 1.0;
392 double radFactorDelta = 0.0093;
393 double radFactorPair = 1.0+radFactorDelta;
394 //THE RADIATIVE CORRECTION FOR TRIPLETS IS THE SAME MAG AS FOR PAIRS
395 //THIS MEANS WE HAVE TO GET THE FRACTION OF TRIPLETS TO PAIRS AND
396 //SCALE THE radFactorDelta USING THE TRIPLET TO PAIR FRACTION
397 //int eFracBin = hcsFraction->GetXaxis()->FindBin(beam.E());
398 double xs_p = grXS_pair->Eval(beam.E() * 1e3);
399 double xs_t = grXS_trip->Eval(beam.E() * 1e3);
400 if (beam.E() > 100.0) {
31
←
Taking false branch
→
401 xs_p = grXS_pair->Eval(100.0 * 1e3);
402 xs_t = grXS_trip->Eval(100.0 * 1e3);
403 }
404 double csFrac = xs_t / xs_p;//hcsFraction->GetBinContent(eFracBin);
405 double radFactorTrip = 1.0 + radFactorDelta/csFrac;
406
407 //GET THE CROSS SECTION
408 if (genSettings.reaction == 2) { //PAIR PRODUCTION
32
←
Taking false branch
→
409 radFactor = radFactorPair;
410 sHfactor = sHfactorPair;
411 crossSection = ampSqPT(2,genSettings.polDir,target,beam,recoil,q1,q2)*hbarcSqr*pow(alphaQED,3)
412 / fluxFactor * rhoFactor * piFactor;
413 }
414 if (genSettings.reaction == 3) { //TRIPLET PRODUCTION
33
←
Taking false branch
→
415 radFactor = radFactorTrip;
416 sHfactor = sHfactorTrip;
417 crossSection = ampSqPT(3,genSettings.polDir,target,beam,recoil,q1,q2)*hbarcSqr*pow(alphaQED,3)
418 / fluxFactor * rhoFactor * piFactor;
419 }
420
421 if (genSettings.corrYes == 1) {
34
←
Taking false branch
→
422 crossSection *= sHfactor;
423 crossSection *= radFactor;
424 if (genSettings.reaction == 2) crossSection *= pow(Z, 2);
425 if (genSettings.reaction == 3) crossSection *= Z;
426 }
427
428 //GET INITIAL fullWeight
429 fullWeight = crossSection*phaseSpaceWeight;
430
431 //CHECK FOR SOMETHING BAD
432 if (fullWeight >= 0 || fullWeight <=0){
433 //do nothing
434 } else {
435 cout<<"!!!!Something Bad!!!!"<<endl;
436 }
437 //CORRECT fullWeight FOR UNPHYSICAL PHASE SPACE EVENTS THAT ARE REMOVED
438 Double_t phaseSpaceCorrection = nGen*1.0/(1.0*nTest);
439 fullWeight *= phaseSpaceCorrection;
440
441 //FILL THE TREE
442 devilTree.eGamma = eGamma;
443 devilTree.weight = fullWeight;
444 devilTree.recoilE = recoil.E();
445 devilTree.recoilPx = recoil.Px();
446 devilTree.recoilPy = recoil.Py();
447 devilTree.recoilPz = recoil.Pz();
448 devilTree.electronE = q1.E();
449 devilTree.electronPx = q1.Px();
450 devilTree.electronPy = q1.Py();
451 devilTree.electronPz = q1.Pz();
452 devilTree.positronE = q2.E();
453 devilTree.positronPx = q2.Px();
454 devilTree.positronPy = q2.Py();
455 devilTree.positronPz = q2.Pz();
456 t1->Fill();
457
458 //CALCULATING CROSS SECTION ERROR (SAME WAY THAT RICHARD DOES) AND PRINT RESULT
459 sum += fullWeight;
460 sum2 += pow(fullWeight,2);
461
462 if (nGen/genSettings.prescale*genSettings.prescale == nGen) {
35
←
Taking true branch
→
463 cout <<"Integrated cross section after " << nGen << " events : "
36
←
Function call argument is an uninitialized value
464 << sum/nGen << " +/- " <<sqrt(sum2-pow(sum,2)/nGen)/nGen<< sigmaVal <<" micro barns"<< endl;
465 }
466 //HDDM STUFF
467 tmpEvt_t tmpEvt;
468 tmpEvt.beam = beam;
469 tmpEvt.target = target;
470 tmpEvt.q1 = q1;
471 tmpEvt.q2 = q2;
472 tmpEvt.recoil = recoil;
473 tmpEvt.nGen = 3;
474 tmpEvt.rxn = genSettings.reaction;
475 tmpEvt.weight = fullWeight;
476 evtNumber++;
477 if (genSettings.tOut == 2) hddmGo.write(tmpEvt,genSettings.runNum,evtNumber);
478 }
479 //WRITE THE TREE AND HISTOGRAMS TO FILE
480 if (genSettings.tOut == 1){
481 fout->cd();
482 t1->Write();
483 hPhaseSpaceR->Write();
484 hGvsEout->Write();
485 }
486 fout->Close();
487 cout<<"All done. Bye"<<endl;
488}
489
490void printUsage(genSettings_t genSettings, int goYes){
491 if (goYes == 0){
492 fprintf(stderrstderr,"\nSWITCHES:\n");
493 fprintf(stderrstderr,"-h\tPrint this message\n");
494 fprintf(stderrstderr,"-z<arg>\trun number\n");
495 fprintf(stderrstderr,"-c<arg>\tbeam configuration file\n");
496 fprintf(stderrstderr,"-d<arg>\tgenerator configuration file\n");
497 fprintf(stderrstderr,"-R<arg>\tReaction:\n");
498 fprintf(stderrstderr,"\t\t-R2 = Pair production off of proton\n");
499 fprintf(stderrstderr,"\t\t-R3 = Triplet production\n");
500 fprintf(stderrstderr,"-n<arg>\tNumber of events to generate\n");
501 fprintf(stderrstderr,"-r<arg>\tUser defined random number seed\n");
502
503 fprintf(stderrstderr,"-t<arg>\tType of output file\n");
504 fprintf(stderrstderr,"\t\t-t1 = ROOT file\n");
505 fprintf(stderrstderr,"\t\t-t2 = HDDM file\n");
506
507 fprintf(stderrstderr,"-p<arg>\tPrescale factor (number of events between printing to terminal)\n");
508 fprintf(stderrstderr,"-P<arg>\tPhoton beam polarization direction:\n");
509 fprintf(stderrstderr,"\t\t-P0 = Unpolarized\n");
510 fprintf(stderrstderr,"\t\t-P2 = Polarized in x-direction (100 percent)\n");
511 fprintf(stderrstderr,"\t\t-P3 = Polarized in y-direction (100 percent)\n");
512 fprintf(stderrstderr,"-b<arg>\tBeam type:\n");
513 fprintf(stderrstderr,"\t\t-b1 = Single photon energy\n");
514 fprintf(stderrstderr,"\t\t-b2 = Bremstrahlung spectra\n");
515 fprintf(stderrstderr,"-e<arg>\tPhoton energy in GeV. ONLY USED IF -b1\n");
516 fprintf(stderrstderr,"-l<arg>\tMinimum incident photon energy in GeV. ONLY USED IF -b2\n");
517 fprintf(stderrstderr,"-u<arg>\tMaximum incident photon energy in GeV. ONLY USED IF -b2\n");
518 //fprintf(stderr,"-s<arg>\tfile with histogram cobrem_vs_E. ONLY USED IF -b2\n");
519 fprintf(stderrstderr,"-o<arg>\tOutFile name\n");
520
521 cout<<""<<endl;
522 cout<<"The above switches overide the default setting."<<endl;
523 cout<<""<<endl;
524
525 if (genSettings.beamType == 1){
526 cout<<"The current operation is equivalent to the command:"<<endl;
527 cout<<"genDevilPT"
528 <<" -z"<<genSettings.runNum
529 <<" -c"<<genSettings.beamconfigfile
530 <<" -d"<<genSettings.genconfigfile
531 <<" -n"<<genSettings.nToGen
532 <<" -R"<<genSettings.reaction
533#ifdef WITHHDDM
534 <<" -t"<<genSettings.tOut
535#endif
536 <<" -r"<<genSettings.rSeed
537 <<" -p"<<genSettings.prescale
538 <<" -b"<<genSettings.beamType
539 <<" -P"<<genSettings.polDir
540 <<" -e"<<genSettings.eGammaInit
541 //<<" -l"<<genSettings.eLower
542 //<<" -u"<<genSettings.eUpper
543 <<" -o"<<genSettings.outFile
544 //<<" -s"<<genSettings.inFileBrem
545 <<"\n"<<endl;
546 } else {
547 cout<<"The current operation is equivalent to the command:"<<endl;
548 cout<<"genDevilPT"
549 <<" -z"<<genSettings.runNum
550 <<" -c"<<genSettings.beamconfigfile
551 <<" -d"<<genSettings.genconfigfile
552 <<" -n"<<genSettings.nToGen
553 <<" -R"<<genSettings.reaction
554 <<" -t"<<genSettings.tOut
555 <<" -r"<<genSettings.rSeed
556 <<" -p"<<genSettings.prescale
557 <<" -b"<<genSettings.beamType
558 <<" -P"<<genSettings.polDir
559 //<<" -e"<<genSettings.eGammaInit
560 <<" -l"<<genSettings.eLower
561 <<" -u"<<genSettings.eUpper
562 <<" -o"<<genSettings.outFile
563 //<<" -s"<<genSettings.inFileBrem
564 <<"\n"<<endl;
565 }
566 cout<<""<<endl;
567 cout<<"NOTE 1-> If providing custom photon spectrum histogram:"<<endl;
568 cout<<" * The file containing the histogram is provided"<<endl;
569 //cout<<" as an argument to the -s switch"<<endl;
570 cout<<" * The histogram MUST be of type TH1D"<<endl;
571 cout<<" * The histogram MUST have energy units of GeV for x-axis"<<endl;
572 cout<<" * The histogram MUST be named cobrem_vs_E"<<endl;
573 cout<<" * The energy bin widths in the histogram are not important"<<endl;
574 cout<<""<<endl;
575//
576// cout<<"NOTE 2-> To be able to output in HDDM mode you must compile like:"<<endl;
577// cout<<" make HDDM=yes"<<endl;
578// cout<<""<<endl;
579
580 }
581
582 if (goYes == 1){
583 if(genSettings.reaction == 2) cout<<"\nReaction: Pair production off of proton"<<endl;
584 if(genSettings.reaction == 3) cout<<"\nReaction: Triplet production"<<endl;
585 cout<<"Generating "<<genSettings.nToGen<<" events"
586 <<" with output file "<<genSettings.outFile<<endl;
587 if(genSettings.beamType == 1){
588 cout<<"Beam type: Single photon energy = "<<genSettings.eGammaInit<<" GeV"<<endl;
589 }
590 if(genSettings.beamType == 2){
591 cout<<"Beam type: Bremstrahlung spectrum from histogram "<<"cobrem_vs_E"<<endl;
592 cout<<" within file "<<genSettings.beamconfigfile<<endl;
593 cout<<" Using energy range from "
594 <<genSettings.eLower<<" GeV to "<<genSettings.eUpper<<" GeV"<<endl;
595 }
596 cout<<"Random number seed = "<<genSettings.rSeed<<endl;
597 cout<<""<<endl;
598 }
599
600 if (goYes == 0) exit(0);
601}