Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm2/src/RunAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/extended/electromagnetic/TestEm2/src/RunAction.cc (Version 11.3.0) and /examples/extended/electromagnetic/TestEm2/src/RunAction.cc (Version 9.1.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 /// \file electromagnetic/TestEm2/src/RunActio <<  26 // $Id: RunAction.cc,v 1.21 2006/06/29 16:50:47 gunter Exp $
 27 /// \brief Implementation of the RunAction cla <<  27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
 28 //                                             <<  28 // 
 29 //                                             << 
 30 //....oooOO0OOooo........oooOO0OOooo........oo     29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oo     30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32                                                    31 
 33 #include "RunAction.hh"                            32 #include "RunAction.hh"
 34                                                    33 
 35 #include "DetectorConstruction.hh"                 34 #include "DetectorConstruction.hh"
 36 #include "EmAcceptance.hh"                     << 
 37 #include "PrimaryGeneratorAction.hh"               35 #include "PrimaryGeneratorAction.hh"
 38 #include "Run.hh"                              << 
 39 #include "RunActionMessenger.hh"                   36 #include "RunActionMessenger.hh"
                                                   >>  37 #include "EmAcceptance.hh"
 40                                                    38 
 41 #include "G4Run.hh"                                39 #include "G4Run.hh"
 42 #include "G4RunManager.hh"                         40 #include "G4RunManager.hh"
 43 #include "G4SystemOfUnits.hh"                  <<  41 #include "G4UnitsTable.hh"
                                                   >>  42 
                                                   >>  43 #include <iomanip>
                                                   >>  44 
 44 #include "Randomize.hh"                            45 #include "Randomize.hh"
 45                                                    46 
                                                   >>  47 #ifdef G4ANALYSIS_USE
                                                   >>  48 #include "AIDA/AIDA.h"
                                                   >>  49 #endif
                                                   >>  50 
 46 //....oooOO0OOooo........oooOO0OOooo........oo     51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 47                                                    52 
 48 RunAction::RunAction(DetectorConstruction* det <<  53 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin)
                                                   >>  54 :Det(det),Kin(kin),af(0),tree(0)
 49 {                                                  55 {
 50   fRunMessenger = new RunActionMessenger(this) <<  56   runMessenger = new RunActionMessenger(this);
 51                                                <<  57   
 52   // Create analysis manager                   <<  58   nLbin = nRbin = MaxBin;
 53   // The choice of analysis technology is done <<  59     
 54   fAnalysisManager = G4AnalysisManager::Instan <<  60   dEdL.resize(nLbin, 0.0);
 55   fAnalysisManager->SetDefaultFileType("root") <<  61   sumELongit.resize(nLbin, 0.0);
 56                                                <<  62   sumELongitCumul.resize(nLbin, 0.0);
 57   // Set the default file name "testem2"       <<  63   sumE2Longit.resize(nLbin, 0.0);
 58   // which can be then redefine in a macro via <<  64   sumE2LongitCumul.resize(nLbin, 0.0);
 59   fAnalysisManager->SetFileName("testem2");    <<  65 
                                                   >>  66   dEdR.resize(nRbin, 0.0);
                                                   >>  67   sumERadial.resize(nRbin, 0.0);
                                                   >>  68   sumERadialCumul.resize(nRbin, 0.0);
                                                   >>  69   sumE2Radial.resize(nRbin, 0.0);
                                                   >>  70   sumE2RadialCumul.resize(nRbin, 0.0);
                                                   >>  71     
                                                   >>  72   edeptrue = 1.;
                                                   >>  73   rmstrue  = 1.;
                                                   >>  74   limittrue = DBL_MAX;
                                                   >>  75   
                                                   >>  76   verbose = 0;
                                                   >>  77   
                                                   >>  78 #ifdef G4ANALYSIS_USE
                                                   >>  79   // Creating the analysis factory
                                                   >>  80   af = AIDA_createAnalysisFactory();
                                                   >>  81   if(!af) {
                                                   >>  82     G4cout << " RunAction::RunAction() :" 
                                                   >>  83            << " problem creating the AIDA analysis factory."
                                                   >>  84            << G4endl;
                                                   >>  85   }      
                                                   >>  86 #endif
                                                   >>  87     
                                                   >>  88   histoName[0] = "testem2";
                                                   >>  89   histoType    = "hbook";  
 60 }                                                  90 }
 61                                                    91 
 62 //....oooOO0OOooo........oooOO0OOooo........oo     92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 63                                                    93 
 64 RunAction::~RunAction()                            94 RunAction::~RunAction()
 65 {                                                  95 {
 66   delete fRunMessenger;                        <<  96   delete runMessenger;
                                                   >>  97   
                                                   >>  98 #ifdef G4ANALYSIS_USE
                                                   >>  99   delete af;
                                                   >> 100 #endif  
 67 }                                                 101 }
 68                                                   102 
 69 //....oooOO0OOooo........oooOO0OOooo........oo    103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 70                                                   104 
 71 void RunAction::BookHisto()                    << 105 void RunAction::bookHisto()
 72 {                                                 106 {
 73   // Get analysis manager                      << 107 #ifdef G4ANALYSIS_USE
 74   fAnalysisManager = G4AnalysisManager::Instan << 108   if(!af) return;
 75                                                << 109   
 76   // Open an output file                       << 110   // Creating a tree mapped to an hbook file.
 77   fAnalysisManager->OpenFile();                << 111   histoName[1] = histoName[0] + "." + histoType;  
 78   fAnalysisManager->SetVerboseLevel(1);        << 112   G4bool readOnly  = false;
 79                                                << 113   G4bool createNew = true;
 80   // Creating histograms                       << 114   G4String options = "--noErrors uncompress";
 81   //                                           << 115   AIDA::ITreeFactory* tf  = af->createTreeFactory();  
 82   G4double Ekin = fKin->GetParticleGun()->GetP << 116   tree = tf->create(histoName[1], histoType, readOnly, createNew, options);
 83   G4int nLbin = fDet->GetnLtot();              << 117   delete tf;
 84   G4int nRbin = fDet->GetnRtot();              << 118   if(!tree) {
 85   G4double dLradl = fDet->GetdLradl();         << 119     G4cout << "RunAction::bookHisto() :" 
 86   G4double dRradl = fDet->GetdRradl();         << 120            << " problem creating the AIDA tree with "
 87                                                << 121            << " storeName = " << histoName[1]
 88   fAnalysisManager->SetFirstHistoId(1);        << 122            << " readOnly = "  << readOnly
 89   fAnalysisManager->CreateH1("h1", "total ener << 123            << " createNew = " << createNew
 90                                                << 124            << G4endl;
 91   fAnalysisManager->CreateH1("h2", "total char << 125     return;
 92                                                << 126   }
 93   fAnalysisManager->CreateH1("h3", "total neut << 127   // Creating a histogram factory, whose histograms will be handled by the tree
 94                                                << 128   AIDA::IHistogramFactory* hf = af->createHistogramFactory(*tree);
 95   fAnalysisManager->CreateH1("h4", "longit ene << 129   
 96                                                << 130   G4double Ekin = Kin->GetParticleGun()->GetParticleEnergy();
 97   fAnalysisManager->CreateP1("p4", "longit ene << 131   G4double dLradl = Det->GetdLradl();
 98                              0., 1000.);       << 132   G4double dRradl = Det->GetdRradl();
                                                   >> 133 
                                                   >> 134   histo[0] = hf->createHistogram1D( "1","total energy deposit(percent of Einc)",
                                                   >> 135                                     110,0.,110.);
                                                   >> 136 
                                                   >> 137   histo[1] = hf->createHistogram1D( "2","total charged tracklength (radl)",
                                                   >> 138                                     110,0.,110.*Ekin/GeV);
                                                   >> 139 
                                                   >> 140   histo[2] = hf->createHistogram1D( "3","total neutral tracklength (radl)",
                                                   >> 141                                     110,0.,1100.*Ekin/GeV);
                                                   >> 142 
                                                   >> 143   histo[3] = hf->createHistogram1D( "4","longit energy profile (% of E inc)",
                                                   >> 144                                     nLbin,0.,nLbin*dLradl);
                                                   >> 145 
                                                   >> 146   G4double Zmin=0.5*dLradl, Zmax=Zmin+nLbin*dLradl;
                                                   >> 147   histo[4] = hf->createHistogram1D( "5","cumul longit energy dep (% of E inc)",
                                                   >> 148                                     nLbin,Zmin,Zmax);
                                                   >> 149 
                                                   >> 150   histo[5] = hf->createHistogram1D( "6","rms on cumul longit Edep (% of E inc)",
                                                   >> 151                                     nLbin,Zmin,Zmax);
                                                   >> 152 
                                                   >> 153   histo[6] = hf->createHistogram1D("7","radial energy profile (% of E inc)",
                                                   >> 154                                     nRbin,0.,nRbin*dRradl);
                                                   >> 155 
                                                   >> 156   G4double Rmin=0.5*dRradl, Rmax=Rmin+nRbin*dRradl;
                                                   >> 157   histo[7]= hf->createHistogram1D("8","cumul radial energy dep (% of E inc)",
                                                   >> 158                                     nRbin,Rmin,Rmax);
                                                   >> 159 
                                                   >> 160   histo[8]= hf->createHistogram1D("9","rms on cumul radial Edep (% of E inc)",
                                                   >> 161                                     nRbin,Rmin,Rmax);
                                                   >> 162             
                                                   >> 163  delete hf;
                                                   >> 164  G4cout << "\n----> Histogram Tree is opened in " << histoName[1] << G4endl;
                                                   >> 165 #endif
                                                   >> 166 }
 99                                                   167 
100   fAnalysisManager->CreateH1("h5", "rms on lon << 168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101                                                   169 
102   G4double Zmin = 0.5 * dLradl, Zmax = Zmin +  << 170 void RunAction::cleanHisto()
103   fAnalysisManager->CreateH1("h6", "cumul long << 171 {
                                                   >> 172 #ifdef G4ANALYSIS_USE
                                                   >> 173   if(tree) {
                                                   >> 174     tree->commit();       // Writing the histograms to the file
                                                   >> 175     tree->close();        // and closing the tree (and the file)
                                                   >> 176     G4cout << "\n----> Histogram Tree is saved in " << histoName[1] << G4endl; 
                                                   >> 177    
                                                   >> 178     delete tree;
                                                   >> 179     tree = 0;
                                                   >> 180   }
                                                   >> 181 #endif
                                                   >> 182 }
104                                                   183 
105   fAnalysisManager->CreateH1("h7", "rms on cum << 184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
106                                                   185 
107   fAnalysisManager->CreateH1("h8", "radial ene << 186 void RunAction::BeginOfRunAction(const G4Run* aRun)
                                                   >> 187 {
                                                   >> 188   G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
108                                                   189 
109   fAnalysisManager->CreateP1("p8", "radial ene << 190   // save Rndm status
110                              0., 1000.);       << 191   G4RunManager::GetRunManager()->SetRandomNumberStore(true);
                                                   >> 192   CLHEP::HepRandom::showEngineStatus();
111                                                   193 
112   fAnalysisManager->CreateH1("h9", "rms on rad << 194   nLbin = Det->GetnLtot();
                                                   >> 195   nRbin = Det->GetnRtot();
113                                                   196 
114   G4double Rmin = 0.5 * dRradl, Rmax = Rmin +  << 197   //initialize arrays of cumulative energy deposition
115   fAnalysisManager->CreateH1("h10", "cumul rad << 198   //
                                                   >> 199   for (G4int i=0; i<nLbin; i++) 
                                                   >> 200      sumELongit[i]=sumE2Longit[i]=sumELongitCumul[i]=sumE2LongitCumul[i]=0.;
                                                   >> 201   
                                                   >> 202   for (G4int j=0; j<nRbin; j++)
                                                   >> 203      sumERadial[j]=sumE2Radial[j]=sumERadialCumul[j]=sumE2RadialCumul[j]=0.;
116                                                   204 
117   fAnalysisManager->CreateH1("h11", "rms on cu << 205   //initialize track length
                                                   >> 206   sumChargTrLength=sum2ChargTrLength=sumNeutrTrLength=sum2NeutrTrLength=0.;
118                                                   207 
119   G4String fileName = fAnalysisManager->GetFil << 208   //histograms
120   G4String extension = fAnalysisManager->GetFi << 209   //
121   G4String fullFileName = fileName + "." + ext << 210   if (aRun->GetRunID() == 0) bookHisto();
122   G4cout << "\n----> Histogram file is opened  << 
123 }                                                 211 }
124                                                   212 
125 //....oooOO0OOooo........oooOO0OOooo........oo    213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
126                                                   214 
127 G4Run* RunAction::GenerateRun()                << 215 void RunAction::fillPerEvent()
128 {                                                 216 {
129   fRun = new Run(fDet, fKin);                  << 217   //accumulate statistic
130   fRun->SetVerbose(fVerbose);                  << 218   //
131   return fRun;                                 << 219   G4double dLCumul = 0.;
                                                   >> 220   for (G4int i=0; i<nLbin; i++)
                                                   >> 221      {
                                                   >> 222       sumELongit[i]  += dEdL[i];
                                                   >> 223       sumE2Longit[i] += dEdL[i]*dEdL[i];
                                                   >> 224       dLCumul        += dEdL[i];
                                                   >> 225       sumELongitCumul[i]  += dLCumul;
                                                   >> 226       sumE2LongitCumul[i] += dLCumul*dLCumul;
                                                   >> 227      }
                                                   >> 228 
                                                   >> 229   G4double dRCumul = 0.;
                                                   >> 230   for (G4int j=0; j<nRbin; j++)
                                                   >> 231      {
                                                   >> 232       sumERadial[j]  += dEdR[j];
                                                   >> 233       sumE2Radial[j] += dEdR[j]*dEdR[j];
                                                   >> 234       dRCumul        += dEdR[j];
                                                   >> 235       sumERadialCumul[j]  += dRCumul;
                                                   >> 236       sumE2RadialCumul[j] += dRCumul*dRCumul;
                                                   >> 237      }
                                                   >> 238 
                                                   >> 239   sumChargTrLength  += ChargTrLength;
                                                   >> 240   sum2ChargTrLength += ChargTrLength*ChargTrLength;
                                                   >> 241   sumNeutrTrLength  += NeutrTrLength;
                                                   >> 242   sum2NeutrTrLength += NeutrTrLength*NeutrTrLength;
                                                   >> 243 
                                                   >> 244 #ifdef G4ANALYSIS_USE
                                                   >> 245   //fill histograms
                                                   >> 246   //
                                                   >> 247   if(tree) {
                                                   >> 248     G4double Ekin=Kin->GetParticleGun()->GetParticleEnergy();
                                                   >> 249     G4double mass=Kin->GetParticleGun()->GetParticleDefinition()->GetPDGMass();
                                                   >> 250     G4double radl=Det->GetMaterial()->GetRadlen();
                                                   >> 251 
                                                   >> 252     histo[0]->fill(100.*dLCumul/(Ekin+mass));
                                                   >> 253     histo[1]->fill(ChargTrLength/radl);
                                                   >> 254     histo[2]->fill(NeutrTrLength/radl);
                                                   >> 255   }
                                                   >> 256 #endif
132 }                                                 257 }
133                                                   258 
134 //....oooOO0OOooo........oooOO0OOooo........oo    259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135                                                   260 
136 void RunAction::BeginOfRunAction(const G4Run*) << 261 void RunAction::EndOfRunAction(const G4Run* aRun)
137 {                                                 262 {
138   // show Rndm status                          << 263   //compute and print statistic
139   if (isMaster) G4Random::showEngineStatus();  << 264   //
                                                   >> 265   G4int NbOfEvents = aRun->GetNumberOfEvent();
                                                   >> 266   G4double kinEnergy = Kin->GetParticleGun()->GetParticleEnergy();
                                                   >> 267   assert(NbOfEvents*kinEnergy > 0);
                                                   >> 268 
                                                   >> 269   G4double mass=Kin->GetParticleGun()->GetParticleDefinition()->GetPDGMass();
                                                   >> 270   G4double norme = 100./(NbOfEvents*(kinEnergy+mass));
140                                                   271 
141   // histograms                                << 272   //longitudinal
142   //                                              273   //
143   BookHisto();                                 << 274   G4double dLradl = Det->GetdLradl();
144 }                                              << 
145                                                   275 
146 //....oooOO0OOooo........oooOO0OOooo........oo << 276   MyVector MeanELongit(nLbin),      rmsELongit(nLbin);
                                                   >> 277   MyVector MeanELongitCumul(nLbin), rmsELongitCumul(nLbin);
147                                                   278 
148 void RunAction::EndOfRunAction(const G4Run*)   << 279   G4int i;
149 {                                              << 280   for (i=0; i<nLbin; i++)
150   // compute and print statistic               << 281    {
                                                   >> 282     MeanELongit[i] = norme*sumELongit[i];
                                                   >> 283      rmsELongit[i] = norme*std::sqrt(std::fabs(NbOfEvents*sumE2Longit[i]
                                                   >> 284                                 - sumELongit[i]*sumELongit[i]));
                                                   >> 285 
                                                   >> 286     MeanELongitCumul[i] = norme*sumELongitCumul[i];
                                                   >> 287      rmsELongitCumul[i] = norme*std::sqrt(std::fabs(NbOfEvents*sumE2LongitCumul[i]
                                                   >> 288                                     - sumELongitCumul[i]*sumELongitCumul[i]));
                                                   >> 289 
                                                   >> 290 
                                                   >> 291 #ifdef G4ANALYSIS_USE
                                                   >> 292     if(tree) {
                                                   >> 293       G4double bin = i*dLradl;
                                                   >> 294       histo[3]->fill(bin,MeanELongit[i]/dLradl);
                                                   >> 295       bin = (i+1)*dLradl;
                                                   >> 296       histo[4]->fill(bin,MeanELongitCumul[i]);
                                                   >> 297       histo[5]->fill(bin, rmsELongitCumul[i]);
                                                   >> 298     }
                                                   >> 299 #endif
                                                   >> 300    }
                                                   >> 301 
                                                   >> 302   //radial
151   //                                              303   //
152   if (isMaster) fRun->EndOfRun(fEdeptrue, fRms << 304   G4double dRradl = Det->GetdRradl();
                                                   >> 305 
                                                   >> 306   MyVector MeanERadial(nRbin),      rmsERadial(nRbin);
                                                   >> 307   MyVector MeanERadialCumul(nRbin), rmsERadialCumul(nRbin);
153                                                   308 
                                                   >> 309   for (i=0; i<nRbin; i++)
                                                   >> 310    {
                                                   >> 311     MeanERadial[i] = norme*sumERadial[i];
                                                   >> 312      rmsERadial[i] = norme*std::sqrt(std::fabs(NbOfEvents*sumE2Radial[i]
                                                   >> 313                                 - sumERadial[i]*sumERadial[i]));
                                                   >> 314 
                                                   >> 315     MeanERadialCumul[i] = norme*sumERadialCumul[i];
                                                   >> 316      rmsERadialCumul[i] = norme*std::sqrt(std::fabs(NbOfEvents*sumE2RadialCumul[i]
                                                   >> 317                                     - sumERadialCumul[i]*sumERadialCumul[i]));
                                                   >> 318 
                                                   >> 319 #ifdef G4ANALYSIS_USE
                                                   >> 320     if(tree) {
                                                   >> 321       G4double bin = i*dRradl;
                                                   >> 322       histo[6]->fill(bin,MeanERadial[i]/dRradl);
                                                   >> 323       bin = (i+1)*dRradl;
                                                   >> 324       histo[7]->fill(bin,MeanERadialCumul[i]);
                                                   >> 325       histo[8]->fill(bin, rmsERadialCumul[i]);
                                                   >> 326     }
                                                   >> 327 #endif
                                                   >> 328    }
                                                   >> 329 
                                                   >> 330   //find Moliere confinement
                                                   >> 331   //
                                                   >> 332   const G4double EMoliere = 90.;
                                                   >> 333   G4double iMoliere = 0.;
                                                   >> 334   if ((MeanERadialCumul[0]       <= EMoliere) &&
                                                   >> 335       (MeanERadialCumul[nRbin-1] >= EMoliere)) {
                                                   >> 336     G4int imin = 0;
                                                   >> 337     while( (imin < nRbin-1) && (MeanERadialCumul[imin] < EMoliere) ) imin++;
                                                   >> 338     G4double ratio = (EMoliere - MeanERadialCumul[imin]) /
                                                   >> 339                      (MeanERadialCumul[imin+1] - MeanERadialCumul[imin]);
                                                   >> 340     iMoliere = 1. + imin + ratio;
                                                   >> 341   }          
                                                   >> 342       
                                                   >> 343   //track length
                                                   >> 344   //
                                                   >> 345   norme = 1./(NbOfEvents*(Det->GetMaterial()->GetRadlen()));
                                                   >> 346   G4double MeanChargTrLength = norme*sumChargTrLength;
                                                   >> 347   G4double  rmsChargTrLength = 
                                                   >> 348             norme*std::sqrt(std::fabs(NbOfEvents*sum2ChargTrLength
                                                   >> 349                                          - sumChargTrLength*sumChargTrLength));
                                                   >> 350 
                                                   >> 351   G4double MeanNeutrTrLength = norme*sumNeutrTrLength;
                                                   >> 352   G4double  rmsNeutrTrLength = 
                                                   >> 353             norme*std::sqrt(std::fabs(NbOfEvents*sum2NeutrTrLength
                                                   >> 354                                          - sumNeutrTrLength*sumNeutrTrLength));
                                                   >> 355 
                                                   >> 356   //print
                                                   >> 357   //
                                                   >> 358 
                                                   >> 359   std::ios::fmtflags mode = G4cout.flags();
                                                   >> 360   G4cout.setf(std::ios::fixed,std::ios::floatfield);
                                                   >> 361   G4int  prec = G4cout.precision(2);
                                                   >> 362   
                                                   >> 363   if (verbose) {
                                                   >> 364 
                                                   >> 365     G4cout << "                 LATERAL PROFILE                   "
                                                   >> 366            << "      CUMULATIVE LATERAL PROFILE" << G4endl << G4endl;
                                                   >> 367 
                                                   >> 368     G4cout << "        bin   " << "           Mean         rms         "
                                                   >> 369            << "        bin "   << "           Mean      rms \n" << G4endl;
                                                   >> 370 
                                                   >> 371     for (i=0; i<nLbin; i++)
                                                   >> 372      {
                                                   >> 373        G4double inf=i*dLradl, sup=inf+dLradl;
                                                   >> 374 
                                                   >> 375        G4cout << std::setw(8) << inf << "->"
                                                   >> 376               << std::setw(5) << sup << " radl: "
                                                   >> 377               << std::setw(7) << MeanELongit[i] << "%  "
                                                   >> 378               << std::setw(9) << rmsELongit[i] << "%       "
                                                   >> 379               << "      0->" << std::setw(5) << sup << " radl: "
                                                   >> 380               << std::setw(7) << MeanELongitCumul[i] << "%  "
                                                   >> 381               << std::setw(7) << rmsELongitCumul[i] << "% "
                                                   >> 382               <<G4endl;
                                                   >> 383      }
                                                   >> 384 
                                                   >> 385     G4cout << G4endl << G4endl << G4endl;
                                                   >> 386 
                                                   >> 387     G4cout << "                  RADIAL PROFILE                   "
                                                   >> 388            << "      CUMULATIVE  RADIAL PROFILE" << G4endl << G4endl;
                                                   >> 389 
                                                   >> 390     G4cout << "        bin   " << "           Mean         rms         "
                                                   >> 391            << "        bin "   << "           Mean      rms \n" << G4endl;
                                                   >> 392 
                                                   >> 393     for (i=0; i<nRbin; i++)
                                                   >> 394      {
                                                   >> 395        G4double inf=i*dRradl, sup=inf+dRradl;
                                                   >> 396 
                                                   >> 397        G4cout << std::setw(8) << inf << "->"
                                                   >> 398               << std::setw(5) << sup << " radl: "
                                                   >> 399               << std::setw(7) << MeanERadial[i] << "%  "
                                                   >> 400               << std::setw(9) << rmsERadial[i] << "%       "
                                                   >> 401               << "      0->" << std::setw(5) << sup << " radl: "
                                                   >> 402               << std::setw(7) << MeanERadialCumul[i] << "%  "
                                                   >> 403               << std::setw(7) << rmsERadialCumul[i] << "% "
                                                   >> 404               <<G4endl;
                                                   >> 405      }
                                                   >> 406   }
                                                   >> 407   
                                                   >> 408   G4cout << "\n SUMMARY \n" << G4endl;
                                                   >> 409   G4cout << " energy deposit : "
                                                   >> 410          << std::setw(7)  << MeanELongitCumul[nLbin-1] << " % E0 +- "
                                                   >> 411          << std::setw(7)  <<  rmsELongitCumul[nLbin-1] << " % E0" << G4endl;
                                                   >> 412   G4cout << " charged traklen: "
                                                   >> 413          << std::setw(7)  << MeanChargTrLength << " radl +- "
                                                   >> 414          << std::setw(7)  <<  rmsChargTrLength << " radl" << G4endl;
                                                   >> 415   G4cout << " neutral traklen: "
                                                   >> 416          << std::setw(7)  << MeanNeutrTrLength << " radl +- "
                                                   >> 417          << std::setw(7)  <<  rmsNeutrTrLength << " radl" << G4endl;
                                                   >> 418    
                                                   >> 419   if (iMoliere > 0. ) {
                                                   >> 420     G4double RMoliere1 = iMoliere*Det->GetdRradl();
                                                   >> 421     G4double RMoliere2 = iMoliere*Det->GetdRlength();    
                                                   >> 422     G4cout << "\n " << EMoliere << " % confinement: radius = "
                                                   >> 423            << RMoliere1 << " radl  ("
                                                   >> 424      << G4BestUnit( RMoliere2, "Length") << ")" << G4endl;
                                                   >> 425   }    
                                                   >> 426 
                                                   >> 427   G4cout.setf(mode,std::ios::floatfield);
                                                   >> 428   G4cout.precision(prec);
                                                   >> 429 
                                                   >> 430   // save histos and close AnalysisFactory
                                                   >> 431   cleanHisto();
                                                   >> 432   
154   // show Rndm status                             433   // show Rndm status
155   if (isMaster) G4Random::showEngineStatus();  << 434   CLHEP::HepRandom::showEngineStatus();
156                                                   435 
157   // save histos and close analysis            << 436   // Acceptance
158   fAnalysisManager->Write();                   << 437   if(limittrue < DBL_MAX) {
159   fAnalysisManager->CloseFile();               << 438     EmAcceptance acc;
160   fAnalysisManager->Clear();                   << 439     acc.BeginOfAcceptance("Total Energy in Absorber",NbOfEvents);
                                                   >> 440     G4double e = MeanELongitCumul[nLbin-1]/100.;
                                                   >> 441     G4double r = rmsELongitCumul[nLbin-1]/100.;
                                                   >> 442     acc.EmAcceptanceGauss("Edep",NbOfEvents,e,edeptrue,rmstrue,limittrue);
                                                   >> 443     acc.EmAcceptanceGauss("Erms",NbOfEvents,r,rmstrue,rmstrue,2.0*limittrue);
                                                   >> 444     acc.EndOfAcceptance();
                                                   >> 445   }
161 }                                                 446 }
162                                                   447 
163 //....oooOO0OOooo........oooOO0OOooo........oo    448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
164                                                   449 
165 void RunAction::SetEdepAndRMS(G4ThreeVector Va    450 void RunAction::SetEdepAndRMS(G4ThreeVector Value)
166 {                                                 451 {
167   fEdeptrue = Value(0);                        << 452   edeptrue = Value(0);
168   fRmstrue = Value(1);                         << 453   rmstrue  = Value(1);
169   fLimittrue = Value(2);                       << 454   limittrue= Value(2);
170 }                                              << 
171 //....oooOO0OOooo........oooOO0OOooo........oo << 
172                                                << 
173 void RunAction::SetVerbose(G4int val)          << 
174 {                                              << 
175   fVerbose = val;                              << 
176   if (fRun) fRun->SetVerbose(val);             << 
177 }                                                 455 }
178                                                   456 
179 //....oooOO0OOooo........oooOO0OOooo........oo    457 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180                                                   458