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

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

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