source: trunk/examples/extended/electromagnetic/TestEm8/src/Em8RunAction.cc @ 1279

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

update to geant4.9.3
File size: 26.8 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: Em8RunAction.cc,v 1.15 2007/11/12 10:54:49 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30//
31
32
33#include "Em8RunAction.hh"
34#include "Em8RunMessenger.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//////////////////////////////////////////////////////////////////////////////
46
47Em8RunAction::Em8RunAction()
48  :histName("histfile"),nbinStep(0),nbinEn(0),nbinTt(0),nbinTb(0),
49   nbinTsec(0),nbinTh(0),nbinThback(0),nbinR(0),nbinGamma(0),
50   nbinvertexz(0)
51{
52  runMessenger = new Em8RunMessenger(this);
53  saveRndm = 1; 
54}
55
56////////////////////////////////////////////////////////////////////////////
57
58Em8RunAction::~Em8RunAction()
59{
60  delete runMessenger;
61}
62
63
64/////////////////////////////////////////////////////////////////////////////
65
66void Em8RunAction::BeginOfRunAction(const G4Run* aRun)
67{ 
68  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
69 
70  // save Rndm status
71  if (saveRndm > 0)
72  { 
73      CLHEP::HepRandom::showEngineStatus();
74      CLHEP::HepRandom::saveEngineStatus("beginOfRun.rndm");
75  } 
76  G4UImanager* UI = G4UImanager::GetUIpointer();
77   
78  G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance();
79
80  if(pVVisManager)    UI->ApplyCommand("/vis/scene/notifyHandlers");
81
82     
83  EnergySumAbs = 0. ;
84  EnergySquareSumAbs = 0.;
85  tlSumAbs = 0. ;
86  tlsquareSumAbs = 0. ;
87  nStepSumCharged = 0. ;
88  nStepSum2Charged= 0. ;
89  nStepSumNeutral = 0. ;
90  nStepSum2Neutral= 0. ;
91  TotNbofEvents = 0. ;
92  SumCharged=0.;
93  SumNeutral=0.;
94  Sum2Charged=0.;
95  Sum2Neutral=0.;
96  Selectron=0.;
97  Spositron=0.;
98
99  Transmitted=0.;
100  Reflected  =0.;
101
102//  plot definitions
103   
104  if(nbinStep>0)
105  {
106    dStep=(Stephigh-Steplow)/nbinStep;
107    entryStep=0.;
108    underStep=0.;
109    overStep=0.;
110    for(G4int ist=0; ist<nbinStep; ist++)
111    {
112      distStep[ist]=0.;
113    }
114  }     
115  if(nbinEn>0)
116  {
117    dEn = (Enhigh-Enlow)/nbinEn ;
118    entryEn=0.;
119    underEn=0.;
120    overEn=0.;
121
122    for (G4int ien=0; ien<nbinEn; ien++)   distEn[ien]=0.;
123  }
124  if(nbinTt>0)
125  {
126    dTt = (Tthigh-Ttlow)/nbinTt ;
127    entryTt=0.;
128    underTt=0.;
129    overTt=0.;
130
131    for (G4int itt=0; itt<nbinTt; itt++)  distTt[itt]=0.;
132
133    Ttmean=0.;
134    Tt2mean=0.;
135  }
136  if(nbinTb>0)
137  {
138    dTb = (Tbhigh-Tblow)/nbinTb ;
139    entryTb=0.;
140    underTb=0.;
141    overTb=0.;
142    for (G4int itt=0; itt<nbinTb; itt++)
143    {
144      distTb[itt]=0.;
145    }
146    Tbmean=0.;
147    Tb2mean=0.;
148  }
149  if(nbinTsec>0)
150  {
151    dTsec = (Tsechigh-Tseclow)/nbinTsec ;
152    entryTsec=0.;
153    underTsec=0.;
154    overTsec=0.;
155    for (G4int its=0; its<nbinTsec; its++)
156    {
157      distTsec[its]=0.;
158    }
159  }
160  if(nbinTh>0)
161  {
162    dTh = (Thhigh-Thlow)/nbinTh ;
163    entryTh=0.;
164    underTh=0.;
165    overTh=0.;
166    for (G4int ith=0; ith<nbinTh; ith++)
167    {
168      distTh[ith]=0.;
169    }
170  }
171
172  if(nbinThback>0)
173  {
174    dThback = (Thhighback-Thlowback)/nbinThback ;
175    entryThback=0.;
176    underThback=0.;
177    overThback=0.;
178    for (G4int ithback=0; ithback<nbinThback; ithback++)
179    {
180      distThback[ithback]=0.;
181    }
182  }
183
184
185  if(nbinR >0)
186  {
187    dR  = (Rhigh-Rlow)/nbinR  ;
188    entryR =0.;
189    underR =0.;
190    overR =0.;
191    for (G4int ir =0; ir <nbinR ; ir++)
192    {
193      distR[ir]=0.;
194    }
195    Rmean=0.;
196    R2mean=0.;
197  }
198
199  if(nbinGamma>0)
200  {
201    dEGamma = std::log(EhighGamma/ElowGamma)/nbinGamma ;
202    entryGamma = 0.;
203    underGamma=0.;
204    overGamma=0.;
205    for (G4int ig=0; ig<nbinGamma; ig++)
206    {
207      distGamma[ig]=0.;
208    }
209  } 
210  if(nbinvertexz>0)
211  {
212    dz=(zhigh-zlow)/nbinvertexz;
213    entryvertexz=0.;
214    undervertexz=0.;
215    oververtexz=0.;
216    for(G4int iz=0; iz<nbinvertexz; iz++)
217    {
218      distvertexz[iz]=0.;
219    }
220  }
221}
222
223/////////////////////////////////////////////////////////////////////////////
224
225void Em8RunAction::EndOfRunAction(const G4Run*)
226{
227  G4double sAbs,sigAbs,sigstep,sigcharged,signeutral;
228
229  tlSumAbs /= TotNbofEvents ;
230  sAbs = tlsquareSumAbs/TotNbofEvents-tlSumAbs*tlSumAbs ;
231  if(sAbs>0.)
232    sAbs = std::sqrt(sAbs/TotNbofEvents) ;
233  else
234    sAbs = 0. ;
235 
236  EnergySumAbs /= TotNbofEvents ;
237  sigAbs = EnergySquareSumAbs/TotNbofEvents-EnergySumAbs*EnergySumAbs;
238  if(sigAbs>0.)
239    sigAbs = std::sqrt(sigAbs/TotNbofEvents);
240  else
241    sigAbs = 0.;
242
243  nStepSumCharged /= TotNbofEvents ;
244  sigstep = nStepSum2Charged/TotNbofEvents-nStepSumCharged*nStepSumCharged;
245  if(sigstep>0.)
246    sigstep = std::sqrt(sigstep/TotNbofEvents);
247  else
248    sigstep = 0.;
249  G4double sigch=sigstep ;
250 
251  nStepSumNeutral /= TotNbofEvents ;
252  sigstep = nStepSum2Neutral/TotNbofEvents-nStepSumNeutral*nStepSumNeutral;
253  if(sigstep>0.)
254    sigstep = std::sqrt(sigstep/TotNbofEvents);
255  else
256    sigstep = 0.;
257  G4double signe=sigstep ;
258 
259  SumCharged /= TotNbofEvents;
260  sigcharged = Sum2Charged/TotNbofEvents-SumCharged*SumCharged; 
261  if(sigcharged>0.)
262    sigcharged = std::sqrt(sigcharged/TotNbofEvents);
263  else
264    sigcharged = 0. ;
265 
266  SumNeutral /= TotNbofEvents;
267  signeutral = Sum2Neutral/TotNbofEvents-SumNeutral*SumNeutral; 
268  if(signeutral>0.)
269    signeutral = std::sqrt(signeutral/TotNbofEvents);
270  else
271    signeutral = 0. ;
272 
273  Selectron /= TotNbofEvents ;
274  Spositron /= TotNbofEvents ;
275
276  Transmitted /=TotNbofEvents ;
277  Reflected   /=TotNbofEvents ;
278  G4cout << " ================== run summary =====================" << G4endl;
279  G4int prec = G4cout.precision(6);
280  G4cout << " end of Run TotNbofEvents = " << 
281           TotNbofEvents << G4endl ;
282  G4cout << "    mean charged track length   in absorber=" <<
283           tlSumAbs/mm      << " +- " << sAbs/mm    <<
284          "  mm  " << G4endl; 
285  G4cout << G4endl;
286  G4cout << "            mean energy deposit in absorber=" <<
287           EnergySumAbs/MeV << " +- " << sigAbs/MeV <<
288          "  MeV " << G4endl ;
289  G4cout << G4endl ;
290  G4cout << " mean number of steps in absorber (charged) =" <<
291           nStepSumCharged         << " +- " << sigch     <<
292          "      " << G4endl ;
293  G4cout << " mean number of steps in absorber (neutral) =" <<
294           nStepSumNeutral         << " +- " << signe     <<
295          "      " << G4endl ;
296  G4cout << G4endl ;
297  G4cout << "   mean number of charged secondaries = " <<
298           SumCharged << " +- " << sigcharged << G4endl; 
299  G4cout << G4endl ;
300  G4cout << "   mean number of neutral secondaries = " <<
301           SumNeutral << " +- " << signeutral << G4endl; 
302  G4cout << G4endl ;
303 
304  G4cout << "   mean number of e-s =" << Selectron << 
305            "  and e+s =" << Spositron << G4endl;
306  G4cout << G4endl; 
307 
308  G4cout << "(number) transmission coeff=" << Transmitted <<
309            "  reflection coeff=" << Reflected << G4endl;
310  G4cout << G4endl; 
311
312  if(nbinStep>0)
313  {G4double E , dnorm, norm ;
314   G4cout << "   step number/event distribution " << G4endl ;
315   G4cout << "#entries=" << entryStep << "    #underflows=" << underStep <<
316             "    #overflows=" << overStep << G4endl ;
317   if( entryStep>0.)
318   {
319     E = Steplow - dStep ;
320     norm = TotNbofEvents ;
321     G4cout << " bin nb   nsteplow     entries     normalized " << G4endl ;
322     for(G4int iss=0; iss<nbinStep; iss++)
323     {
324      E += dStep ;
325      dnorm = distStep[iss]/norm;
326      G4cout << std::setw(5) << iss << std::setw(10) << E << 
327                std::setw(12) << distStep[iss] <<
328                std::setw(12) << dnorm << G4endl ;
329     }
330     G4cout << G4endl;
331   }     
332  }
333  if(nbinEn > 0)
334  {
335    std::ofstream fileOut("distribution.out", std::ios::out ) ;
336    fileOut.setf( std::ios::scientific, std::ios::floatfield );
337
338    std::ofstream normOut("normDist.out", std::ios::out ) ;
339    normOut.setf( std::ios::scientific, std::ios::floatfield );
340
341    G4double E , dnorm, norm,fmax,Emp; // width ;
342    Emp=-999.999 ;
343    G4cout << " energy deposit distribution " << G4endl ;
344    G4cout << "#entries=" << entryEn << "    #underflows=" << underEn <<
345             "    #overflows=" << overEn << G4endl ;
346    if( entryEn>0.)
347    {
348      E = Enlow - dEn ;
349      norm = TotNbofEvents*1.0 ;   // *dEn ;
350      G4cout << " bin nb      Elow      entries     normalized " << G4endl ;
351      fmax = 0. ;
352
353      for(G4int ien=0; ien<nbinEn; ien++)
354      {
355        E += dEn ;
356
357        if(distEn[ien]>fmax)
358        {
359          fmax = distEn[ien] ;
360          Emp = E ;                // most probable roughly
361        }
362        dnorm = distEn[ien]/norm;
363
364        G4cout << std::setw(5) << ien << std::setw(10) << E/keV << 
365                  std::setw(12) << distEn[ien] <<
366                  std::setw(12) << dnorm << G4endl ;
367
368        fileOut << E/keV << "\t"<< distEn[ien] << G4endl ;
369        normOut << E/keV << "\t"<< dnorm << G4endl ;
370      }
371      G4cout << G4endl;
372      G4int ii ;
373      G4double E1,E2 ;
374      E1=-1.e6 ;
375      E2=+1.e6 ;
376      E = Enlow -dEn ;
377      ii = -1;
378
379      for(G4int i1=0; i1<nbinEn; i1++)
380      {
381        E += dEn ;
382        if(ii<0)
383        {
384          if(distEn[i1] >= 0.5*fmax)
385          {
386            E1=E ;
387            ii=i1 ;
388          }
389        }
390      }
391      E = Enlow -dEn ;
392
393      for(G4int i2=0; i2<nbinEn; i2++)
394      {
395        E += dEn ;
396
397        if(distEn[i2] >= 0.5*fmax)   E2=E ;
398      }
399      G4cout << " Emp = " << std::setw(15) << Emp/MeV << "   width="
400            << std::setw(15) << (E2-E1)/MeV <<   "  MeV " << G4endl;
401      G4cout << G4endl ;
402    }     
403  }
404  if(nbinTt>0)
405  {
406     G4double E , dnorm, norm ,sig;
407   G4cout << " transmitted energy distribution " << G4endl ;
408   G4cout << "#entries=" << entryTt << "    #underflows=" << underTt <<
409             "    #overflows=" << overTt << G4endl ;
410   if( entryTt>0.)
411   {
412     Ttmean /= entryTt;
413     sig=Tt2mean/entryTt-Ttmean*Ttmean ;
414     if(sig<=0.)
415       sig=0.;
416     else
417       sig=std::sqrt(sig/entryTt) ;
418     G4cout << " mean energy of transmitted particles=" << Ttmean/keV << 
419               " +- " << sig/keV << "  keV." << G4endl;
420     E = Ttlow - dTt ;
421     norm = TotNbofEvents*dTt ;
422     G4cout << " bin nb      Elow      entries     normalized " << G4endl ;
423     for(G4int itt=0; itt<nbinTt; itt++)
424     {
425      E += dTt ;
426      dnorm = distTt[itt]/norm;
427      G4cout << std::setw(5) << itt << std::setw(10) << E << 
428                std::setw(12) << distTt[itt] <<
429                std::setw(12) << dnorm << G4endl ;
430     }
431     G4cout << G4endl;
432   }     
433  }
434  if(nbinTb>0)
435  {
436     G4double E , dnorm, norm ,sig;
437   G4cout << " backscattered energy distribution " << G4endl ;
438   G4cout << "#entries=" << entryTb << "    #underflows=" << underTb <<
439             "    #overflows=" << overTb << G4endl ;
440   if( entryTb>0.)
441   {
442     Tbmean /= entryTb;
443     sig=Tb2mean/entryTb-Tbmean*Tbmean ;
444     if(sig<=0.)
445       sig=0.;
446     else
447       sig=std::sqrt(sig/entryTb) ;
448     G4cout << " mean energy of backscattered particles=" << Tbmean/keV << 
449               " +- " << sig/keV << "  keV." << G4endl;
450     E = Tblow - dTb ;
451     norm = TotNbofEvents*dTb ;
452     G4cout << " bin nb      Elow      entries     normalized " << G4endl ;
453     for(G4int itt=0; itt<nbinTb; itt++)
454     {
455      E += dTb ;
456      dnorm = distTb[itt]/norm;
457      G4cout << std::setw(5) << itt << std::setw(10) << E << 
458                std::setw(12) << distTb[itt] <<
459                std::setw(12) << dnorm << G4endl ;
460     }
461     G4cout << G4endl;
462   }     
463  }
464  if(nbinTsec>0)
465  {G4double E , dnorm, norm ;
466   G4cout << " energy distribution of charged secondaries " << G4endl ;
467   G4cout << "#entries=" << entryTsec << "    #underflows=" << underTsec <<
468             "    #overflows=" << overTsec << G4endl ;
469   if( entryTsec>0.)
470   {
471     E = Tseclow - dTsec ;
472     norm = TotNbofEvents*dTsec ;
473     G4cout << " bin nb      Elow      entries     normalized " << G4endl ;
474     for(G4int itt=0; itt<nbinTsec; itt++)
475     {
476      E += dTsec ;
477      dnorm = distTsec[itt]/norm;
478      G4cout << std::setw(5) << itt << std::setw(10) << E << 
479                std::setw(12) << distTsec[itt] <<
480                std::setw(12) << dnorm << G4endl ;
481     }
482     G4cout << G4endl;
483   }     
484  }
485
486  if(nbinR >0)
487  {G4double R , dnorm, norm,sig  ;
488   G4cout << "  R  distribution " << G4endl ;
489   G4cout << "#entries=" << entryR  << "    #underflows=" << underR  <<
490             "    #overflows=" << overR  << G4endl ;
491   if( entryR >0.)
492   {
493     Rmean /= entryR;
494     sig = R2mean/entryR - Rmean*Rmean;
495     if(sig<=0.) sig=0. ;
496     else        sig = std::sqrt(sig/entryR) ;
497     G4cout << " mean lateral displacement at exit=" << Rmean/mm << " +- "
498            << sig/mm << "  mm." << G4endl ; 
499     R = Rlow - dR  ;
500     norm = TotNbofEvents*dR  ;
501     G4cout << " bin nb      Rlow      entries     normalized " << G4endl ;
502     for(G4int ier=0; ier<nbinR ; ier++)
503     {
504      R+= dR  ;
505      dnorm = distR[ier]/norm;
506      G4cout << std::setw(5) << ier << std::setw(10) << R  <<
507                std::setw(12) << distR[ier] <<
508                std::setw(12) << dnorm << G4endl ;
509     }
510     G4cout << G4endl;
511   }
512  }
513
514  if(nbinTh>0)
515  {G4double Th,Thdeg, dnorm, norm,fac0,fnorm,pere,Thpere,Thmean,sum;
516   G4cout << "      angle   distribution " << G4endl ;
517   G4cout << "#entries=" << entryTh << "    #underflows=" << underTh <<
518             "    #overflows=" << overTh << G4endl ;
519   if( entryTh>0.)
520   {
521     Th= Thlow - dTh ;
522     norm = TotNbofEvents ;
523     if(distTh[0] == 0.)
524       fac0 = 1. ;
525     else
526       fac0 = 1./distTh[0] ;
527     pere = 1./std::exp(1.) ;
528
529     G4cout << " bin nb  Thlowdeg      Thlowrad      " <<
530               " entries         normalized " << G4endl ;
531     Thpere = 0. ;
532     sum = 0. ;
533     Thmean = 0. ;
534     for(G4int ien=0; ien<nbinTh; ien++)
535     {
536      Th+= dTh ;
537      Thdeg = Th*180./pi ;
538      sum += distTh[ien] ;
539      Thmean += distTh[ien]*(Th+0.5*dTh) ;
540      dnorm = distTh[ien]/norm;
541      fnorm = fac0*distTh[ien] ;
542      if( fnorm > pere)
543        Thpere = Th ; 
544      G4cout << std::setw(5) << ien << std::setw(10) << Thdeg << "   " <<
545                std::setw(10) << Th << "  " <<   
546                std::setw(12) << distTh[ien] << "  " <<
547                std::setw(12) << dnorm << "  " << std::setw(12) << fnorm <<G4endl ;
548     }
549     Thmean /= sum ;
550     G4cout << G4endl;
551     G4cout << " mean = " << Thmean << "  rad  or " << 180.*Thmean/pi <<
552               " deg." << G4endl;
553     G4cout << " theta(1/e)=" << Thpere << " - " << Thpere+dTh << " rad   "
554            << " or " << 180.*Thpere/pi << " - " << 180.*(Thpere+dTh)/pi
555            << " deg." << G4endl;
556     G4cout << G4endl;
557   }
558  }
559
560  if(nbinThback>0)
561  {G4double Thb,Thdegb, dnormb, normb,fac0b,fnormb,pereb,Thpereb,Thmeanb,sumb;
562   G4cout << " backscattering angle   distribution " << G4endl ;
563   G4cout << "#entries=" << entryThback << "    #underflows=" << underThback <<
564             "    #overflows=" << overThback << G4endl ;
565   if( entryThback>0.)
566   {
567     Thb= Thlowback - dThback ;
568     normb = TotNbofEvents ;
569     if(distThback[0] == 0.)
570       fac0b = 1. ;
571     else
572       fac0b = 1./distThback[0] ;
573     pereb = 1./std::exp(1.) ;
574
575     G4cout << " bin nb  Thlowdeg      Thlowrad      " <<
576               " entries         normalized " << G4endl ;
577     Thpereb = 0. ;
578     sumb = 0. ;
579     Thmeanb = 0. ;
580     for(G4int ien=0; ien<nbinThback; ien++)
581     {
582      Thb+= dThback ;
583      Thdegb = Thb*180./pi ;
584      sumb += distThback[ien] ;
585      Thmeanb += distThback[ien]*(Thb+0.5*dThback) ;
586      dnormb = distThback[ien]/normb;
587      fnormb = fac0b*distThback[ien] ;
588      if( fnormb > pereb)
589        Thpereb = Thb ;
590      G4cout << std::setw(5) << ien << std::setw(10) << Thdegb << "   " <<
591                std::setw(10) << Thb << "  " <<
592                std::setw(12) << distThback[ien] << "  " <<
593                std::setw(12) << dnormb << "  " << std::setw(12) << fnormb <<G4endl ;
594     }
595     Thmeanb /= sumb ;
596     G4cout << G4endl;
597     G4cout << " mean = " << Thmeanb << "  rad  or " << 180.*Thmeanb/pi <<
598               " deg." << G4endl;
599     G4cout << " theta(1/e)=" << Thpereb << " - " << Thpereb+dThback << " rad   "
600            << " or " << 180.*Thpereb/pi << " - " << 180.*(Thpereb+dThback)/pi
601            << " deg." << G4endl;
602     G4cout << G4endl;
603   }
604  }
605
606  if(nbinGamma>0)
607  {G4double E , fact,dnorm, norm  ;
608   G4cout << " gamma energy distribution " << G4endl ;
609   G4cout << "#entries=" << entryGamma << "    #underflows=" << underGamma <<
610             "    #overflows=" << overGamma << G4endl ;
611   if( entryGamma>0.)
612   {
613     fact=std::exp(dEGamma) ;
614     E = ElowGamma/fact  ;
615     norm = TotNbofEvents*dEGamma;
616     G4cout << " bin nb         Elow      entries       normalized " << G4endl ;
617     for(G4int itt=0; itt<nbinGamma; itt++)
618     {
619      E *= fact ;
620      dnorm = distGamma[itt]/norm;
621      G4cout << std::setw(5) << itt << std::setw(13) << E << 
622                std::setw(12) << distGamma[itt] <<
623                std::setw(15) << dnorm << G4endl ;
624     }
625     G4cout << G4endl;
626   }     
627  }
628
629  if(nbinvertexz >0)
630  {G4double z , dnorm, norm  ;
631   G4cout << " vertex Z  distribution " << G4endl ;
632   G4cout << "#entries=" << entryvertexz  << "    #underflows=" << undervertexz  <<
633             "    #overflows=" << oververtexz  << G4endl ;
634   if( entryvertexz >0.)
635   {
636     z =zlow - dz  ;
637     norm = TotNbofEvents*dz  ;
638     G4cout << " bin nb      zlow      entries     normalized " << G4endl ;
639     for(G4int iez=0; iez<nbinvertexz ; iez++)
640     {
641      z+= dz  ;
642      if(std::fabs(z)<1.e-12) z=0.;
643      dnorm = distvertexz[iez]/norm;
644      G4cout << std::setw(5) << iez << std::setw(10) << z  <<
645                std::setw(12) << distvertexz[iez] <<
646                std::setw(12) << dnorm << G4endl ;
647     }
648     G4cout << G4endl;
649   }
650  }
651 
652 G4cout.precision(prec);
653 
654  if (G4VVisManager::GetConcreteInstance())
655  {
656    G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/update");
657  }
658
659  // save Rndm status
660
661  if (saveRndm == 1)
662  { 
663    CLHEP::HepRandom::showEngineStatus();
664    CLHEP::HepRandom::saveEngineStatus("endOfRun.rndm");
665  }     
666}
667
668///////////////////////////////////////////////////////////////////////////
669
670void Em8RunAction::CountEvent()
671{
672  TotNbofEvents += 1. ;
673}
674
675/////////////////////////////////////////////////////////////////////////
676
677void Em8RunAction::AddnStepsCharged(G4double ns)
678{
679  nStepSumCharged += ns;
680  nStepSum2Charged += ns*ns;
681}
682
683////////////////////////////////////////////////////////////////////////
684
685void Em8RunAction::AddnStepsNeutral(G4double ns)
686{
687  nStepSumNeutral += ns;
688  nStepSum2Neutral += ns*ns;
689}
690
691////////////////////////////////////////////////////////////////////////////
692
693void Em8RunAction::AddEdeps(G4double Eabs)
694{
695  EnergySumAbs += Eabs;
696  EnergySquareSumAbs += Eabs*Eabs;
697}
698
699/////////////////////////////////////////////////////////////////////////////
700
701void Em8RunAction::AddTrackLength(G4double tlabs)
702{
703  tlSumAbs += tlabs;
704  tlsquareSumAbs += tlabs*tlabs ;
705}
706
707/////////////////////////////////////////////////////////////////////////////
708
709void Em8RunAction::AddTrRef(G4double tr,G4double ref)
710{
711  Transmitted += tr ;
712  Reflected   += ref;
713}
714
715/////////////////////////////////////////////////////////////////////////////
716
717void Em8RunAction::FillNbOfSteps(G4double)
718{}
719
720//////////////////////////////////////////////////////////////////////////////
721
722void Em8RunAction::FillEn(G4double En)
723{
724
725  G4double bin ;
726  G4int ibin;
727
728  {
729    entryEn += 1. ;
730 
731    if(En < Enlow)          underEn += 1. ;
732    else if( En >= Enhigh)  overEn  += 1. ;
733    else
734    {
735      bin = (En-Enlow)/dEn ;
736      ibin= (G4int)bin ;
737      distEn[ibin] += 1. ;
738    }
739  }
740}
741
742////////////////////////////////////////////////////////////////////////////
743
744void Em8RunAction::FillTt(G4double)
745{
746}
747
748//////////////////////////////////////////////////////////////////////////////
749
750void Em8RunAction::FillTb(G4double)
751{
752}
753
754///////////////////////////////////////////////////////////////////////////////
755
756void Em8RunAction::FillTsec(G4double)
757{
758}
759
760/////////////////////////////////////////////////////////////////////////////
761
762void Em8RunAction::FillGammaSpectrum(G4double)
763{
764}
765
766////////////////////////////////////////////////////////////////////////////////
767
768void Em8RunAction::FillTh(G4double)
769{
770}
771
772//////////////////////////////////////////////////////////////////////////
773
774void Em8RunAction::FillThBack(G4double)
775{
776}
777
778//////////////////////////////////////////////////////////////////////////
779
780void Em8RunAction::FillR(G4double)
781{
782}
783
784/////////////////////////////////////////////////////////////////////////////
785
786void Em8RunAction::Fillvertexz(G4double)
787{
788}
789
790//////////////////////////////////////////////////////////////////////////////
791
792void Em8RunAction::SethistName(G4String name)
793{
794  histName = name ;
795  G4cout << " hist file = " << histName << G4endl;
796}
797
798void Em8RunAction::SetnbinStep(G4int nbin)
799{
800  nbinStep = nbin ;
801  if(nbinStep>0)
802  G4cout << " Nb of bins in #step plot = " << nbinStep << G4endl ;
803}
804
805void Em8RunAction::SetSteplow(G4double low)
806{
807  Steplow = low ;
808  if(nbinStep>0)
809  G4cout << " low  in the #step plot = " << Steplow << G4endl ;
810}
811void Em8RunAction::SetStephigh(G4double high)
812{
813  Stephigh = high ;
814  if(nbinStep>0)
815  G4cout << " high in the #step plot = " << Stephigh << G4endl ;
816}
817
818////////////////////////////////////////////////////////////////////////
819
820void Em8RunAction::SetnbinEn(G4int nbin)
821{
822  nbinEn = nbin ;
823
824  if(nbinEn > 0) G4cout << " Nb of bins in Edep plot = " << nbinEn << G4endl ;
825}
826
827void Em8RunAction::SetEnlow(G4double Elow)
828{
829  Enlow = Elow ;
830  if(nbinEn>0) G4cout << " Elow  in the  Edep plot = " << Enlow << G4endl ;
831}
832
833void Em8RunAction::SetEnhigh(G4double Ehigh)
834{
835  Enhigh = Ehigh ;
836
837  if(nbinEn>0) G4cout << " Ehigh in the  Edep plot = " << Enhigh << G4endl ;
838}
839
840/////////////////////////////////////////////////////////////////////////
841
842void Em8RunAction::SetnbinGamma(G4int nbin)
843{
844  nbinGamma = nbin ;
845  if(nbinGamma>0)
846  G4cout << " Nb of bins in gamma spectrum plot = " << nbinGamma << G4endl ;
847}
848
849void Em8RunAction::SetElowGamma(G4double Elow)
850{
851  ElowGamma = Elow ;
852  if(nbinGamma>0)
853  G4cout << " Elow  in the gamma spectrum plot = " << ElowGamma << G4endl ;
854}
855
856void Em8RunAction::SetEhighGamma(G4double Ehigh)
857{
858  EhighGamma = Ehigh ;
859  if(nbinGamma>0)
860  G4cout << " Ehigh in the gamma spectrum plot = " << EhighGamma << G4endl ;
861}
862
863void Em8RunAction::SetnbinTt(G4int nbin)
864{
865  nbinTt = nbin ;
866  if(nbinTt>0)
867  G4cout << " Nb of bins in Etransmisssion plot = " << nbinTt << G4endl ;
868}
869
870void Em8RunAction::SetTtlow(G4double Elow)
871{
872  Ttlow = Elow ;
873  if(nbinTt>0)
874  G4cout << " Elow  in the  Etransmission plot = " << Ttlow << G4endl ;
875}
876
877void Em8RunAction::SetTthigh(G4double Ehigh)
878{
879  Tthigh = Ehigh ;
880  if(nbinTt>0)
881  G4cout << " Ehigh in the  Etransmission plot = " << Tthigh << G4endl ;
882}
883
884void Em8RunAction::SetnbinTb(G4int nbin)
885{
886  nbinTb = nbin ;
887  if(nbinTb>0)
888  G4cout << " Nb of bins in Ebackscattered plot = " << nbinTb << G4endl ;
889}
890
891void Em8RunAction::SetTblow(G4double Elow)
892{
893  Tblow = Elow ;
894  if(nbinTb>0)
895  G4cout << " Elow  in the  Ebackscattered plot = " << Tblow << G4endl ;
896}
897
898void Em8RunAction::SetTbhigh(G4double Ehigh)
899{
900  Tbhigh = Ehigh ;
901  if(nbinTb>0)
902  G4cout << " Ehigh in the  Ebackscattered plot = " << Tbhigh << G4endl ;
903}
904
905void Em8RunAction::SetnbinTsec(G4int nbin)
906{
907  nbinTsec = nbin ;
908  if(nbinTsec>0)
909  G4cout << " Nb of bins in Tsecondary  plot = " << nbinTsec << G4endl ;
910}
911
912void Em8RunAction::SetTseclow(G4double Elow)
913{
914  Tseclow = Elow ;
915  if(nbinTsec>0)
916  G4cout << " Elow  in the  Tsecondary plot = " << Tseclow << G4endl ;
917}
918
919void Em8RunAction::SetTsechigh(G4double Ehigh)
920{
921  Tsechigh = Ehigh ;
922  if(nbinTsec>0)
923  G4cout << " Ehigh in the  Tsecondary plot = " << Tsechigh << G4endl ;
924}
925 
926void Em8RunAction::SetnbinR(G4int nbin)
927{
928  nbinR  = nbin ;
929  if(nbinR>0)
930  G4cout << " Nb of bins in R plot = " << nbinR << G4endl ;
931}
932
933void Em8RunAction::SetRlow(G4double rlow)
934{
935  Rlow = rlow ;
936  if(nbinR>0)
937  G4cout << " Rlow  in the  R plot = " << Rlow << G4endl ;
938}
939
940void Em8RunAction::SetRhigh(G4double rhigh)
941{
942  Rhigh = rhigh ;
943  if(nbinR>0)
944  G4cout << " Rhigh in the R plot = " << Rhigh << G4endl ;
945}
946
947void Em8RunAction::Setnbinzvertex(G4int nbin)
948{
949  nbinvertexz  = nbin ;
950  if(nbinvertexz>0)
951  G4cout << " Nb of bins in Z plot = " << nbinvertexz << G4endl ;
952}
953
954void Em8RunAction::Setzlow(G4double z)
955{
956  zlow = z ;
957  if(nbinvertexz>0)
958  G4cout << " zlow  in the  Z plot = " << zlow << G4endl ;
959}
960
961void Em8RunAction::Setzhigh(G4double z)
962{
963  zhigh = z ;
964  if(nbinvertexz>0)
965  G4cout << " zhigh in the Z plot = " << zhigh << G4endl ;
966}
967
968void Em8RunAction::SetnbinTh(G4int nbin)
969{
970  nbinTh = nbin ;
971  if(nbinTh>0)
972  G4cout << " Nb of bins in Theta plot = " << nbinTh << G4endl ;
973}
974
975void Em8RunAction::SetThlow(G4double Tlow)
976{
977  Thlow = Tlow ;
978  if(nbinTh>0)
979  G4cout << " Tlow  in the  Theta plot = " << Thlow << G4endl ;
980}
981
982void Em8RunAction::SetThhigh(G4double Thigh)
983{
984  Thhigh = Thigh ;
985  if(nbinTh>0)
986  G4cout << " Thigh in the Theta plot = " << Thhigh << G4endl ;
987}
988
989void Em8RunAction::SetnbinThBack(G4int nbin)
990{
991  nbinThback = nbin ;
992  if(nbinThback>0)
993  G4cout << " Nb of bins in Theta plot = " << nbinThback << G4endl ;
994}
995
996void Em8RunAction::SetThlowBack(G4double Tlow)
997{
998  Thlowback = Tlow ;
999  if(nbinThback>0)
1000  G4cout << " Tlow  in the  Theta plot = " << Thlowback << G4endl ;
1001}
1002
1003void Em8RunAction::SetThhighBack(G4double Thigh)
1004{
1005  Thhighback = Thigh ;
1006  if(nbinThback>0)
1007  G4cout << " Thigh in the Theta plot = " << Thhighback << G4endl ;
1008}
1009
1010void Em8RunAction::CountParticles(G4double nch,G4double nne)
1011{
1012  SumCharged += nch ;
1013  SumNeutral += nne ;
1014  Sum2Charged += nch*nch ;
1015  Sum2Neutral += nne*nne ;
1016}
1017
1018void Em8RunAction::AddEP(G4double nele,G4double npos) 
1019{
1020  Selectron += nele;
1021  Spositron += npos;
1022}
1023
1024//
1025//
1026////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.