Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm3/src/Run.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/TestEm3/src/Run.cc (Version 11.3.0) and /examples/extended/electromagnetic/TestEm3/src/Run.cc (Version 11.1.1)


  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/TestEm3/src/Run.cc       26 /// \file electromagnetic/TestEm3/src/Run.cc
 27 /// \brief Implementation of the Run class         27 /// \brief Implementation of the Run class
 28 //                                                 28 //
 29 //                                             <<  29 // 
 30 //....oooOO0OOooo........oooOO0OOooo........oo     30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oo     31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32                                                    32 
 33 #include "Run.hh"                                  33 #include "Run.hh"
 34                                                << 
 35 #include "DetectorConstruction.hh"                 34 #include "DetectorConstruction.hh"
 36 #include "EmAcceptance.hh"                     << 
 37 #include "HistoManager.hh"                     << 
 38 #include "PrimaryGeneratorAction.hh"               35 #include "PrimaryGeneratorAction.hh"
                                                   >>  36 #include "HistoManager.hh"
                                                   >>  37 #include "EmAcceptance.hh"
 39                                                    38 
 40 #include "G4Electron.hh"                       << 
 41 #include "G4Gamma.hh"                          << 
 42 #include "G4ParticleDefinition.hh"             << 
 43 #include "G4ParticleTable.hh"                      39 #include "G4ParticleTable.hh"
 44 #include "G4Positron.hh"                       <<  40 #include "G4ParticleDefinition.hh"
 45 #include "G4SystemOfUnits.hh"                  << 
 46 #include "G4Track.hh"                              41 #include "G4Track.hh"
                                                   >>  42 #include "G4Gamma.hh"
                                                   >>  43 #include "G4Electron.hh"
                                                   >>  44 #include "G4Positron.hh"
                                                   >>  45 
 47 #include "G4UnitsTable.hh"                         46 #include "G4UnitsTable.hh"
                                                   >>  47 #include "G4SystemOfUnits.hh"
 48                                                    48 
 49 #include <iomanip>                                 49 #include <iomanip>
 50                                                    50 
 51 //....oooOO0OOooo........oooOO0OOooo........oo     51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 52                                                    52 
 53 Run::Run(DetectorConstruction* det) : fDetecto <<  53 Run::Run(DetectorConstruction* det)
                                                   >>  54 : G4Run(),
                                                   >>  55   fDetector(det), 
                                                   >>  56   fParticle(nullptr), fEkin(0.),
                                                   >>  57   fChargedStep(0), fNeutralStep(0),
                                                   >>  58   fN_gamma(0), fN_elec(0), fN_pos(0),
                                                   >>  59   fApplyLimit(false)
 54 {                                                  60 {
 55   // initialize cumulative quantities          <<  61   //initialize cumulative quantities
 56   //                                               62   //
 57   for (G4int k = 0; k < kMaxAbsor; k++) {      <<  63   for (G4int k=0; k<kMaxAbsor; k++) {
 58     fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = <<  64     fSumEAbs[k] = fSum2EAbs[k]  = fSumLAbs[k] = fSum2LAbs[k] = 0.;
 59     fEnergyDeposit[k].clear();                     65     fEnergyDeposit[k].clear();
 60     fEdeptrue[k] = fRmstrue[k] = 1.;               66     fEdeptrue[k] = fRmstrue[k] = 1.;
 61     fLimittrue[k] = DBL_MAX;                   <<  67     fLimittrue[k] = DBL_MAX;  
 62   }                                                68   }
 63                                                <<  69   
 64   fEdepTot = fEdepTot2 = 0.;                   <<  70   //initialize Eflow
 65   fEleakTot = fEleakTot2 = 0.;                 << 
 66   fEtotal = fEtotal2 = 0.;                     << 
 67                                                << 
 68   // initialize Eflow                          << 
 69   //                                               71   //
 70   G4int nbPlanes = (fDetector->GetNbOfLayers() <<  72   G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2;
 71   fEnergyFlow.resize(nbPlanes);                    73   fEnergyFlow.resize(nbPlanes);
 72   fLateralEleak.resize(nbPlanes);                  74   fLateralEleak.resize(nbPlanes);
 73   for (G4int k = 0; k < nbPlanes; k++) {       <<  75   for (G4int k=0; k<nbPlanes; k++) {fEnergyFlow[k] = fLateralEleak[k] = 0.; }  
 74     fEnergyFlow[k] = fLateralEleak[k] = 0.;    << 
 75   }                                            << 
 76 }                                                  76 }
 77                                                    77 
 78 //....oooOO0OOooo........oooOO0OOooo........oo     78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 79                                                    79 
                                                   >>  80 Run::~Run()
                                                   >>  81 { }
                                                   >>  82 
                                                   >>  83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >>  84 
 80 void Run::SetPrimary(G4ParticleDefinition* par     85 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
 81 {                                              <<  86 { 
 82   fParticle = particle;                            87   fParticle = particle;
 83   fEkin = energy;                                  88   fEkin = energy;
 84 }                                                  89 }
 85                                                    90 
 86 //....oooOO0OOooo........oooOO0OOooo........oo     91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 87                                                    92 
 88 void Run::FillPerEvent(G4int kAbs, G4double EA     93 void Run::FillPerEvent(G4int kAbs, G4double EAbs, G4double LAbs)
 89 {                                                  94 {
 90   // accumulate statistic with restriction     <<  95   //accumulate statistic with restriction
 91   //                                               96   //
 92   if (fApplyLimit) fEnergyDeposit[kAbs].push_b <<  97   if(fApplyLimit) fEnergyDeposit[kAbs].push_back(EAbs);
 93   fSumEAbs[kAbs] += EAbs;                      <<  98   fSumEAbs[kAbs]  += EAbs;  fSum2EAbs[kAbs]  += EAbs*EAbs;
 94   fSum2EAbs[kAbs] += EAbs * EAbs;              <<  99   fSumLAbs[kAbs]  += LAbs;  fSum2LAbs[kAbs]  += LAbs*LAbs;
 95   fSumLAbs[kAbs] += LAbs;                      << 
 96   fSum2LAbs[kAbs] += LAbs * LAbs;              << 
 97 }                                              << 
 98                                                << 
 99 //....oooOO0OOooo........oooOO0OOooo........oo << 
100                                                << 
101 void Run::SumEnergies(G4double edeptot, G4doub << 
102 {                                              << 
103   fEdepTot += edeptot;                         << 
104   fEdepTot2 += edeptot * edeptot;              << 
105   fEleakTot += eleak;                          << 
106   fEleakTot2 += eleak * eleak;                 << 
107                                                << 
108   G4double etotal = edeptot + eleak;           << 
109   fEtotal += etotal;                           << 
110   fEtotal2 += etotal * etotal;                 << 
111 }                                                 100 }
112                                                   101 
113 //....oooOO0OOooo........oooOO0OOooo........oo    102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114                                                   103 
115 void Run::SumEnergyFlow(G4int plane, G4double     104 void Run::SumEnergyFlow(G4int plane, G4double Eflow)
116 {                                                 105 {
117   fEnergyFlow[plane] += Eflow;                    106   fEnergyFlow[plane] += Eflow;
118 }                                                 107 }
119                                                   108 
120 //....oooOO0OOooo........oooOO0OOooo........oo    109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
121                                                << 110                                             
122 void Run::SumLateralEleak(G4int cell, G4double    111 void Run::SumLateralEleak(G4int cell, G4double Eflow)
123 {                                                 112 {
124   fLateralEleak[cell] += Eflow;                   113   fLateralEleak[cell] += Eflow;
125 }                                                 114 }
126                                                   115 
127 //....oooOO0OOooo........oooOO0OOooo........oo    116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128                                                << 117                                             
129 void Run::AddChargedStep()                     << 118 void Run::AddChargedStep() 
130 {                                                 119 {
131   fChargedStep += 1.0;                         << 120   fChargedStep += 1.0; 
132 }                                                 121 }
133                                                   122 
134 //....oooOO0OOooo........oooOO0OOooo........oo    123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135                                                << 124     
136 void Run::AddNeutralStep()                     << 125 void Run::AddNeutralStep() 
137 {                                                 126 {
138   fNeutralStep += 1.0;                         << 127   fNeutralStep += 1.0; 
139 }                                                 128 }
140                                                << 129     
141 //....oooOO0OOooo........oooOO0OOooo........oo    130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
142                                                   131 
143 void Run::AddSecondaryTrack(const G4Track* tra    132 void Run::AddSecondaryTrack(const G4Track* track)
144 {                                                 133 {
145   const G4ParticleDefinition* d = track->GetDe    134   const G4ParticleDefinition* d = track->GetDefinition();
146   if (d == G4Gamma::Gamma()) {                 << 135   if(d == G4Gamma::Gamma()) { ++fN_gamma; }
147     ++fN_gamma;                                << 136   else if (d == G4Electron::Electron()) { ++fN_elec; }
148   }                                            << 137   else if (d == G4Positron::Positron()) { ++fN_pos; }
149   else if (d == G4Electron::Electron()) {      << 
150     ++fN_elec;                                 << 
151   }                                            << 
152   else if (d == G4Positron::Positron()) {      << 
153     ++fN_pos;                                  << 
154   }                                            << 
155 }                                                 138 }
156                                                << 139  
157 //....oooOO0OOooo........oooOO0OOooo........oo    140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
158                                                   141 
159 void Run::Merge(const G4Run* run)                 142 void Run::Merge(const G4Run* run)
160 {                                                 143 {
161   const Run* localRun = static_cast<const Run*    144   const Run* localRun = static_cast<const Run*>(run);
162                                                   145 
163   // pass information about primary particle      146   // pass information about primary particle
164   fParticle = localRun->fParticle;                147   fParticle = localRun->fParticle;
165   fEkin = localRun->fEkin;                     << 148   fEkin     = localRun->fEkin;
166                                                   149 
167   // accumulate sums                              150   // accumulate sums
168   //                                              151   //
169   for (G4int k = 0; k < kMaxAbsor; k++) {      << 152   for (G4int k=0; k<kMaxAbsor; k++) {
170     fSumEAbs[k] += localRun->fSumEAbs[k];      << 153     fSumEAbs[k]  += localRun->fSumEAbs[k]; 
171     fSum2EAbs[k] += localRun->fSum2EAbs[k];    << 154     fSum2EAbs[k] += localRun->fSum2EAbs[k]; 
172     fSumLAbs[k] += localRun->fSumLAbs[k];      << 155     fSumLAbs[k]  += localRun->fSumLAbs[k]; 
173     fSum2LAbs[k] += localRun->fSum2LAbs[k];       156     fSum2LAbs[k] += localRun->fSum2LAbs[k];
174   }                                               157   }
175                                                << 158    
176   fEdepTot += localRun->fEdepTot;              << 159   G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2;
177   fEdepTot2 += localRun->fEdepTot2;            << 160   for (G4int k=0; k<nbPlanes; k++) {
178                                                << 161     fEnergyFlow[k]   += localRun->fEnergyFlow[k];
179   fEleakTot += localRun->fEleakTot;            << 
180   fEleakTot2 += localRun->fEleakTot2;          << 
181                                                << 
182   fEtotal += localRun->fEtotal;                << 
183   fEtotal2 += localRun->fEtotal2;              << 
184                                                << 
185   G4int nbPlanes = (fDetector->GetNbOfLayers() << 
186   for (G4int k = 0; k < nbPlanes; k++) {       << 
187     fEnergyFlow[k] += localRun->fEnergyFlow[k] << 
188     fLateralEleak[k] += localRun->fLateralElea    162     fLateralEleak[k] += localRun->fLateralEleak[k];
189   }                                               163   }
190                                                << 164   
191   for (G4int k=0; k<kMaxAbsor; k++) {          << 165        
192     fEnergyDeposit[k].insert(fEnergyDeposit[k] << 166   fChargedStep += localRun->fChargedStep;  
193     localRun->fEnergyDeposit[k].begin(), local << 
194   }                                            << 
195                                                << 
196   fChargedStep += localRun->fChargedStep;      << 
197   fNeutralStep += localRun->fNeutralStep;         167   fNeutralStep += localRun->fNeutralStep;
198                                                << 168   
199   fN_gamma += localRun->fN_gamma;              << 169   fN_gamma += localRun->fN_gamma;    
200   fN_elec += localRun->fN_elec;                << 170   fN_elec  += localRun->fN_elec;
201   fN_pos += localRun->fN_pos;                  << 171   fN_pos   += localRun->fN_pos;
202                                                << 172      
203   fApplyLimit = localRun->fApplyLimit;            173   fApplyLimit = localRun->fApplyLimit;
204                                                << 174   
205   for (G4int k = 0; k < kMaxAbsor; k++) {      << 175   for (G4int k=0; k<kMaxAbsor; k++) {
206     fEdeptrue[k] = localRun->fEdeptrue[k];     << 176     fEdeptrue[k]  = localRun->fEdeptrue[k]; 
207     fRmstrue[k] = localRun->fRmstrue[k];       << 177     fRmstrue[k]   = localRun->fRmstrue[k]; 
208     fLimittrue[k] = localRun->fLimittrue[k];   << 178     fLimittrue[k] = localRun->fLimittrue[k]; 
209   }                                            << 179   }
210                                                << 180     
211   G4Run::Merge(run);                           << 181   G4Run::Merge(run); 
212 }                                              << 182 } 
213                                                   183 
214 //....oooOO0OOooo........oooOO0OOooo........oo    184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
215                                                   185 
216 void Run::EndOfRun()                              186 void Run::EndOfRun()
217 {                                                 187 {
218   G4int nEvt = numberOfEvent;                     188   G4int nEvt = numberOfEvent;
219   G4double norm = G4double(nEvt);              << 189   G4double  norm = G4double(nEvt);
220   if (norm > 0) norm = 1. / norm;              << 190   if(norm > 0) norm = 1./norm;
221   G4double qnorm = std::sqrt(norm);               191   G4double qnorm = std::sqrt(norm);
222                                                   192 
223   fChargedStep *= norm;                           193   fChargedStep *= norm;
224   fNeutralStep *= norm;                           194   fNeutralStep *= norm;
225                                                   195 
226   // compute and print statistic               << 196   //compute and print statistic
227   //                                              197   //
228   G4double beamEnergy = fEkin;                    198   G4double beamEnergy = fEkin;
229   G4double sqbeam = std::sqrt(beamEnergy / GeV << 199   G4double sqbeam = std::sqrt(beamEnergy/GeV);
230                                                   200 
231   G4double MeanEAbs, MeanEAbs2, rmsEAbs, resol << 201   G4double MeanEAbs,MeanEAbs2,rmsEAbs,resolution,rmsres;
232   G4double MeanLAbs, MeanLAbs2, rmsLAbs;       << 202   G4double MeanLAbs,MeanLAbs2,rmsLAbs;
233                                                   203 
234   std::ios::fmtflags mode = G4cout.flags();       204   std::ios::fmtflags mode = G4cout.flags();
235   G4int prec = G4cout.precision(2);            << 205   G4int  prec = G4cout.precision(2);
236   G4cout << "\n-------------------------------    206   G4cout << "\n------------------------------------------------------------\n";
237   G4cout << std::setw(14) << "material" << std << 207   G4cout << std::setw(14) << "material"
238          << "sqrt(E0(GeV))*rmsE/Emean" << std: << 208          << std::setw(17) << "Edep       RMS"
239                                                << 209          << std::setw(33) << "sqrt(E0(GeV))*rmsE/Emean"
240   for (G4int k = 1; k <= fDetector->GetNbOfAbs << 210          << std::setw(23) << "total tracklen \n \n";
241     MeanEAbs = fSumEAbs[k] * norm;             << 211 
242     MeanEAbs2 = fSum2EAbs[k] * norm;           << 212   for (G4int k=1; k<=fDetector->GetNbOfAbsor(); k++)
243     rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - M << 213     {
244     // G4cout << "k= " << k << "  RMS= " <<  r << 214       MeanEAbs  = fSumEAbs[k]*norm;
245     //      << "  fApplyLimit: " << fApplyLimi << 215       MeanEAbs2 = fSum2EAbs[k]*norm;
246     if (fApplyLimit) {                         << 216       rmsEAbs  = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
247       G4int nn = 0;                            << 217       //G4cout << "k= " << k << "  RMS= " <<  rmsEAbs 
248       G4double sume = 0.0;                     << 218       //     << "  fApplyLimit: " << fApplyLimit << G4endl;
249       G4double sume2 = 0.0;                    << 219       if(fApplyLimit) {
250       // compute trancated means               << 220         G4int    nn    = 0;
251       G4double lim = rmsEAbs * 2.5;            << 221         G4double sume  = 0.0;
252       for (G4int i = 0; i < nEvt; i++) {       << 222         G4double sume2 = 0.0;
253         G4double e = (fEnergyDeposit[k])[i];   << 223         // compute trancated means  
254         if (std::abs(e - MeanEAbs) < lim) {    << 224         G4double lim   = rmsEAbs * 2.5;
255           sume += e;                           << 225         for(G4int i=0; i<nEvt; i++) {
256           sume2 += e * e;                      << 226           G4double e = (fEnergyDeposit[k])[i];
257           nn++;                                << 227           if(std::abs(e - MeanEAbs) < lim) {
                                                   >> 228             sume  += e;
                                                   >> 229             sume2 += e*e;
                                                   >> 230             nn++;
                                                   >> 231           }
258         }                                         232         }
                                                   >> 233         G4double norm1 = G4double(nn);
                                                   >> 234         if(norm1 > 0.0) norm1 = 1.0/norm1;
                                                   >> 235         MeanEAbs  = sume*norm1;
                                                   >> 236         MeanEAbs2 = sume2*norm1;
                                                   >> 237         rmsEAbs  = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
259       }                                           238       }
260       G4double norm1 = G4double(nn);           << 
261       if (norm1 > 0.0) norm1 = 1.0 / norm1;    << 
262       MeanEAbs = sume * norm1;                 << 
263       MeanEAbs2 = sume2 * norm1;               << 
264       rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - << 
265     }                                          << 
266                                                << 
267     resolution = (MeanEAbs > 0.) ? 100. * sqbe << 
268     rmsres = resolution * qnorm;               << 
269                                                << 
270     // Save mean and RMS                       << 
271     fSumEAbs[k] = MeanEAbs;                    << 
272     fSum2EAbs[k] = rmsEAbs;                    << 
273                                                << 
274     MeanLAbs = fSumLAbs[k] * norm;             << 
275     MeanLAbs2 = fSum2LAbs[k] * norm;           << 
276     rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - M << 
277                                                << 
278     // print                                   << 
279     //                                         << 
280     G4cout << std::setw(14) << fDetector->GetA << 
281            << std::setprecision(5) << std::set << 
282            << std::setprecision(4) << std::set << 
283            << resolution << " +- " << std::set << 
284            << std::setw(10) << G4BestUnit(Mean << 
285            << G4BestUnit(rmsLAbs, "Length") << << 
286   }                                            << 
287                                                   239 
288   // total energy deposited                    << 240       resolution= 100.*sqbeam*rmsEAbs/MeanEAbs;
289   //                                           << 241       rmsres    = resolution*qnorm;
290   fEdepTot *= norm;                            << 
291   fEdepTot2 *= norm;                           << 
292   G4double rmsEdep = std::sqrt(std::abs(fEdepT << 
293                                                << 
294   G4cout << "\n Total energy deposited = " <<  << 
295          << " +- " << G4BestUnit(rmsEdep, "Ene << 
296                                                << 
297   // Energy leakage                            << 
298   //                                           << 
299   fEleakTot *= norm;                           << 
300   fEleakTot2 *= norm;                          << 
301   G4double rmsEleak = std::sqrt(std::abs(fElea << 
302                                                << 
303   G4cout << " Energy leakage = " << G4BestUnit << 
304          << G4BestUnit(rmsEleak, "Energy") <<  << 
305                                                   242 
306   // total energy                              << 243       // Save mean and RMS
307   //                                           << 244       fSumEAbs[k] = MeanEAbs;
308   fEtotal *= norm;                             << 245       fSum2EAbs[k] = rmsEAbs;
309   fEtotal2 *= norm;                            << 246 
310   G4double rmsEtotal = std::sqrt(std::abs(fEto << 247       MeanLAbs  = fSumLAbs[k]*norm;
311                                                << 248       MeanLAbs2 = fSum2LAbs[k]*norm;
312   G4cout << " Total energy :  Edep + Eleak = " << 249       rmsLAbs  = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs*MeanLAbs));
313          << G4BestUnit(rmsEtotal, "Energy") << << 250 
314                                                << 251       //print
315   G4cout << "--------------------------------- << 252       //
                                                   >> 253       G4cout
                                                   >> 254        << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName() << ": "
                                                   >> 255        << std::setprecision(5)
                                                   >> 256        << std::setw(6) << G4BestUnit(MeanEAbs,"Energy") << " :  "
                                                   >> 257        << std::setprecision(4)
                                                   >> 258        << std::setw(5) << G4BestUnit( rmsEAbs,"Energy")  
                                                   >> 259        << std::setw(10) << resolution  << " +- " 
                                                   >> 260        << std::setw(5) << rmsres << " %"
                                                   >> 261        << std::setprecision(3)
                                                   >> 262        << std::setw(10) << G4BestUnit(MeanLAbs,"Length")  << " +- "
                                                   >> 263        << std::setw(4) << G4BestUnit( rmsLAbs,"Length")
                                                   >> 264        << G4endl;
                                                   >> 265     }
                                                   >> 266   G4cout << "\n------------------------------------------------------------\n";
316                                                   267 
317   G4cout << " Beam particle " << fParticle->Ge << 268   G4cout << " Beam particle " 
318          << "  E = " << G4BestUnit(beamEnergy, << 269          << fParticle->GetParticleName()
319   G4cout << " Mean number of gamma       " <<  << 270          << "  E = " << G4BestUnit(beamEnergy,"Energy") << G4endl;
320   G4cout << " Mean number of e-          " <<  << 271   G4cout << " Mean number of gamma       " << (G4double)fN_gamma*norm << G4endl;
321   G4cout << " Mean number of e+          " <<  << 272   G4cout << " Mean number of e-          " << (G4double)fN_elec*norm << G4endl;
322   G4cout << std::setprecision(6) << " Mean num << 273   G4cout << " Mean number of e+          " << (G4double)fN_pos*norm << G4endl;
                                                   >> 274   G4cout << std::setprecision(6)
                                                   >> 275          << " Mean number of charged steps  " << fChargedStep << G4endl;
323   G4cout << " Mean number of neutral steps  "     276   G4cout << " Mean number of neutral steps  " << fNeutralStep << G4endl;
324   G4cout << "---------------------------------    277   G4cout << "------------------------------------------------------------\n";
325                                                << 278   
326   // Energy flow                               << 279   //Energy flow
327   //                                              280   //
328   G4AnalysisManager* analysis = G4AnalysisMana    281   G4AnalysisManager* analysis = G4AnalysisManager::Instance();
329   G4int Idmax = (fDetector->GetNbOfLayers()) * << 282   G4int Idmax = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor());
330   for (G4int Id = 1; Id <= Idmax + 1; Id++) {  << 283   for (G4int Id=1; Id<=Idmax+1; Id++) {
331     analysis->FillH1(2 * kMaxAbsor + 1, (G4dou << 284     analysis->FillH1(2*kMaxAbsor+1, (G4double)Id, fEnergyFlow[Id]);
332     analysis->FillH1(2 * kMaxAbsor + 2, (G4dou << 285     analysis->FillH1(2*kMaxAbsor+2, (G4double)Id, fLateralEleak[Id]);
333   }                                               286   }
334                                                << 287   
335   // Energy deposit from energy flow balance   << 288   //Energy deposit from energy flow balance
336   //                                              289   //
337   G4double EdepTot[kMaxAbsor];                    290   G4double EdepTot[kMaxAbsor];
338   for (G4int k = 0; k < kMaxAbsor; k++)        << 291   for (G4int k=0; k<kMaxAbsor; k++) EdepTot[k] = 0.;
339     EdepTot[k] = 0.;                           << 292   
340                                                << 
341   G4int nbOfAbsor = fDetector->GetNbOfAbsor();    293   G4int nbOfAbsor = fDetector->GetNbOfAbsor();
342   for (G4int Id = 1; Id <= Idmax; Id++) {      << 294   for (G4int Id=1; Id<=Idmax; Id++) {
343     G4int iAbsor = Id % nbOfAbsor;             << 295     G4int iAbsor = Id%nbOfAbsor; if (iAbsor==0) iAbsor = nbOfAbsor;
344     if (iAbsor == 0) iAbsor = nbOfAbsor;       << 296     EdepTot[iAbsor] += (fEnergyFlow[Id]-fEnergyFlow[Id+1]-fLateralEleak[Id]);
345     EdepTot[iAbsor] += (fEnergyFlow[Id] - fEne << 297   }
346   }                                            << 298   
347                                                << 299   G4cout << std::setprecision(3)
348   G4cout << std::setprecision(3) << "\n Energy << 300          << "\n Energy deposition from Energy flow balance : \n"
349          << std::setw(10) << "  material \t To    301          << std::setw(10) << "  material \t Total Edep \n \n";
350   G4cout.precision(6);                            302   G4cout.precision(6);
351                                                << 303   
352   for (G4int k = 1; k <= nbOfAbsor; k++) {     << 304   for (G4int k=1; k<=nbOfAbsor; k++) {
353     EdepTot[k] *= norm;                        << 305     EdepTot [k] *= norm;
354     G4cout << std::setw(10) << fDetector->GetA    306     G4cout << std::setw(10) << fDetector->GetAbsorMaterial(k)->GetName() << ":"
355            << "\t " << G4BestUnit(EdepTot[k],  << 307            << "\t " << G4BestUnit(EdepTot [k],"Energy") << "\n";
356   }                                               308   }
357                                                << 309   
358   G4cout << "\n------------------------------- << 310   G4cout << "\n------------------------------------------------------------\n" 
                                                   >> 311          << G4endl;
359                                                   312 
360   // Acceptance                                   313   // Acceptance
361   EmAcceptance acc;                               314   EmAcceptance acc;
362   G4bool isStarted = false;                       315   G4bool isStarted = false;
363   for (G4int j = 1; j <= fDetector->GetNbOfAbs << 316   for (G4int j=1; j<=fDetector->GetNbOfAbsor(); j++) {
364     if (fLimittrue[j] < DBL_MAX) {                317     if (fLimittrue[j] < DBL_MAX) {
365       if (!isStarted) {                           318       if (!isStarted) {
366         acc.BeginOfAcceptance("Sampling Calori << 319         acc.BeginOfAcceptance("Sampling Calorimeter",nEvt);
367         isStarted = true;                         320         isStarted = true;
368       }                                           321       }
369       MeanEAbs = fSumEAbs[j];                     322       MeanEAbs = fSumEAbs[j];
370       rmsEAbs = fSum2EAbs[j];                  << 323       rmsEAbs  = fSum2EAbs[j];
371       G4String mat = fDetector->GetAbsorMateri    324       G4String mat = fDetector->GetAbsorMaterial(j)->GetName();
372       acc.EmAcceptanceGauss("Edep" + mat, nEvt << 325       acc.EmAcceptanceGauss("Edep"+mat, nEvt, MeanEAbs,
373       acc.EmAcceptanceGauss("Erms" + mat, nEvt << 326                              fEdeptrue[j], fRmstrue[j], fLimittrue[j]);
374                             2.0 * fLimittrue[j << 327       acc.EmAcceptanceGauss("Erms"+mat, nEvt, rmsEAbs,
                                                   >> 328                              fRmstrue[j], fRmstrue[j], 2.0*fLimittrue[j]);
375     }                                             329     }
376   }                                               330   }
377   if (isStarted) acc.EndOfAcceptance();        << 331   if(isStarted) acc.EndOfAcceptance();
378                                                   332 
379   // normalize histograms                      << 333   //normalize histograms
380   //                                              334   //
381   for (G4int ih = kMaxAbsor + 1; ih < kMaxHist << 335   for (G4int ih = kMaxAbsor+1; ih < kMaxHisto; ih++) {
382     analysis->ScaleH1(ih, norm / MeV);         << 336     analysis->ScaleH1(ih,norm/MeV);
383   }                                               337   }
384                                                << 338   
385   G4cout.setf(mode, std::ios::floatfield);     << 339   G4cout.setf(mode,std::ios::floatfield);
386   G4cout.precision(prec);                         340   G4cout.precision(prec);
387 }                                                 341 }
388                                                   342 
389 //....oooOO0OOooo........oooOO0OOooo........oo    343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
390                                                   344 
391 void Run::SetEdepAndRMS(G4int i, G4double edep    345 void Run::SetEdepAndRMS(G4int i, G4double edep, G4double rms, G4double lim)
392 {                                                 346 {
393   if (i >= 0 && i < kMaxAbsor) {               << 347   if (i>=0 && i<kMaxAbsor) {
394     fEdeptrue[i] = edep;                       << 348     fEdeptrue [i] = edep;
395     fRmstrue[i] = rms;                         << 349     fRmstrue  [i] = rms;
396     fLimittrue[i] = lim;                          350     fLimittrue[i] = lim;
397   }                                               351   }
398 }                                                 352 }
399                                                   353 
400 //....oooOO0OOooo........oooOO0OOooo........oo    354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
401                                                   355 
402 void Run::SetApplyLimit(G4bool val)               356 void Run::SetApplyLimit(G4bool val)
403 {                                                 357 {
404   fApplyLimit = val;                              358   fApplyLimit = val;
405 }                                                 359 }
406                                                   360 
407 //....oooOO0OOooo........oooOO0OOooo........oo    361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
408                                                   362