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.4)


  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.26 2010/11/09 21:25:15 asaim Exp $
 27 /// \brief Implementation of the RunAction cla <<  27 // GEANT4 tag $Name: geant4-09-04 $
 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    = "root";
 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 = "";
 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   histo[4] = hf->createHistogram1D( "5","rms on longit Edep (% of E inc)",
                                                   >> 147                                     nLbin,0.,nLbin*dLradl);
                                                   >> 148 
                                                   >> 149   G4double Zmin=0.5*dLradl, Zmax=Zmin+nLbin*dLradl;
                                                   >> 150   histo[5] = hf->createHistogram1D( "6","cumul longit energy dep (% of E inc)",
                                                   >> 151                                     nLbin,Zmin,Zmax);
                                                   >> 152             
                                                   >> 153   histo[6] = hf->createHistogram1D( "7","rms on cumul longit Edep (% of E inc)",
                                                   >> 154                                     nLbin,Zmin,Zmax);
                                                   >> 155 
                                                   >> 156   histo[7] = hf->createHistogram1D( "8","radial energy profile (% of E inc)",
                                                   >> 157                                     nRbin,0.,nRbin*dRradl);
                                                   >> 158                         
                                                   >> 159   histo[8] = hf->createHistogram1D( "9","rms on radial Edep (% of E inc)",
                                                   >> 160                                     nRbin,0.,nRbin*dRradl);     
                                                   >> 161 
                                                   >> 162   G4double Rmin=0.5*dRradl, Rmax=Rmin+nRbin*dRradl;
                                                   >> 163   histo[9] = hf->createHistogram1D("10","cumul radial energy dep (% of E inc)",
                                                   >> 164                                     nRbin,Rmin,Rmax);
                                                   >> 165 
                                                   >> 166   histo[10]= hf->createHistogram1D("11","rms on cumul radial Edep (% of E inc)",
                                                   >> 167                                     nRbin,Rmin,Rmax);       
                                                   >> 168             
                                                   >> 169  delete hf;
                                                   >> 170  G4cout << "\n----> Histogram Tree is opened in " << histoName[1] << G4endl;
                                                   >> 171 #endif
                                                   >> 172 }
 99                                                   173 
100   fAnalysisManager->CreateH1("h5", "rms on lon << 174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101                                                   175 
102   G4double Zmin = 0.5 * dLradl, Zmax = Zmin +  << 176 void RunAction::cleanHisto()
103   fAnalysisManager->CreateH1("h6", "cumul long << 177 {
                                                   >> 178 #ifdef G4ANALYSIS_USE
                                                   >> 179   if(tree) {
                                                   >> 180     tree->commit();       // Writing the histograms to the file
                                                   >> 181     tree->close();        // and closing the tree (and the file)
                                                   >> 182     G4cout << "\n----> Histogram Tree is saved in " << histoName[1] << G4endl; 
                                                   >> 183    
                                                   >> 184     delete tree;
                                                   >> 185     tree = 0;
                                                   >> 186   }
                                                   >> 187 #endif
                                                   >> 188 }
104                                                   189 
105   fAnalysisManager->CreateH1("h7", "rms on cum << 190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
106                                                   191 
107   fAnalysisManager->CreateH1("h8", "radial ene << 192 void RunAction::BeginOfRunAction(const G4Run* aRun)
                                                   >> 193 {
                                                   >> 194   G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
108                                                   195 
109   fAnalysisManager->CreateP1("p8", "radial ene << 196   // save Rndm status
110                              0., 1000.);       << 197   G4RunManager::GetRunManager()->SetRandomNumberStore(true);
                                                   >> 198   CLHEP::HepRandom::showEngineStatus();
111                                                   199 
112   fAnalysisManager->CreateH1("h9", "rms on rad << 200   nLbin = Det->GetnLtot();
                                                   >> 201   nRbin = Det->GetnRtot();
113                                                   202 
114   G4double Rmin = 0.5 * dRradl, Rmax = Rmin +  << 203   //initialize arrays of cumulative energy deposition
115   fAnalysisManager->CreateH1("h10", "cumul rad << 204   //
                                                   >> 205   for (G4int i=0; i<nLbin; i++) 
                                                   >> 206      sumELongit[i]=sumE2Longit[i]=sumELongitCumul[i]=sumE2LongitCumul[i]=0.;
                                                   >> 207   
                                                   >> 208   for (G4int j=0; j<nRbin; j++)
                                                   >> 209      sumERadial[j]=sumE2Radial[j]=sumERadialCumul[j]=sumE2RadialCumul[j]=0.;
116                                                   210 
117   fAnalysisManager->CreateH1("h11", "rms on cu << 211   //initialize track length
                                                   >> 212   sumChargTrLength=sum2ChargTrLength=sumNeutrTrLength=sum2NeutrTrLength=0.;
118                                                   213 
119   G4String fileName = fAnalysisManager->GetFil << 214   //histograms
120   G4String extension = fAnalysisManager->GetFi << 215   //
121   G4String fullFileName = fileName + "." + ext << 216   if (aRun->GetRunID() == 0) bookHisto();
122   G4cout << "\n----> Histogram file is opened  << 
123 }                                                 217 }
124                                                   218 
125 //....oooOO0OOooo........oooOO0OOooo........oo    219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
126                                                   220 
127 G4Run* RunAction::GenerateRun()                << 221 void RunAction::fillPerEvent()
128 {                                                 222 {
129   fRun = new Run(fDet, fKin);                  << 223   //accumulate statistic
130   fRun->SetVerbose(fVerbose);                  << 224   //
131   return fRun;                                 << 225   G4double dLCumul = 0.;
                                                   >> 226   for (G4int i=0; i<nLbin; i++)
                                                   >> 227      {
                                                   >> 228       sumELongit[i]  += dEdL[i];
                                                   >> 229       sumE2Longit[i] += dEdL[i]*dEdL[i];
                                                   >> 230       dLCumul        += dEdL[i];
                                                   >> 231       sumELongitCumul[i]  += dLCumul;
                                                   >> 232       sumE2LongitCumul[i] += dLCumul*dLCumul;
                                                   >> 233      }
                                                   >> 234 
                                                   >> 235   G4double dRCumul = 0.;
                                                   >> 236   for (G4int j=0; j<nRbin; j++)
                                                   >> 237      {
                                                   >> 238       sumERadial[j]  += dEdR[j];
                                                   >> 239       sumE2Radial[j] += dEdR[j]*dEdR[j];
                                                   >> 240       dRCumul        += dEdR[j];
                                                   >> 241       sumERadialCumul[j]  += dRCumul;
                                                   >> 242       sumE2RadialCumul[j] += dRCumul*dRCumul;
                                                   >> 243      }
                                                   >> 244 
                                                   >> 245   sumChargTrLength  += ChargTrLength;
                                                   >> 246   sum2ChargTrLength += ChargTrLength*ChargTrLength;
                                                   >> 247   sumNeutrTrLength  += NeutrTrLength;
                                                   >> 248   sum2NeutrTrLength += NeutrTrLength*NeutrTrLength;
                                                   >> 249 
                                                   >> 250 #ifdef G4ANALYSIS_USE
                                                   >> 251   //fill histograms
                                                   >> 252   //
                                                   >> 253   if(tree) {
                                                   >> 254     G4double Ekin=Kin->GetParticleGun()->GetParticleEnergy();
                                                   >> 255     G4double mass=Kin->GetParticleGun()->GetParticleDefinition()->GetPDGMass();
                                                   >> 256     G4double radl=Det->GetMaterial()->GetRadlen();
                                                   >> 257 
                                                   >> 258     histo[0]->fill(100.*dLCumul/(Ekin+mass));
                                                   >> 259     histo[1]->fill(ChargTrLength/radl);
                                                   >> 260     histo[2]->fill(NeutrTrLength/radl);
                                                   >> 261   }
                                                   >> 262 #endif
132 }                                                 263 }
133                                                   264 
134 //....oooOO0OOooo........oooOO0OOooo........oo    265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135                                                   266 
136 void RunAction::BeginOfRunAction(const G4Run*) << 267 void RunAction::EndOfRunAction(const G4Run* aRun)
137 {                                                 268 {
138   // show Rndm status                          << 269   //compute and print statistic
139   if (isMaster) G4Random::showEngineStatus();  << 270   //
                                                   >> 271   G4int NbOfEvents = aRun->GetNumberOfEvent();
                                                   >> 272   G4double kinEnergy = Kin->GetParticleGun()->GetParticleEnergy();
                                                   >> 273   assert(NbOfEvents*kinEnergy > 0);
                                                   >> 274 
                                                   >> 275   G4double mass=Kin->GetParticleGun()->GetParticleDefinition()->GetPDGMass();
                                                   >> 276   G4double norme = 100./(NbOfEvents*(kinEnergy+mass));
140                                                   277 
141   // histograms                                << 278   //longitudinal
142   //                                              279   //
143   BookHisto();                                 << 280   G4double dLradl = Det->GetdLradl();
144 }                                              << 
145                                                   281 
146 //....oooOO0OOooo........oooOO0OOooo........oo << 282   MyVector MeanELongit(nLbin),      rmsELongit(nLbin);
                                                   >> 283   MyVector MeanELongitCumul(nLbin), rmsELongitCumul(nLbin);
147                                                   284 
148 void RunAction::EndOfRunAction(const G4Run*)   << 285   G4int i;
149 {                                              << 286   for (i=0; i<nLbin; i++)
150   // compute and print statistic               << 287    {
                                                   >> 288     MeanELongit[i] = norme*sumELongit[i];
                                                   >> 289      rmsELongit[i] = norme*std::sqrt(std::fabs(NbOfEvents*sumE2Longit[i]
                                                   >> 290                                 - sumELongit[i]*sumELongit[i]));
                                                   >> 291 
                                                   >> 292     MeanELongitCumul[i] = norme*sumELongitCumul[i];
                                                   >> 293      rmsELongitCumul[i] = norme*std::sqrt(std::fabs(NbOfEvents*sumE2LongitCumul[i]
                                                   >> 294                                     - sumELongitCumul[i]*sumELongitCumul[i]));
                                                   >> 295 
                                                   >> 296 
                                                   >> 297 #ifdef G4ANALYSIS_USE
                                                   >> 298     if(tree) {
                                                   >> 299       G4double bin = (i+0.5)*dLradl;
                                                   >> 300       histo[3]->fill(bin,MeanELongit[i]/dLradl);
                                                   >> 301       histo[4]->fill(bin, rmsELongit[i]/dLradl);      
                                                   >> 302       bin = (i+1)*dLradl;
                                                   >> 303       histo[5]->fill(bin,MeanELongitCumul[i]);
                                                   >> 304       histo[6]->fill(bin, rmsELongitCumul[i]);
                                                   >> 305     }
                                                   >> 306 #endif
                                                   >> 307    }
                                                   >> 308 
                                                   >> 309   //radial
151   //                                              310   //
152   if (isMaster) fRun->EndOfRun(fEdeptrue, fRms << 311   G4double dRradl = Det->GetdRradl();
                                                   >> 312 
                                                   >> 313   MyVector MeanERadial(nRbin),      rmsERadial(nRbin);
                                                   >> 314   MyVector MeanERadialCumul(nRbin), rmsERadialCumul(nRbin);
153                                                   315 
                                                   >> 316   for (i=0; i<nRbin; i++)
                                                   >> 317    {
                                                   >> 318     MeanERadial[i] = norme*sumERadial[i];
                                                   >> 319      rmsERadial[i] = norme*std::sqrt(std::fabs(NbOfEvents*sumE2Radial[i]
                                                   >> 320                                 - sumERadial[i]*sumERadial[i]));
                                                   >> 321 
                                                   >> 322     MeanERadialCumul[i] = norme*sumERadialCumul[i];
                                                   >> 323      rmsERadialCumul[i] = norme*std::sqrt(std::fabs(NbOfEvents*sumE2RadialCumul[i]
                                                   >> 324                                     - sumERadialCumul[i]*sumERadialCumul[i]));
                                                   >> 325 
                                                   >> 326 #ifdef G4ANALYSIS_USE
                                                   >> 327     if(tree) {
                                                   >> 328       G4double bin = (i+0.5)*dRradl;
                                                   >> 329       histo[7]->fill(bin,MeanERadial[i]/dRradl);
                                                   >> 330       histo[8]->fill(bin, rmsERadial[i]/dRradl);      
                                                   >> 331       bin = (i+1)*dRradl;
                                                   >> 332       histo[9] ->fill(bin,MeanERadialCumul[i]);
                                                   >> 333       histo[10]->fill(bin, rmsERadialCumul[i]);
                                                   >> 334     }
                                                   >> 335 #endif
                                                   >> 336    }
                                                   >> 337 
                                                   >> 338   //find Moliere confinement
                                                   >> 339   //
                                                   >> 340   const G4double EMoliere = 90.;
                                                   >> 341   G4double iMoliere = 0.;
                                                   >> 342   if ((MeanERadialCumul[0]       <= EMoliere) &&
                                                   >> 343       (MeanERadialCumul[nRbin-1] >= EMoliere)) {
                                                   >> 344     G4int imin = 0;
                                                   >> 345     while( (imin < nRbin-1) && (MeanERadialCumul[imin] < EMoliere) ) imin++;
                                                   >> 346     G4double ratio = (EMoliere - MeanERadialCumul[imin]) /
                                                   >> 347                      (MeanERadialCumul[imin+1] - MeanERadialCumul[imin]);
                                                   >> 348     iMoliere = 1. + imin + ratio;
                                                   >> 349   }          
                                                   >> 350       
                                                   >> 351   //track length
                                                   >> 352   //
                                                   >> 353   norme = 1./(NbOfEvents*(Det->GetMaterial()->GetRadlen()));
                                                   >> 354   G4double MeanChargTrLength = norme*sumChargTrLength;
                                                   >> 355   G4double  rmsChargTrLength = 
                                                   >> 356             norme*std::sqrt(std::fabs(NbOfEvents*sum2ChargTrLength
                                                   >> 357                                          - sumChargTrLength*sumChargTrLength));
                                                   >> 358 
                                                   >> 359   G4double MeanNeutrTrLength = norme*sumNeutrTrLength;
                                                   >> 360   G4double  rmsNeutrTrLength = 
                                                   >> 361             norme*std::sqrt(std::fabs(NbOfEvents*sum2NeutrTrLength
                                                   >> 362                                          - sumNeutrTrLength*sumNeutrTrLength));
                                                   >> 363 
                                                   >> 364   //print
                                                   >> 365   //
                                                   >> 366 
                                                   >> 367   std::ios::fmtflags mode = G4cout.flags();
                                                   >> 368   G4cout.setf(std::ios::fixed,std::ios::floatfield);
                                                   >> 369   G4int  prec = G4cout.precision(2);
                                                   >> 370   
                                                   >> 371   if (verbose) {
                                                   >> 372 
                                                   >> 373     G4cout << "                 LATERAL PROFILE                   "
                                                   >> 374            << "      CUMULATIVE LATERAL PROFILE" << G4endl << G4endl;
                                                   >> 375 
                                                   >> 376     G4cout << "        bin   " << "           Mean         rms         "
                                                   >> 377            << "        bin "   << "           Mean      rms \n" << G4endl;
                                                   >> 378 
                                                   >> 379     for (i=0; i<nLbin; i++)
                                                   >> 380      {
                                                   >> 381        G4double inf=i*dLradl, sup=inf+dLradl;
                                                   >> 382 
                                                   >> 383        G4cout << std::setw(8) << inf << "->"
                                                   >> 384               << std::setw(5) << sup << " radl: "
                                                   >> 385               << std::setw(7) << MeanELongit[i] << "%  "
                                                   >> 386               << std::setw(9) << rmsELongit[i] << "%       "
                                                   >> 387               << "      0->" << std::setw(5) << sup << " radl: "
                                                   >> 388               << std::setw(7) << MeanELongitCumul[i] << "%  "
                                                   >> 389               << std::setw(7) << rmsELongitCumul[i] << "% "
                                                   >> 390               <<G4endl;
                                                   >> 391      }
                                                   >> 392 
                                                   >> 393     G4cout << G4endl << G4endl << G4endl;
                                                   >> 394 
                                                   >> 395     G4cout << "                  RADIAL PROFILE                   "
                                                   >> 396            << "      CUMULATIVE  RADIAL PROFILE" << G4endl << G4endl;
                                                   >> 397 
                                                   >> 398     G4cout << "        bin   " << "           Mean         rms         "
                                                   >> 399            << "        bin "   << "           Mean      rms \n" << G4endl;
                                                   >> 400 
                                                   >> 401     for (i=0; i<nRbin; i++)
                                                   >> 402      {
                                                   >> 403        G4double inf=i*dRradl, sup=inf+dRradl;
                                                   >> 404 
                                                   >> 405        G4cout << std::setw(8) << inf << "->"
                                                   >> 406               << std::setw(5) << sup << " radl: "
                                                   >> 407               << std::setw(7) << MeanERadial[i] << "%  "
                                                   >> 408               << std::setw(9) << rmsERadial[i] << "%       "
                                                   >> 409               << "      0->" << std::setw(5) << sup << " radl: "
                                                   >> 410               << std::setw(7) << MeanERadialCumul[i] << "%  "
                                                   >> 411               << std::setw(7) << rmsERadialCumul[i] << "% "
                                                   >> 412               <<G4endl;
                                                   >> 413      }
                                                   >> 414   }
                                                   >> 415   
                                                   >> 416   G4cout << "\n SUMMARY \n" << G4endl;
                                                   >> 417   G4cout << " energy deposit : "
                                                   >> 418          << std::setw(7)  << MeanELongitCumul[nLbin-1] << " % E0 +- "
                                                   >> 419          << std::setw(7)  <<  rmsELongitCumul[nLbin-1] << " % E0" << G4endl;
                                                   >> 420   G4cout << " charged traklen: "
                                                   >> 421          << std::setw(7)  << MeanChargTrLength << " radl +- "
                                                   >> 422          << std::setw(7)  <<  rmsChargTrLength << " radl" << G4endl;
                                                   >> 423   G4cout << " neutral traklen: "
                                                   >> 424          << std::setw(7)  << MeanNeutrTrLength << " radl +- "
                                                   >> 425          << std::setw(7)  <<  rmsNeutrTrLength << " radl" << G4endl;
                                                   >> 426    
                                                   >> 427   if (iMoliere > 0. ) {
                                                   >> 428     G4double RMoliere1 = iMoliere*Det->GetdRradl();
                                                   >> 429     G4double RMoliere2 = iMoliere*Det->GetdRlength();    
                                                   >> 430     G4cout << "\n " << EMoliere << " % confinement: radius = "
                                                   >> 431            << RMoliere1 << " radl  ("
                                                   >> 432      << G4BestUnit( RMoliere2, "Length") << ")" << G4endl;
                                                   >> 433   }    
                                                   >> 434 
                                                   >> 435   G4cout.setf(mode,std::ios::floatfield);
                                                   >> 436   G4cout.precision(prec);
                                                   >> 437 
                                                   >> 438   // save histos and close AnalysisFactory
                                                   >> 439   cleanHisto();
                                                   >> 440   
154   // show Rndm status                             441   // show Rndm status
155   if (isMaster) G4Random::showEngineStatus();  << 442   CLHEP::HepRandom::showEngineStatus();
156                                                   443 
157   // save histos and close analysis            << 444   // Acceptance
158   fAnalysisManager->Write();                   << 445   if(limittrue < DBL_MAX) {
159   fAnalysisManager->CloseFile();               << 446     EmAcceptance acc;
160   fAnalysisManager->Clear();                   << 447     acc.BeginOfAcceptance("Total Energy in Absorber",NbOfEvents);
                                                   >> 448     G4double e = MeanELongitCumul[nLbin-1]/100.;
                                                   >> 449     G4double r = rmsELongitCumul[nLbin-1]/100.;
                                                   >> 450     acc.EmAcceptanceGauss("Edep",NbOfEvents,e,edeptrue,rmstrue,limittrue);
                                                   >> 451     acc.EmAcceptanceGauss("Erms",NbOfEvents,r,rmstrue,rmstrue,2.0*limittrue);
                                                   >> 452     acc.EndOfAcceptance();
                                                   >> 453   }
161 }                                                 454 }
162                                                   455 
163 //....oooOO0OOooo........oooOO0OOooo........oo    456 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
164                                                   457 
165 void RunAction::SetEdepAndRMS(G4ThreeVector Va    458 void RunAction::SetEdepAndRMS(G4ThreeVector Value)
166 {                                                 459 {
167   fEdeptrue = Value(0);                        << 460   edeptrue = Value(0);
168   fRmstrue = Value(1);                         << 461   rmstrue  = Value(1);
169   fLimittrue = Value(2);                       << 462   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 }                                                 463 }
178                                                   464 
179 //....oooOO0OOooo........oooOO0OOooo........oo    465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180                                                   466