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 | |
18 | struct genSettings_t { |
19 | int beamType; |
20 | int runNum; |
21 | int polDir; |
22 | double eGammaInit; |
23 | double eLower; |
24 | double eUpper; |
25 | int corrYes; |
26 | int nToGen; |
27 | int prescale; |
28 | int rSeed; |
29 | int tOut; |
30 | int reaction; |
31 | char outFile[80]; |
32 | |
33 | string beamconfigfile = "beam.cfg"; |
34 | TString genconfigfile = "gen.cfg"; |
35 | }; |
36 | |
37 | |
38 | double ampSqPT(int type, int polDir, TLorentzVector target, TLorentzVector beam, |
39 | TLorentzVector recoil,TLorentzVector q1,TLorentzVector q2); |
40 | void printUsage(genSettings_t genSettings, int goYes); |
41 | |
42 | int main(int argc, char **argv){ |
43 | |
44 | char *argptr; |
45 | |
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 | |
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 | |
70 | for (int i=1; i<argc; i++) { |
71 | argptr = argv[i]; |
72 | if (*argptr == '-') { |
73 | argptr++; |
74 | switch (*argptr) { |
75 | case 'c': |
76 | |
77 | genSettings.beamconfigfile = ++argptr; |
78 | break; |
79 | case 'd': |
80 | |
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 | |
121 | |
122 | |
123 | case 'h': |
124 | |
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++) { |
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 == "") { |
147 | cout << "No beam configuration file: run gen_ee -h for help " << endl; |
148 | exit(1); |
149 | } |
150 | |
151 | if (genSettings.genconfigfile == "") { |
152 | cout << "No generator configuration file: run gen_ee -h for help " << endl; |
153 | exit(1); |
154 | } |
155 | |
156 | |
157 | TH1D * cobrem_vs_E = 0; |
158 | if (genSettings.beamconfigfile != "") { |
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 | |
167 | if (genSettings.genconfigfile != "") { |
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]; |
| Value stored to 'A' is never read |
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){ |
192 | if (outFileSet == 1) sprintf(rootFile,genSettings.outFile); |
193 | sprintf(hddmFile,"/dev/null"); |
194 | } |
195 | if (genSettings.tOut == 2){ |
196 | if (outFileSet == 1) sprintf(hddmFile,genSettings.outFile); |
197 | sprintf(rootFile,"/dev/null"); |
198 | } |
199 | |
200 | double eGamma = genSettings.eGammaInit; |
201 | |
202 | |
203 | TH1D* hGvsE; |
204 | |
205 | |
206 | hGvsE=(TH1D*)cobrem_vs_E; |
207 | TH1D* hGvsEout = (TH1D*)hGvsE->Clone("hGvsEout"); |
208 | |
209 | |
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 | |
216 | |
217 | |
218 | |
219 | |
220 | |
221 | |
222 | HddmOut hddmGo(hddmFile); |
223 | int evtNumber = 0; |
224 | |
225 | TFile *fout = new TFile(rootFile,"RECREATE"); |
226 | |
227 | |
228 | TTree *t1 = new TTree("t1","genDevilPairs"); |
229 | devilTreePT_t devilTree; |
230 | setBranchesT1(t1, &devilTree); |
231 | |
232 | |
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 | |
243 | |
244 | |
245 | double Mcut= 5.0e-3; |
246 | if (genSettings.reaction == 3) Mcut= 20.0e-3; |
247 | if (genSettings.reaction == 2) Mcut= 1.0; |
248 | double qRcut = 1.0e-3; |
249 | if (genSettings.reaction == 2) qRcut= 2.0; |
250 | |
251 | |
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); |
255 | if (genSettings.reaction == 3) target.SetPxPyPzE(0.0,0.0,0.0,mElectron); |
256 | TLorentzVector wVec = beam + target; |
257 | TLorentzVector q1; |
258 | TLorentzVector q2; |
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); |
265 | if (genSettings.reaction == 3) recoil.SetPxPyPzE(0.0,0.0,0.0,mElectron); |
266 | |
267 | |
268 | double wMax = sqrt(target.Mag2() + 2*eGamma*target.Mag()); |
269 | if (genSettings.beamType == 2) { |
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; |
283 | PhasePT phaseGenR; |
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 | |
305 | for (int nGenTmp = 1; nGenTmp <= genSettings.nToGen; nGenTmp++) { |
306 | genVal = -1; |
307 | |
308 | if (genSettings.beamType == 2){ |
309 | yMax = gMax*1.02; |
310 | yVal = 0.0; |
311 | testValY = yMax + 10.0; |
312 | while(testValY>yVal){ |
313 | eGamma = random->Uniform(genSettings.eLower,genSettings.eUpper); |
314 | testValY = random->Uniform(0.0,yMax); |
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) { |
324 | if (nGenTmp >= nSkip) nTest++; |
325 | |
326 | genVal = phaseGenR.Gen(random); |
327 | if (genVal == -1) sumTest1 += 1.0; |
328 | if (genVal == -2) sumTest2 += 1.0; |
329 | if (genVal == -3) sumTest3 += 1.0; |
330 | if (genVal == -4) sumTest4 += 1.0; |
331 | if (genVal == -5) sumTest5 += 1.0; |
332 | |
333 | } |
334 | if (nGenTmp < nSkip) continue; |
335 | nGen++; |
336 | |
337 | |
338 | phaseSpaceWeight = phaseGenR.GetWeight(); |
339 | |
340 | |
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 | |
351 | hPhaseSpaceR->Fill(q12.Mag2(),q23.Mag2(),phaseSpaceWeight); |
352 | |
353 | |
354 | |
355 | |
356 | |
357 | |
358 | |
359 | |
360 | |
361 | |
362 | |
363 | |
364 | |
365 | |
366 | |
367 | |
368 | |
369 | |
370 | |
371 | if (genSettings.reaction == 3 && recoil.P() > q1.P()) continue; |
372 | |
373 | |
374 | |
375 | |
376 | |
377 | |
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 | |
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); |
388 | double sHfactorTrip = 1.0 - pow(fH,2); |
389 | |
390 | |
391 | double radFactor = 1.0; |
392 | double radFactorDelta = 0.0093; |
393 | double radFactorPair = 1.0+radFactorDelta; |
394 | |
395 | |
396 | |
397 | |
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) { |
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; |
405 | double radFactorTrip = 1.0 + radFactorDelta/csFrac; |
406 | |
407 | |
408 | if (genSettings.reaction == 2) { |
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) { |
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) { |
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 | |
429 | fullWeight = crossSection*phaseSpaceWeight; |
430 | |
431 | |
432 | if (fullWeight >= 0 || fullWeight <=0){ |
433 | |
434 | } else { |
435 | cout<<"!!!!Something Bad!!!!"<<endl; |
436 | } |
437 | |
438 | Double_t phaseSpaceCorrection = nGen*1.0/(1.0*nTest); |
439 | fullWeight *= phaseSpaceCorrection; |
440 | |
441 | |
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 | |
459 | sum += fullWeight; |
460 | sum2 += pow(fullWeight,2); |
461 | |
462 | if (nGen/genSettings.prescale*genSettings.prescale == nGen) { |
463 | cout <<"Integrated cross section after " << nGen << " events : " |
464 | << sum/nGen << " +/- " <<sqrt(sum2-pow(sum,2)/nGen)/nGen<< sigmaVal <<" micro barns"<< endl; |
465 | } |
466 | |
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 | |
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 | |
490 | void 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 | |
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 | |
542 | |
543 | <<" -o"<<genSettings.outFile |
544 | |
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 | |
560 | <<" -l"<<genSettings.eLower |
561 | <<" -u"<<genSettings.eUpper |
562 | <<" -o"<<genSettings.outFile |
563 | |
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 | |
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 | |
577 | |
578 | |
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 | } |