source: trunk/examples/extended/electromagnetic/TestEm10/src/Em10RunAction.cc@ 1036

Last change on this file since 1036 was 807, checked in by garnier, 17 years ago

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