source: trunk/source/geometry/magneticfield/test/field03/src/F03RunAction.cc@ 1199

Last change on this file since 1199 was 1199, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 34.1 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: F03RunAction.cc,v 1.4 2006/06/29 18:30:23 gunter Exp $
28// GEANT4 tag $Name: HEAD $
29//
30//
31
32
33#include "F03RunAction.hh"
34#include "F03RunMessenger.hh"
35
36#include "G4Run.hh"
37#include "G4UImanager.hh"
38#include "G4VVisManager.hh"
39#include "G4ios.hh"
40#include <iomanip>
41
42#include "Randomize.hh"
43#ifndef G4NOHIST
44 #include "CLHEP/Hist/HBookFile.h"
45#endif
46#include <assert.h>
47
48//////////////////////////////////////////////////////////////////////////////
49
50F03RunAction::F03RunAction()
51 :histName("histfile"),nbinStep(0),nbinEn(0),nbinTt(0),nbinTb(0),
52 nbinTsec(0),nbinTh(0),nbinThback(0),nbinR(0),nbinGamma(0),
53 nbinvertexz(0)
54{
55 runMessenger = new F03RunMessenger(this);
56 saveRndm = 1;
57#ifndef G4NOHIST
58 histo1=0; histo2=0; histo3=0; histo4=0; histo5=0;
59 histo6=0; histo7=0; histo8=0; histo9=0; histo10=0;
60#endif
61}
62
63////////////////////////////////////////////////////////////////////////////
64
65F03RunAction::~F03RunAction()
66{
67 delete runMessenger;
68#ifndef G4NOHIST
69 if(histo1) delete histo1 ;
70 if(histo2) delete histo2 ;
71 if(histo3) delete histo3 ;
72 if(histo4) delete histo4 ;
73 if(histo5) delete histo5 ;
74 if(histo6) delete histo6 ;
75 if(histo7) delete histo7 ;
76 if(histo8) delete histo8 ;
77 if(histo9) delete histo9 ;
78 if(histo10) delete histo10 ;
79 delete hbookManager;
80#endif
81}
82
83////////////////////////////////////////////////////////////////////////////////
84
85void F03RunAction::bookHisto()
86{
87#ifndef G4NOHIST
88 // init hbook
89
90 hbookManager = new HBookFile(histName, 68);
91 assert (hbookManager != 0);
92
93 // book histograms
94
95 if(nbinStep>0)
96 {
97 histo1 = hbookManager->histogram("number of steps/event"
98 ,nbinStep,Steplow,Stephigh) ;
99 assert (histo1 != 0);
100 }
101 if(nbinEn>0)
102 {
103 histo2 = hbookManager->histogram("Energy Loss (keV)"
104 ,nbinEn,Enlow/keV,Enhigh/keV) ;
105 assert (histo2 != 0);
106 }
107 if(nbinTh>0)
108 {
109 histo3 = hbookManager->histogram("angle distribution at exit(deg)"
110 ,nbinTh,Thlow/deg,Thhigh/deg) ;
111 assert (histo3 != 0);
112 }
113 if(nbinR>0)
114 {
115 histo4 = hbookManager->histogram("Lateral distribution at exit (mm)"
116 ,nbinR ,Rlow,Rhigh) ;
117 assert (histo4 != 0);
118 }
119 if(nbinTt>0)
120 {
121 histo5 = hbookManager->histogram("kinetic energy of the primary at exit(MeV)"
122 ,nbinTt,Ttlow,Tthigh) ;
123 assert (histo5 != 0);
124 }
125 if(nbinThback>0)
126 {
127 histo6 = hbookManager->histogram("angle distribution of backscattered primaries(deg)"
128 ,nbinThback,Thlowback/deg,Thhighback/deg) ;
129 assert (histo6 != 0);
130 }
131 if(nbinTb>0)
132 {
133 histo7 = hbookManager->histogram("kinetic energy of the backscattered primaries (MeV)"
134 ,nbinTb,Tblow,Tbhigh) ;
135 assert (histo7 != 0);
136 }
137 if(nbinTsec>0)
138 {
139 histo8 = hbookManager->histogram("kinetic energy of the charged secondaries (MeV)"
140 ,nbinTsec,Tseclow,Tsechigh) ;
141 assert (histo8 != 0);
142 }
143 if(nbinvertexz>0)
144 {
145 histo9 = hbookManager->histogram("z of secondary charged vertices(mm)"
146 ,nbinvertexz ,zlow,zhigh) ;
147 assert (histo9 != 0);
148 }
149 if(nbinGamma>0)
150 {
151 histo10= hbookManager->histogram("kinetic energy of gammas escaping the absorber (MeV)"
152 // ,nbinGamma,ElowGamma,EhighGamma) ;
153 ,nbinGamma,std::log10(ElowGamma),std::log10(EhighGamma)) ;
154 assert (histo10 != 0);
155 }
156#endif
157}
158
159/////////////////////////////////////////////////////////////////////////////
160
161void F03RunAction::BeginOfRunAction(const G4Run* aRun)
162{
163 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
164
165 // save Rndm status
166 if (saveRndm > 0)
167 {
168 HepRandom::showEngineStatus();
169 HepRandom::saveEngineStatus("beginOfRun.rndm");
170 }
171 G4UImanager* UI = G4UImanager::GetUIpointer();
172
173 G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance();
174
175 if(pVVisManager) UI->ApplyCommand("/vis/scene/notifyHandlers");
176
177
178 EnergySumAbs = 0. ;
179 EnergySquareSumAbs = 0.;
180 tlSumAbs = 0. ;
181 tlsquareSumAbs = 0. ;
182 nStepSumCharged = 0. ;
183 nStepSum2Charged= 0. ;
184 nStepSumNeutral = 0. ;
185 nStepSum2Neutral= 0. ;
186 TotNbofEvents = 0. ;
187 SumCharged=0.;
188 SumNeutral=0.;
189 Sum2Charged=0.;
190 Sum2Neutral=0.;
191 Selectron=0.;
192 Spositron=0.;
193
194 Transmitted=0.;
195 Reflected =0.;
196
197// plot definitions
198
199 if(nbinStep>0)
200 {
201 dStep=(Stephigh-Steplow)/nbinStep;
202 entryStep=0.;
203 underStep=0.;
204 overStep=0.;
205 for(G4int ist=0; ist<nbinStep; ist++)
206 {
207 distStep[ist]=0.;
208 }
209 }
210 if(nbinEn>0)
211 {
212 dEn = (Enhigh-Enlow)/nbinEn ;
213 entryEn=0.;
214 underEn=0.;
215 overEn=0.;
216
217 for (G4int ien=0; ien<nbinEn; ien++) distEn[ien]=0.;
218 }
219 if(nbinTt>0)
220 {
221 dTt = (Tthigh-Ttlow)/nbinTt ;
222 entryTt=0.;
223 underTt=0.;
224 overTt=0.;
225
226 for (G4int itt=0; itt<nbinTt; itt++) distTt[itt]=0.;
227
228 Ttmean=0.;
229 Tt2mean=0.;
230 }
231 if(nbinTb>0)
232 {
233 dTb = (Tbhigh-Tblow)/nbinTb ;
234 entryTb=0.;
235 underTb=0.;
236 overTb=0.;
237 for (G4int itt=0; itt<nbinTb; itt++)
238 {
239 distTb[itt]=0.;
240 }
241 Tbmean=0.;
242 Tb2mean=0.;
243 }
244 if(nbinTsec>0)
245 {
246 dTsec = (Tsechigh-Tseclow)/nbinTsec ;
247 entryTsec=0.;
248 underTsec=0.;
249 overTsec=0.;
250 for (G4int its=0; its<nbinTsec; its++)
251 {
252 distTsec[its]=0.;
253 }
254 }
255 if(nbinTh>0)
256 {
257 dTh = (Thhigh-Thlow)/nbinTh ;
258 entryTh=0.;
259 underTh=0.;
260 overTh=0.;
261 for (G4int ith=0; ith<nbinTh; ith++)
262 {
263 distTh[ith]=0.;
264 }
265 }
266
267 if(nbinThback>0)
268 {
269 dThback = (Thhighback-Thlowback)/nbinThback ;
270 entryThback=0.;
271 underThback=0.;
272 overThback=0.;
273 for (G4int ithback=0; ithback<nbinThback; ithback++)
274 {
275 distThback[ithback]=0.;
276 }
277 }
278
279
280 if(nbinR >0)
281 {
282 dR = (Rhigh-Rlow)/nbinR ;
283 entryR =0.;
284 underR =0.;
285 overR =0.;
286 for (G4int ir =0; ir <nbinR ; ir++)
287 {
288 distR[ir]=0.;
289 }
290 Rmean=0.;
291 R2mean=0.;
292 }
293
294 if(nbinGamma>0)
295 {
296 dEGamma = std::log(EhighGamma/ElowGamma)/nbinGamma ;
297 entryGamma = 0.;
298 underGamma=0.;
299 overGamma=0.;
300 for (G4int ig=0; ig<nbinGamma; ig++)
301 {
302 distGamma[ig]=0.;
303 }
304 }
305 if(nbinvertexz>0)
306 {
307 dz=(zhigh-zlow)/nbinvertexz;
308 entryvertexz=0.;
309 undervertexz=0.;
310 oververtexz=0.;
311 for(G4int iz=0; iz<nbinvertexz; iz++)
312 {
313 distvertexz[iz]=0.;
314 }
315 }
316
317 bookHisto();
318}
319
320/////////////////////////////////////////////////////////////////////////////
321
322void F03RunAction::EndOfRunAction(const G4Run* aRun)
323{
324 G4double sAbs,sigAbs,sigstep,sigcharged,signeutral;
325
326 tlSumAbs /= TotNbofEvents ;
327 sAbs = tlsquareSumAbs/TotNbofEvents-tlSumAbs*tlSumAbs ;
328 if(sAbs>0.)
329 sAbs = std::sqrt(sAbs/TotNbofEvents) ;
330 else
331 sAbs = 0. ;
332
333 EnergySumAbs /= TotNbofEvents ;
334 sigAbs = EnergySquareSumAbs/TotNbofEvents-EnergySumAbs*EnergySumAbs;
335 if(sigAbs>0.)
336 sigAbs = std::sqrt(sigAbs/TotNbofEvents);
337 else
338 sigAbs = 0.;
339
340 nStepSumCharged /= TotNbofEvents ;
341 sigstep = nStepSum2Charged/TotNbofEvents-nStepSumCharged*nStepSumCharged;
342 if(sigstep>0.)
343 sigstep = std::sqrt(sigstep/TotNbofEvents);
344 else
345 sigstep = 0.;
346 G4double sigch=sigstep ;
347
348 nStepSumNeutral /= TotNbofEvents ;
349 sigstep = nStepSum2Neutral/TotNbofEvents-nStepSumNeutral*nStepSumNeutral;
350 if(sigstep>0.)
351 sigstep = std::sqrt(sigstep/TotNbofEvents);
352 else
353 sigstep = 0.;
354 G4double signe=sigstep ;
355
356 SumCharged /= TotNbofEvents;
357 sigcharged = Sum2Charged/TotNbofEvents-SumCharged*SumCharged;
358 if(sigcharged>0.)
359 sigcharged = std::sqrt(sigcharged/TotNbofEvents);
360 else
361 sigcharged = 0. ;
362
363 SumNeutral /= TotNbofEvents;
364 signeutral = Sum2Neutral/TotNbofEvents-SumNeutral*SumNeutral;
365 if(signeutral>0.)
366 signeutral = std::sqrt(signeutral/TotNbofEvents);
367 else
368 signeutral = 0. ;
369
370 Selectron /= TotNbofEvents ;
371 Spositron /= TotNbofEvents ;
372
373 Transmitted /=TotNbofEvents ;
374 Reflected /=TotNbofEvents ;
375 G4cout << " ================== run summary =====================" << G4endl;
376 G4int prec = G4cout.precision(6);
377 G4cout << " end of Run TotNbofEvents = " <<
378 TotNbofEvents << G4endl ;
379 G4cout << " mean charged track length in absorber=" <<
380 tlSumAbs/mm << " +- " << sAbs/mm <<
381 " mm " << G4endl;
382 G4cout << G4endl;
383 G4cout << " mean energy deposit in absorber=" <<
384 EnergySumAbs/MeV << " +- " << sigAbs/MeV <<
385 " MeV " << G4endl ;
386 G4cout << G4endl ;
387 G4cout << " mean number of steps in absorber (charged) =" <<
388 nStepSumCharged << " +- " << sigch <<
389 " " << G4endl ;
390 G4cout << " mean number of steps in absorber (neutral) =" <<
391 nStepSumNeutral << " +- " << signe <<
392 " " << G4endl ;
393 G4cout << G4endl ;
394 G4cout << " mean number of charged secondaries = " <<
395 SumCharged << " +- " << sigcharged << G4endl;
396 G4cout << G4endl ;
397 G4cout << " mean number of neutral secondaries = " <<
398 SumNeutral << " +- " << signeutral << G4endl;
399 G4cout << G4endl ;
400
401 G4cout << " mean number of e-s =" << Selectron <<
402 " and e+s =" << Spositron << G4endl;
403 G4cout << G4endl;
404
405 G4cout << "(number) transmission coeff=" << Transmitted <<
406 " reflection coeff=" << Reflected << G4endl;
407 G4cout << G4endl;
408
409 if(nbinStep>0)
410 {G4double E , dnorm, norm ;
411 G4cout << " step number/event distribution " << G4endl ;
412 G4cout << "#entries=" << entryStep << " #underflows=" << underStep <<
413 " #overflows=" << overStep << G4endl ;
414 if( entryStep>0.)
415 {
416 E = Steplow - dStep ;
417 norm = TotNbofEvents ;
418 G4cout << " bin nb nsteplow entries normalized " << G4endl ;
419 for(G4int iss=0; iss<nbinStep; iss++)
420 {
421 E += dStep ;
422 dnorm = distStep[iss]/norm;
423 G4cout << std::setw(5) << iss << std::setw(10) << E <<
424 std::setw(12) << distStep[iss] <<
425 std::setw(12) << dnorm << G4endl ;
426 }
427 G4cout << G4endl;
428 }
429 }
430 if(nbinEn > 0)
431 {
432 std::ofstream fileOut("distribution.out", std::ios::out ) ;
433 fileOut.setf( std::ios::scientific, std::ios::floatfield );
434
435 std::ofstream normOut("normDist.out", std::ios::out ) ;
436 normOut.setf( std::ios::scientific, std::ios::floatfield );
437
438 G4double E , dnorm, norm,fmax,Emp,width ;
439 Emp=-999.999 ;
440 G4cout << " energy deposit distribution " << G4endl ;
441 G4cout << "#entries=" << entryEn << " #underflows=" << underEn <<
442 " #overflows=" << overEn << G4endl ;
443 if( entryEn>0.)
444 {
445 E = Enlow - dEn ;
446 norm = TotNbofEvents*1.0 ; // *dEn ;
447 G4cout << " bin nb Elow entries normalized " << G4endl ;
448 fmax = 0. ;
449
450 for(G4int ien=0; ien<nbinEn; ien++)
451 {
452 E += dEn ;
453
454 if(distEn[ien]>fmax)
455 {
456 fmax = distEn[ien] ;
457 Emp = E ; // most probable roughly
458 }
459 dnorm = distEn[ien]/norm;
460
461 G4cout << std::setw(5) << ien << std::setw(10) << E/keV <<
462 std::setw(12) << distEn[ien] <<
463 std::setw(12) << dnorm << G4endl ;
464
465 fileOut << E/keV << "\t"<< distEn[ien] << G4endl ;
466 normOut << E/keV << "\t"<< dnorm << G4endl ;
467 }
468 G4cout << G4endl;
469 G4int ii ;
470 G4double E1,E2 ;
471 E1=-1.e6 ;
472 E2=+1.e6 ;
473 E = Enlow -dEn ;
474 ii = -1;
475
476 for(G4int i1=0; i1<nbinEn; i1++)
477 {
478 E += dEn ;
479 if(ii<0)
480 {
481 if(distEn[i1] >= 0.5*fmax)
482 {
483 E1=E ;
484 ii=i1 ;
485 }
486 }
487 }
488 E = Enlow -dEn ;
489
490 for(G4int i2=0; i2<nbinEn; i2++)
491 {
492 E += dEn ;
493
494 if(distEn[i2] >= 0.5*fmax) E2=E ;
495 }
496 G4cout << " Emp = " << std::setw(15) << Emp/MeV << " width="
497 << std::setw(15) << (E2-E1)/MeV << " MeV " << G4endl;
498 G4cout << G4endl ;
499 }
500 }
501 if(nbinTt>0)
502 {
503 G4double E , dnorm, norm ,sig;
504 G4cout << " transmitted energy distribution " << G4endl ;
505 G4cout << "#entries=" << entryTt << " #underflows=" << underTt <<
506 " #overflows=" << overTt << G4endl ;
507 if( entryTt>0.)
508 {
509 Ttmean /= entryTt;
510 sig=Tt2mean/entryTt-Ttmean*Ttmean ;
511 if(sig<=0.)
512 sig=0.;
513 else
514 sig=std::sqrt(sig/entryTt) ;
515 G4cout << " mean energy of transmitted particles=" << Ttmean/keV <<
516 " +- " << sig/keV << " keV." << G4endl;
517 E = Ttlow - dTt ;
518 norm = TotNbofEvents*dTt ;
519 G4cout << " bin nb Elow entries normalized " << G4endl ;
520 for(G4int itt=0; itt<nbinTt; itt++)
521 {
522 E += dTt ;
523 dnorm = distTt[itt]/norm;
524 G4cout << std::setw(5) << itt << std::setw(10) << E <<
525 std::setw(12) << distTt[itt] <<
526 std::setw(12) << dnorm << G4endl ;
527 }
528 G4cout << G4endl;
529 }
530 }
531 if(nbinTb>0)
532 {
533 G4double E , dnorm, norm ,sig;
534 G4cout << " backscattered energy distribution " << G4endl ;
535 G4cout << "#entries=" << entryTb << " #underflows=" << underTb <<
536 " #overflows=" << overTb << G4endl ;
537 if( entryTb>0.)
538 {
539 Tbmean /= entryTb;
540 sig=Tb2mean/entryTb-Tbmean*Tbmean ;
541 if(sig<=0.)
542 sig=0.;
543 else
544 sig=std::sqrt(sig/entryTb) ;
545 G4cout << " mean energy of backscattered particles=" << Tbmean/keV <<
546 " +- " << sig/keV << " keV." << G4endl;
547 E = Tblow - dTb ;
548 norm = TotNbofEvents*dTb ;
549 G4cout << " bin nb Elow entries normalized " << G4endl ;
550 for(G4int itt=0; itt<nbinTb; itt++)
551 {
552 E += dTb ;
553 dnorm = distTb[itt]/norm;
554 G4cout << std::setw(5) << itt << std::setw(10) << E <<
555 std::setw(12) << distTb[itt] <<
556 std::setw(12) << dnorm << G4endl ;
557 }
558 G4cout << G4endl;
559 }
560 }
561 if(nbinTsec>0)
562 {G4double E , dnorm, norm ;
563 G4cout << " energy distribution of charged secondaries " << G4endl ;
564 G4cout << "#entries=" << entryTsec << " #underflows=" << underTsec <<
565 " #overflows=" << overTsec << G4endl ;
566 if( entryTsec>0.)
567 {
568 E = Tseclow - dTsec ;
569 norm = TotNbofEvents*dTsec ;
570 G4cout << " bin nb Elow entries normalized " << G4endl ;
571 for(G4int itt=0; itt<nbinTsec; itt++)
572 {
573 E += dTsec ;
574 dnorm = distTsec[itt]/norm;
575 G4cout << std::setw(5) << itt << std::setw(10) << E <<
576 std::setw(12) << distTsec[itt] <<
577 std::setw(12) << dnorm << G4endl ;
578 }
579 G4cout << G4endl;
580 }
581 }
582
583 if(nbinR >0)
584 {
585 G4double R , dnorm, norm,sig ;
586
587 std::ofstream fileOut("diameter.out", std::ios::out ) ;
588 fileOut.setf( std::ios::scientific, std::ios::floatfield );
589
590 std::ofstream normOut("normDiameter.out", std::ios::out ) ;
591 normOut.setf( std::ios::scientific, std::ios::floatfield );
592
593 G4cout << " R distribution " << G4endl ;
594 G4cout << "#entries=" << entryR << " #underflows=" << underR <<
595 " #overflows=" << overR << G4endl ;
596 if( entryR >0.)
597 {
598 Rmean /= entryR;
599 sig = R2mean/entryR - Rmean*Rmean;
600 if(sig<=0.) sig=0. ;
601 else sig = std::sqrt(sig/entryR) ;
602 G4cout << " mean lateral displacement at exit=" << Rmean/mm << " +- "
603 << sig/mm << " mm." << G4endl ;
604 R = Rlow - dR ;
605 norm = TotNbofEvents*dR ;
606 G4cout << " bin nb Rlow entries normalized " << G4endl ;
607 for(G4int ier=0; ier<nbinR ; ier++)
608 {
609 R+= dR ;
610 dnorm = distR[ier]/norm;
611 G4cout << std::setw(5) << ier << std::setw(10) << R <<
612 std::setw(12) << distR[ier] <<
613 std::setw(12) << dnorm << G4endl ;
614 fileOut << R/mm << "\t"<< distR[ier] << G4endl ;
615 normOut << R/mm << "\t"<< dnorm << G4endl ;
616 }
617 G4cout << G4endl;
618 }
619 }
620
621 if(nbinTh>0)
622 {G4double Th,Thdeg, dnorm, norm,fac0,fnorm,pere,Thpere,Thmean,sum;
623 G4cout << " angle distribution " << G4endl ;
624 G4cout << "#entries=" << entryTh << " #underflows=" << underTh <<
625 " #overflows=" << overTh << G4endl ;
626 if( entryTh>0.)
627 {
628 Th= Thlow - dTh ;
629 norm = TotNbofEvents ;
630 if(distTh[0] == 0.)
631 fac0 = 1. ;
632 else
633 fac0 = 1./distTh[0] ;
634 pere = 1./std::exp(1.) ;
635
636 G4cout << " bin nb Thlowdeg Thlowrad " <<
637 " entries normalized " << G4endl ;
638 Thpere = 0. ;
639 sum = 0. ;
640 Thmean = 0. ;
641 for(G4int ien=0; ien<nbinTh; ien++)
642 {
643 Th+= dTh ;
644 Thdeg = Th*180./pi ;
645 sum += distTh[ien] ;
646 Thmean += distTh[ien]*(Th+0.5*dTh) ;
647 dnorm = distTh[ien]/norm;
648 fnorm = fac0*distTh[ien] ;
649 if( fnorm > pere)
650 Thpere = Th ;
651 G4cout << std::setw(5) << ien << std::setw(10) << Thdeg << " " <<
652 std::setw(10) << Th << " " <<
653 std::setw(12) << distTh[ien] << " " <<
654 std::setw(12) << dnorm << " " << std::setw(12) << fnorm <<G4endl ;
655 }
656 Thmean /= sum ;
657 G4cout << G4endl;
658 G4cout << " mean = " << Thmean << " rad or " << 180.*Thmean/pi <<
659 " deg." << G4endl;
660 G4cout << " theta(1/e)=" << Thpere << " - " << Thpere+dTh << " rad "
661 << " or " << 180.*Thpere/pi << " - " << 180.*(Thpere+dTh)/pi
662 << " deg." << G4endl;
663 G4cout << G4endl;
664 }
665 }
666
667 if(nbinThback>0)
668 {G4double Thb,Thdegb, dnormb, normb,fac0b,fnormb,pereb,Thpereb,Thmeanb,sumb;
669 G4cout << " backscattering angle distribution " << G4endl ;
670 G4cout << "#entries=" << entryThback << " #underflows=" << underThback <<
671 " #overflows=" << overThback << G4endl ;
672 if( entryThback>0.)
673 {
674 Thb= Thlowback - dThback ;
675 normb = TotNbofEvents ;
676 if(distThback[0] == 0.)
677 fac0b = 1. ;
678 else
679 fac0b = 1./distThback[0] ;
680 pereb = 1./std::exp(1.) ;
681
682 G4cout << " bin nb Thlowdeg Thlowrad " <<
683 " entries normalized " << G4endl ;
684 Thpereb = 0. ;
685 sumb = 0. ;
686 Thmeanb = 0. ;
687 for(G4int ien=0; ien<nbinThback; ien++)
688 {
689 Thb+= dThback ;
690 Thdegb = Thb*180./pi ;
691 sumb += distThback[ien] ;
692 Thmeanb += distThback[ien]*(Thb+0.5*dThback) ;
693 dnormb = distThback[ien]/normb;
694 fnormb = fac0b*distThback[ien] ;
695 if( fnormb > pereb)
696 Thpereb = Thb ;
697 G4cout << std::setw(5) << ien << std::setw(10) << Thdegb << " " <<
698 std::setw(10) << Thb << " " <<
699 std::setw(12) << distThback[ien] << " " <<
700 std::setw(12) << dnormb << " " << std::setw(12) << fnormb <<G4endl ;
701 }
702 Thmeanb /= sumb ;
703 G4cout << G4endl;
704 G4cout << " mean = " << Thmeanb << " rad or " << 180.*Thmeanb/pi <<
705 " deg." << G4endl;
706 G4cout << " theta(1/e)=" << Thpereb << " - " << Thpereb+dThback << " rad "
707 << " or " << 180.*Thpereb/pi << " - " << 180.*(Thpereb+dThback)/pi
708 << " deg." << G4endl;
709 G4cout << G4endl;
710 }
711 }
712
713 if(nbinGamma>0)
714 {G4double E , fact,dnorm, norm ;
715 G4cout << " gamma energy distribution " << G4endl ;
716 G4cout << "#entries=" << entryGamma << " #underflows=" << underGamma <<
717 " #overflows=" << overGamma << G4endl ;
718 if( entryGamma>0.)
719 {
720 fact=std::exp(dEGamma) ;
721 E = ElowGamma/fact ;
722 norm = TotNbofEvents*dEGamma;
723 G4cout << " bin nb Elow entries normalized " << G4endl ;
724 for(G4int itt=0; itt<nbinGamma; itt++)
725 {
726 E *= fact ;
727 dnorm = distGamma[itt]/norm;
728 G4cout << std::setw(5) << itt << std::setw(13) << E <<
729 std::setw(12) << distGamma[itt] <<
730 std::setw(15) << dnorm << G4endl ;
731 }
732 G4cout << G4endl;
733 }
734 }
735
736 if(nbinvertexz >0)
737 {G4double z , dnorm, norm ;
738 G4cout << " vertex Z distribution " << G4endl ;
739 G4cout << "#entries=" << entryvertexz << " #underflows=" << undervertexz <<
740 " #overflows=" << oververtexz << G4endl ;
741 if( entryvertexz >0.)
742 {
743 z =zlow - dz ;
744 norm = TotNbofEvents*dz ;
745 G4cout << " bin nb zlow entries normalized " << G4endl ;
746 for(G4int iez=0; iez<nbinvertexz ; iez++)
747 {
748 z+= dz ;
749 if(std::abs(z)<1.e-12) z=0.;
750 dnorm = distvertexz[iez]/norm;
751 G4cout << std::setw(5) << iez << std::setw(10) << z <<
752 std::setw(12) << distvertexz[iez] <<
753 std::setw(12) << dnorm << G4endl ;
754 }
755 G4cout << G4endl;
756 }
757 }
758
759 G4cout.precision(prec);
760
761 if (G4VVisManager::GetConcreteInstance())
762 {
763 G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/update");
764 }
765 // Write histogram file
766
767#ifndef G4NOHIST
768 hbookManager->write();
769#endif
770
771 // save Rndm status
772
773 if (saveRndm == 1)
774 {
775 HepRandom::showEngineStatus();
776 HepRandom::saveEngineStatus("endOfRun.rndm");
777 }
778}
779
780///////////////////////////////////////////////////////////////////////////
781
782void F03RunAction::CountEvent()
783{
784 TotNbofEvents += 1. ;
785}
786
787/////////////////////////////////////////////////////////////////////////
788
789void F03RunAction::AddnStepsCharged(G4double ns)
790{
791 nStepSumCharged += ns;
792 nStepSum2Charged += ns*ns;
793}
794
795////////////////////////////////////////////////////////////////////////
796
797void F03RunAction::AddnStepsNeutral(G4double ns)
798{
799 nStepSumNeutral += ns;
800 nStepSum2Neutral += ns*ns;
801}
802
803////////////////////////////////////////////////////////////////////////////
804
805void F03RunAction::AddEdeps(G4double Eabs)
806{
807 EnergySumAbs += Eabs;
808 EnergySquareSumAbs += Eabs*Eabs;
809}
810
811/////////////////////////////////////////////////////////////////////////////
812
813void F03RunAction::AddTrackLength(G4double tlabs)
814{
815 tlSumAbs += tlabs;
816 tlsquareSumAbs += tlabs*tlabs ;
817}
818
819/////////////////////////////////////////////////////////////////////////////
820
821void F03RunAction::AddTrRef(G4double tr,G4double ref)
822{
823 Transmitted += tr ;
824 Reflected += ref;
825}
826
827/////////////////////////////////////////////////////////////////////////////
828
829void F03RunAction::FillNbOfSteps(G4double ns)
830{
831#ifndef G4NOHIST
832 const G4double eps = 1.e-10 ;
833 G4double n,bin ;
834 G4int ibin;
835
836 if(histo1)
837 {
838 entryStep += 1. ;
839
840 if(ns<Steplow)
841 underStep += 1. ;
842 else if(ns>=Stephigh)
843 overStep += 1. ;
844 else
845 {
846 n = ns+eps ;
847 bin = (n-Steplow)/dStep ;
848 ibin= (G4int)bin ;
849 distStep[ibin] += 1. ;
850 }
851 histo1->accumulate(ns) ;
852 }
853#endif
854}
855
856//////////////////////////////////////////////////////////////////////////////
857
858void F03RunAction::FillEn(G4double En)
859{
860#ifndef G4NOHIST
861 G4double bin ;
862 G4int ibin;
863
864 if(histo2)
865 {
866 entryEn += 1. ;
867
868 if(En < Enlow) underEn += 1. ;
869 else if( En >= Enhigh) overEn += 1. ;
870 else
871 {
872 bin = (En-Enlow)/dEn ;
873 ibin= (G4int)bin ;
874 distEn[ibin] += 1. ;
875 }
876 histo2->accumulate(En/keV) ; // was /MeV
877 }
878#endif
879}
880
881////////////////////////////////////////////////////////////////////////////
882
883void F03RunAction::FillTt(G4double En)
884{
885#ifndef G4NOHIST
886 G4double bin ;
887 G4int ibin;
888
889 if(histo5)
890 {
891 entryTt += 1. ;
892 Ttmean += En ;
893 Tt2mean += En*En ;
894
895 if(En<Ttlow)
896 underTt += 1. ;
897 else if(En>=Tthigh)
898 overTt += 1. ;
899 else
900 {
901 bin = (En-Ttlow)/dTt ;
902 ibin= (G4int)bin ;
903 distTt[ibin] += 1. ;
904 }
905 histo5->accumulate(En/MeV) ;
906 }
907#endif
908}
909
910//////////////////////////////////////////////////////////////////////////////
911
912void F03RunAction::FillTb(G4double En)
913{
914#ifndef G4NOHIST
915 G4double bin ;
916 G4int ibin;
917
918 if(histo7)
919 {
920 entryTb += 1. ;
921 Tbmean += En ;
922 Tb2mean += En*En ;
923
924 if(En<Tblow)
925 underTb += 1. ;
926 else if(En>=Tbhigh)
927 overTb += 1. ;
928 else
929 {
930 bin = (En-Tblow)/dTb ;
931 ibin= (G4int)bin ;
932 distTb[ibin] += 1. ;
933 }
934 histo7->accumulate(En/MeV) ;
935 }
936#endif
937}
938
939///////////////////////////////////////////////////////////////////////////////
940
941void F03RunAction::FillTsec(G4double En)
942{
943#ifndef G4NOHIST
944 G4double bin ;
945 G4int ibin;
946
947 if(histo8)
948 {
949 entryTsec += 1. ;
950
951 if(En<Tseclow)
952 underTsec += 1. ;
953 else if(En>=Tsechigh)
954 overTsec += 1. ;
955 else
956 {
957 bin = (En-Tseclow)/dTsec ;
958 ibin= (G4int)bin ;
959 distTsec[ibin] += 1. ;
960 }
961 histo8->accumulate(En/MeV) ;
962 }
963#endif
964}
965
966/////////////////////////////////////////////////////////////////////////////
967
968void F03RunAction::FillGammaSpectrum(G4double En)
969{
970#ifndef G4NOHIST
971 G4double bin ;
972 G4int ibin;
973
974 if(histo10)
975 {
976 entryGamma += 1. ;
977
978 if(En<ElowGamma)
979 underGamma += 1. ;
980 else if(En>=EhighGamma)
981 overGamma += 1. ;
982 else
983 {
984 bin = std::log(En/ElowGamma)/dEGamma;
985 ibin= (G4int)bin ;
986 distGamma[ibin] += 1. ;
987 }
988 histo10->accumulate(std::log10(En/MeV)) ;
989 }
990#endif
991}
992
993////////////////////////////////////////////////////////////////////////////////
994
995void F03RunAction::FillTh(G4double Th)
996{
997#ifndef G4NOHIST
998 static const G4double cn=pi/(64800.*dTh) ;
999 static const G4double cs=pi/
1000 (64800.*(std::cos(Thlow)-std::cos(Thlow+dTh)));
1001 G4double bin,Thbin ,wg;
1002 G4int ibin;
1003
1004 if(histo3)
1005 {
1006 entryTh += 1. ;
1007
1008 wg = 0.;
1009
1010 if(Th<Thlow)
1011 underTh += 1. ;
1012 else if(Th>=Thhigh)
1013 overTh += 1. ;
1014 else
1015 {
1016 bin = (Th-Thlow)/dTh ;
1017 ibin= (G4int)bin ;
1018 Thbin = Thlow+ibin*dTh ;
1019 if(Th > 0.001*dTh)
1020 wg=cn/std::sin(Th) ;
1021 else
1022 {
1023 G4double thdeg=Th*180./pi;
1024 G4cout << "theta < 0.001*dth (from plot excluded) theta="
1025 << std::setw(12) << std::setprecision(4) << thdeg << G4endl;
1026 wg=0. ;
1027 }
1028 distTh[ibin] += wg ;
1029 }
1030
1031 histo3->accumulate(Th/deg, wg) ;
1032 }
1033#endif
1034}
1035
1036//////////////////////////////////////////////////////////////////////////
1037
1038void F03RunAction::FillThBack(G4double Th)
1039{
1040#ifndef G4NOHIST
1041 static const G4double cn=pi/(64800.*dThback) ;
1042 static const G4double cs=pi/
1043 (64800.*(std::cos(Thlowback)-std::cos(Thlowback+dThback)));
1044 G4double bin,Thbin,wg ;
1045 G4int ibin;
1046
1047 if(histo6)
1048 {
1049 entryThback += 1. ;
1050
1051 if(Th<Thlowback)
1052 underThback += 1. ;
1053 else if(Th>=Thhighback)
1054 overThback += 1. ;
1055 else
1056 {
1057 bin = (Th-Thlowback)/dThback ;
1058 ibin= (G4int)bin ;
1059 Thbin = Thlowback+ibin*dThback ;
1060 if(Th > 0.001*dThback)
1061 wg=cn/std::sin(Th) ;
1062 else
1063 {
1064 G4double thdeg=Th*180./pi;
1065 G4cout << "theta < 0.001*dth (from plot excluded) theta="
1066 << std::setw(12) << std::setprecision(4) << thdeg << G4endl;
1067 wg=0. ;
1068 }
1069 distThback[ibin] += wg ;
1070 }
1071 histo6->accumulate(Th/deg, wg) ;
1072 }
1073#endif
1074}
1075
1076//////////////////////////////////////////////////////////////////////////
1077
1078void F03RunAction::FillR(G4double R )
1079{
1080#ifndef G4NOHIST
1081 G4double bin ;
1082 G4int ibin;
1083
1084 if(histo4)
1085 {
1086 entryR += 1. ;
1087 Rmean += R ;
1088 R2mean += R*R ;
1089
1090 if(R <Rlow)
1091 underR += 1. ;
1092 else if(R >=Rhigh)
1093 overR += 1. ;
1094 else
1095 {
1096 bin = (R -Rlow)/dR ;
1097 ibin= (G4int)bin ;
1098 distR[ibin] += 1. ;
1099 }
1100 histo4->accumulate(R/mm) ;
1101 }
1102#endif
1103}
1104
1105/////////////////////////////////////////////////////////////////////////////
1106
1107void F03RunAction::Fillvertexz(G4double z )
1108{
1109#ifndef G4NOHIST
1110 G4double bin ;
1111 G4int ibin;
1112
1113 if(histo9)
1114 {
1115 entryvertexz += 1. ;
1116
1117 if(z <zlow)
1118 undervertexz += 1. ;
1119 else if(z >=zhigh)
1120 oververtexz += 1. ;
1121 else
1122 {
1123 bin = (z -zlow)/dz ;
1124 ibin = (G4int)bin ;
1125 distvertexz[ibin] += 1. ;
1126 }
1127 histo9->accumulate(z/mm) ;
1128 }
1129#endif
1130}
1131
1132//////////////////////////////////////////////////////////////////////////////
1133
1134void F03RunAction::SethistName(G4String name)
1135{
1136 histName = name ;
1137 G4cout << " hist file = " << histName << G4endl;
1138}
1139
1140void F03RunAction::SetnbinStep(G4int nbin)
1141{
1142 nbinStep = nbin ;
1143 if(nbinStep>0)
1144 G4cout << " Nb of bins in #step plot = " << nbinStep << G4endl ;
1145}
1146
1147void F03RunAction::SetSteplow(G4double low)
1148{
1149 Steplow = low ;
1150 if(nbinStep>0)
1151 G4cout << " low in the #step plot = " << Steplow << G4endl ;
1152}
1153void F03RunAction::SetStephigh(G4double high)
1154{
1155 Stephigh = high ;
1156 if(nbinStep>0)
1157 G4cout << " high in the #step plot = " << Stephigh << G4endl ;
1158}
1159
1160////////////////////////////////////////////////////////////////////////
1161
1162void F03RunAction::SetnbinEn(G4int nbin)
1163{
1164 nbinEn = nbin ;
1165
1166 if(nbinEn > 0) G4cout << " Nb of bins in Edep plot = " << nbinEn << G4endl ;
1167}
1168
1169void F03RunAction::SetEnlow(G4double Elow)
1170{
1171 Enlow = Elow ;
1172 if(nbinEn>0) G4cout << " Elow in the Edep plot = " << Enlow << G4endl ;
1173}
1174
1175void F03RunAction::SetEnhigh(G4double Ehigh)
1176{
1177 Enhigh = Ehigh ;
1178
1179 if(nbinEn>0) G4cout << " Ehigh in the Edep plot = " << Enhigh << G4endl ;
1180}
1181
1182/////////////////////////////////////////////////////////////////////////
1183
1184void F03RunAction::SetnbinGamma(G4int nbin)
1185{
1186 nbinGamma = nbin ;
1187 if(nbinGamma>0)
1188 G4cout << " Nb of bins in gamma spectrum plot = " << nbinGamma << G4endl ;
1189}
1190
1191void F03RunAction::SetElowGamma(G4double Elow)
1192{
1193 ElowGamma = Elow ;
1194 if(nbinGamma>0)
1195 G4cout << " Elow in the gamma spectrum plot = " << ElowGamma << G4endl ;
1196}
1197
1198void F03RunAction::SetEhighGamma(G4double Ehigh)
1199{
1200 EhighGamma = Ehigh ;
1201 if(nbinGamma>0)
1202 G4cout << " Ehigh in the gamma spectrum plot = " << EhighGamma << G4endl ;
1203}
1204
1205void F03RunAction::SetnbinTt(G4int nbin)
1206{
1207 nbinTt = nbin ;
1208 if(nbinTt>0)
1209 G4cout << " Nb of bins in Etransmisssion plot = " << nbinTt << G4endl ;
1210}
1211
1212void F03RunAction::SetTtlow(G4double Elow)
1213{
1214 Ttlow = Elow ;
1215 if(nbinTt>0)
1216 G4cout << " Elow in the Etransmission plot = " << Ttlow << G4endl ;
1217}
1218
1219void F03RunAction::SetTthigh(G4double Ehigh)
1220{
1221 Tthigh = Ehigh ;
1222 if(nbinTt>0)
1223 G4cout << " Ehigh in the Etransmission plot = " << Tthigh << G4endl ;
1224}
1225
1226void F03RunAction::SetnbinTb(G4int nbin)
1227{
1228 nbinTb = nbin ;
1229 if(nbinTb>0)
1230 G4cout << " Nb of bins in Ebackscattered plot = " << nbinTb << G4endl ;
1231}
1232
1233void F03RunAction::SetTblow(G4double Elow)
1234{
1235 Tblow = Elow ;
1236 if(nbinTb>0)
1237 G4cout << " Elow in the Ebackscattered plot = " << Tblow << G4endl ;
1238}
1239
1240void F03RunAction::SetTbhigh(G4double Ehigh)
1241{
1242 Tbhigh = Ehigh ;
1243 if(nbinTb>0)
1244 G4cout << " Ehigh in the Ebackscattered plot = " << Tbhigh << G4endl ;
1245}
1246
1247void F03RunAction::SetnbinTsec(G4int nbin)
1248{
1249 nbinTsec = nbin ;
1250 if(nbinTsec>0)
1251 G4cout << " Nb of bins in Tsecondary plot = " << nbinTsec << G4endl ;
1252}
1253
1254void F03RunAction::SetTseclow(G4double Elow)
1255{
1256 Tseclow = Elow ;
1257 if(nbinTsec>0)
1258 G4cout << " Elow in the Tsecondary plot = " << Tseclow << G4endl ;
1259}
1260
1261void F03RunAction::SetTsechigh(G4double Ehigh)
1262{
1263 Tsechigh = Ehigh ;
1264 if(nbinTsec>0)
1265 G4cout << " Ehigh in the Tsecondary plot = " << Tsechigh << G4endl ;
1266}
1267
1268void F03RunAction::SetnbinR(G4int nbin)
1269{
1270 nbinR = nbin ;
1271 if(nbinR>0)
1272 G4cout << " Nb of bins in R plot = " << nbinR << G4endl ;
1273}
1274
1275void F03RunAction::SetRlow(G4double rlow)
1276{
1277 Rlow = rlow ;
1278 if(nbinR>0)
1279 G4cout << " Rlow in the R plot = " << Rlow << G4endl ;
1280}
1281
1282void F03RunAction::SetRhigh(G4double rhigh)
1283{
1284 Rhigh = rhigh ;
1285 if(nbinR>0)
1286 G4cout << " Rhigh in the R plot = " << Rhigh << G4endl ;
1287}
1288
1289void F03RunAction::Setnbinzvertex(G4int nbin)
1290{
1291 nbinvertexz = nbin ;
1292 if(nbinvertexz>0)
1293 G4cout << " Nb of bins in Z plot = " << nbinvertexz << G4endl ;
1294}
1295
1296void F03RunAction::Setzlow(G4double z)
1297{
1298 zlow = z ;
1299 if(nbinvertexz>0)
1300 G4cout << " zlow in the Z plot = " << zlow << G4endl ;
1301}
1302
1303void F03RunAction::Setzhigh(G4double z)
1304{
1305 zhigh = z ;
1306 if(nbinvertexz>0)
1307 G4cout << " zhigh in the Z plot = " << zhigh << G4endl ;
1308}
1309
1310void F03RunAction::SetnbinTh(G4int nbin)
1311{
1312 nbinTh = nbin ;
1313 if(nbinTh>0)
1314 G4cout << " Nb of bins in Theta plot = " << nbinTh << G4endl ;
1315}
1316
1317void F03RunAction::SetThlow(G4double Tlow)
1318{
1319 Thlow = Tlow ;
1320 if(nbinTh>0)
1321 G4cout << " Tlow in the Theta plot = " << Thlow << G4endl ;
1322}
1323
1324void F03RunAction::SetThhigh(G4double Thigh)
1325{
1326 Thhigh = Thigh ;
1327 if(nbinTh>0)
1328 G4cout << " Thigh in the Theta plot = " << Thhigh << G4endl ;
1329}
1330
1331void F03RunAction::SetnbinThBack(G4int nbin)
1332{
1333 nbinThback = nbin ;
1334 if(nbinThback>0)
1335 G4cout << " Nb of bins in Theta plot = " << nbinThback << G4endl ;
1336}
1337
1338void F03RunAction::SetThlowBack(G4double Tlow)
1339{
1340 Thlowback = Tlow ;
1341 if(nbinThback>0)
1342 G4cout << " Tlow in the Theta plot = " << Thlowback << G4endl ;
1343}
1344
1345void F03RunAction::SetThhighBack(G4double Thigh)
1346{
1347 Thhighback = Thigh ;
1348 if(nbinThback>0)
1349 G4cout << " Thigh in the Theta plot = " << Thhighback << G4endl ;
1350}
1351
1352void F03RunAction::CountParticles(G4double nch,G4double nne)
1353{
1354 SumCharged += nch ;
1355 SumNeutral += nne ;
1356 Sum2Charged += nch*nch ;
1357 Sum2Neutral += nne*nne ;
1358}
1359
1360void F03RunAction::AddEP(G4double nele,G4double npos)
1361{
1362 Selectron += nele;
1363 Spositron += npos;
1364}
1365
1366//
1367//
1368////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.