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

Last change on this file since 1281 was 1230, checked in by garnier, 14 years ago

update to geant4.9.3
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: geant4-09-03-cand-01 $
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.