source: trunk/source/geometry/magneticfield/test/field02/src/F02RunAction.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

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