Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm18/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/TestEm18/src/RunAction.cc (Version 11.3.0) and /examples/extended/electromagnetic/TestEm18/src/RunAction.cc (Version 9.0.p1)


  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/TestEm18/src/RunActi <<  26 // $Id: RunAction.cc,v 1.2 2007/02/16 11:59:47 maire Exp $
 27 /// \brief Implementation of the RunAction cla <<  27 // GEANT4 tag $Name: geant4-08-03-patch-01 $
 28 //                                             << 
 29 //                                                 28 //
 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                                                << 
 35 #include "DetectorConstruction.hh"                 33 #include "DetectorConstruction.hh"
 36 #include "HistoManager.hh"                     << 
 37 #include "PrimaryGeneratorAction.hh"               34 #include "PrimaryGeneratorAction.hh"
                                                   >>  35 #include "HistoManager.hh"
 38                                                    36 
 39 #include "G4EmCalculator.hh"                   << 
 40 #include "G4Run.hh"                                37 #include "G4Run.hh"
                                                   >>  38 #include "G4RunManager.hh"
 41 #include "G4UnitsTable.hh"                         39 #include "G4UnitsTable.hh"
 42 #include "Randomize.hh"                        <<  40 #include "G4EmCalculator.hh"
 43                                                    41 
                                                   >>  42 #include "Randomize.hh"
 44 #include <iomanip>                                 43 #include <iomanip>
 45                                                    44 
 46 //....oooOO0OOooo........oooOO0OOooo........oo     45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 47                                                    46 
 48 RunAction::RunAction(DetectorConstruction* det <<  47 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin,
 49   : fDetector(det), fPrimary(kin)              <<  48                      HistoManager* histo)
 50 {                                              <<  49 :detector(det), primary(kin), histoManager(histo)
 51   fHistoManager = new HistoManager();          <<  50 { }
 52 }                                              << 
 53                                                    51 
 54 //....oooOO0OOooo........oooOO0OOooo........oo     52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 55                                                    53 
 56 RunAction::~RunAction()                            54 RunAction::~RunAction()
 57 {                                              <<  55 { }
 58   delete fHistoManager;                        << 
 59 }                                              << 
 60                                                    56 
 61 //....oooOO0OOooo........oooOO0OOooo........oo     57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 62                                                    58 
 63 void RunAction::BeginOfRunAction(const G4Run*) <<  59 void RunAction::BeginOfRunAction(const G4Run* run)
 64 {                                                  60 {
 65   // initialisation                            <<  61   G4cout << "### Run " << run->GetRunID() << " start." << G4endl;
 66   //                                           << 
 67   fNbSteps = 0;                                << 
 68   fTrackLength = 0.;                           << 
 69   fStepMin = DBL_MAX;                          << 
 70   fStepMax = 0.;                               << 
 71                                                << 
 72   fEdepPrimary = fEdepSecondary = fEdepTotal = << 
 73   fEdepPrimMin = fEdepSecMin = fEdepTotMin = D << 
 74   fEdepPrimMax = fEdepSecMax = fEdepTotMax = 0 << 
 75                                                << 
 76   fEnergyTransfered = 0.;                      << 
 77   fEtransfMin = DBL_MAX;                       << 
 78   fEtransfMax = 0.;                            << 
 79                                                << 
 80   fEnergyLost = 0.;                            << 
 81   fElostMin = DBL_MAX;                         << 
 82   fElostMax = 0.;                              << 
 83                                                << 
 84   fEnergyBalance = 0.;                         << 
 85   fEbalMin = DBL_MAX;                          << 
 86   fEbalMax = 0.;                               << 
 87                                                    62 
 88   // histograms                                <<  63   //initialisation
 89   //                                               64   //
 90   G4AnalysisManager* analysisManager = G4Analy <<  65   energyDeposit = 0.;
 91   if (analysisManager->IsActive()) {           <<  66   
 92     analysisManager->OpenFile();               <<  67   nbCharged = nbNeutral = 0;
 93   }                                            <<  68   energyCharged = energyNeutral = 0.;  
                                                   >>  69   emin[0] = emin[1] = DBL_MAX;
                                                   >>  70   emax[0] = emax[1] = 0.;    
                                                   >>  71     
                                                   >>  72   nbSteps = 0;
                                                   >>  73   trackLength = 0.; 
 94                                                    74 
 95   // show Rndm status                          <<  75   histoManager->book();
 96   CLHEP::HepRandom::showEngineStatus();        << 
 97 }                                              << 
 98                                                << 
 99 //....oooOO0OOooo........oooOO0OOooo........oo << 
100                                                << 
101 void RunAction::CountProcesses(G4String procNa << 
102 {                                              << 
103   std::map<G4String, G4int>::iterator it = fPr << 
104   if (it == fProcCounter.end()) {              << 
105     fProcCounter[procName] = 1;                << 
106   }                                            << 
107   else {                                       << 
108     fProcCounter[procName]++;                  << 
109   }                                            << 
110 }                                              << 
111                                                << 
112 //....oooOO0OOooo........oooOO0OOooo........oo << 
113                                                << 
114 void RunAction::TrackLength(G4double step)     << 
115 {                                              << 
116   fTrackLength += step;                        << 
117   fNbSteps++;                                  << 
118   if (step < fStepMin) fStepMin = step;        << 
119   if (step > fStepMax) fStepMax = step;        << 
120 }                                              << 
121                                                    76 
122 //....oooOO0OOooo........oooOO0OOooo........oo <<  77   // do not save Rndm status
123                                                <<  78   G4RunManager::GetRunManager()->SetRandomNumberStore(false);
124 void RunAction::EnergyDeposited(G4double edepP <<  79   CLHEP::HepRandom::showEngineStatus();
125 {                                              << 
126   fEdepPrimary += edepPrim;                    << 
127   if (edepPrim < fEdepPrimMin) fEdepPrimMin =  << 
128   if (edepPrim > fEdepPrimMax) fEdepPrimMax =  << 
129                                                << 
130   fEdepSecondary += edepSecond;                << 
131   if (edepSecond < fEdepSecMin) fEdepSecMin =  << 
132   if (edepSecond > fEdepSecMax) fEdepSecMax =  << 
133 }                                              << 
134                                                << 
135 //....oooOO0OOooo........oooOO0OOooo........oo << 
136                                                << 
137 void RunAction::EnergyTransferedByProcess(G4St << 
138 {                                              << 
139   std::map<G4String, MinMaxData>::iterator it  << 
140   if (it == fEtransfByProcess.end()) {         << 
141     fEtransfByProcess[process] = MinMaxData(1, << 
142   }                                            << 
143   else {                                       << 
144     MinMaxData& data = it->second;             << 
145     data.fCount++;                             << 
146     data.fVsum += energy;                      << 
147     // update min max                          << 
148     G4double emin = data.fVmin;                << 
149     if (energy < emin) data.fVmin = energy;    << 
150     G4double emax = data.fVmax;                << 
151     if (energy > emax) data.fVmax = energy;    << 
152   }                                            << 
153 }                                              << 
154                                                << 
155 //....oooOO0OOooo........oooOO0OOooo........oo << 
156                                                << 
157 void RunAction::EnergyTransfered(G4double ener << 
158 {                                              << 
159   fEnergyTransfered += energy;                 << 
160   if (energy < fEtransfMin) fEtransfMin = ener << 
161   if (energy > fEtransfMax) fEtransfMax = ener << 
162 }                                              << 
163                                                << 
164 //....oooOO0OOooo........oooOO0OOooo........oo << 
165                                                << 
166 void RunAction::TotalEnergyLost(G4double energ << 
167 {                                              << 
168   fEnergyLost += energy;                       << 
169   if (energy < fElostMin) fElostMin = energy;  << 
170   if (energy > fElostMax) fElostMax = energy;  << 
171 }                                              << 
172                                                << 
173 //....oooOO0OOooo........oooOO0OOooo........oo << 
174                                                << 
175 void RunAction::EnergyBalance(G4double energy) << 
176 {                                              << 
177   fEnergyBalance += energy;                    << 
178   if (energy < fEbalMin) fEbalMin = energy;    << 
179   if (energy > fEbalMax) fEbalMax = energy;    << 
180 }                                              << 
181                                                << 
182 //....oooOO0OOooo........oooOO0OOooo........oo << 
183                                                << 
184 void RunAction::TotalEnergyDeposit(G4double en << 
185 {                                              << 
186   fEdepTotal += energy;                        << 
187   if (energy < fEdepTotMin) fEdepTotMin = ener << 
188   if (energy > fEdepTotMax) fEdepTotMax = ener << 
189 }                                              << 
190                                                << 
191 //....oooOO0OOooo........oooOO0OOooo........oo << 
192                                                << 
193 void RunAction::EnergySpectrumOfSecondaries(G4 << 
194 {                                              << 
195   std::map<G4String, MinMaxData>::iterator it  << 
196   if (it == fEkinOfSecondaries.end()) {        << 
197     fEkinOfSecondaries[particle] = MinMaxData( << 
198   }                                            << 
199   else {                                       << 
200     MinMaxData& data = it->second;             << 
201     data.fCount++;                             << 
202     data.fVsum += energy;                      << 
203     // update min max                          << 
204     G4double emin = data.fVmin;                << 
205     if (energy < emin) data.fVmin = energy;    << 
206     G4double emax = data.fVmax;                << 
207     if (energy > emax) data.fVmax = energy;    << 
208   }                                            << 
209 }                                                  80 }
210                                                    81 
211 //....oooOO0OOooo........oooOO0OOooo........oo     82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
212                                                    83 
213 void RunAction::EndOfRunAction(const G4Run* aR     84 void RunAction::EndOfRunAction(const G4Run* aRun)
214 {                                                  85 {
215   G4int nbEvents = aRun->GetNumberOfEvent();       86   G4int nbEvents = aRun->GetNumberOfEvent();
216   if (nbEvents == 0) return;                       87   if (nbEvents == 0) return;
217                                                <<  88   
218   G4Material* material = fDetector->GetMateria <<  89   G4Material* material = detector->GetMaterial();
219   G4double length = fDetector->GetSize();      <<  90   G4double length  = detector->GetSize();
220   G4double density = material->GetDensity();       91   G4double density = material->GetDensity();
221                                                <<  92    
222   G4ParticleDefinition* particle = fPrimary->G <<  93   G4ParticleDefinition* particle = primary->GetParticleGun()
                                                   >>  94                                           ->GetParticleDefinition();
223   G4String partName = particle->GetParticleNam     95   G4String partName = particle->GetParticleName();
224   G4double ePrimary = fPrimary->GetParticleGun <<  96   G4double eprimary = primary->GetParticleGun()->GetParticleEnergy();
225                                                <<  97   
226   G4int prec = G4cout.precision(3);                98   G4int prec = G4cout.precision(3);
227   G4cout << "\n ======================== run s <<  99   G4cout << "\n ======================== run summary ======================\n";  
228   G4cout << "\n The run was " << nbEvents << "    100   G4cout << "\n The run was " << nbEvents << " " << partName << " of "
229          << G4BestUnit(ePrimary, "Energy") <<  << 101          << G4BestUnit(eprimary,"Energy") << " through " 
230          << material->GetName() << " (density: << 102    << G4BestUnit(length,"Length") << " of "
                                                   >> 103    << material->GetName() << " (density: " 
                                                   >> 104    << G4BestUnit(density,"Volumic Mass") << ")";
                                                   >> 105   G4cout << "\n ===========================================================\n";
231   G4cout << G4endl;                               106   G4cout << G4endl;
232                                                << 107   
233   if (particle->GetPDGCharge() == 0.) return;     108   if (particle->GetPDGCharge() == 0.) return;
                                                   >> 109    
                                                   >> 110   G4cout.precision(5);
                                                   >> 111   
                                                   >> 112   //track length
                                                   >> 113   //
                                                   >> 114   G4double trackLPerEvent = trackLength/nbEvents;
                                                   >> 115   G4double nbStepPerEvent = double(nbSteps)/nbEvents;
                                                   >> 116   G4double stepSize = trackLength/nbSteps;
                                                   >> 117   
                                                   >> 118   G4cout 
                                                   >> 119     << "\n trackLength= " 
                                                   >> 120     << G4BestUnit(trackLPerEvent, "Length")
                                                   >> 121     << "\t nb of steps= " << nbStepPerEvent
                                                   >> 122     << "  stepSize= " << G4BestUnit(stepSize, "Length")
                                                   >> 123     << G4endl;
                                                   >> 124       
                                                   >> 125   //charged secondaries (ionization, direct pair production)
                                                   >> 126   //
                                                   >> 127   G4double energyPerEvent = energyCharged/nbEvents;
                                                   >> 128   G4double nbPerEvent = double(nbCharged)/nbEvents;
                                                   >> 129   G4double meanEkin = 0.;
                                                   >> 130   if (nbCharged) meanEkin = energyCharged/nbCharged;
                                                   >> 131   
                                                   >> 132   G4cout 
                                                   >> 133     << "\n d-rays  : eLoss/primary= " 
                                                   >> 134     << G4BestUnit(energyPerEvent, "Energy")
                                                   >> 135     << "\t  nb of d-rays= " << nbPerEvent
                                                   >> 136     << "  <Tkin>= " << G4BestUnit(meanEkin, "Energy")
                                                   >> 137     << "  Tmin= "   << G4BestUnit(emin[0],  "Energy")
                                                   >> 138     << "  Tmax= "   << G4BestUnit(emax[0],  "Energy")
                                                   >> 139     << G4endl;
                                                   >> 140          
                                                   >> 141   //neutral secondaries (bremsstrahlung)
                                                   >> 142   //
                                                   >> 143   energyPerEvent = energyNeutral/nbEvents;
                                                   >> 144   nbPerEvent = double(nbNeutral)/nbEvents;
                                                   >> 145   meanEkin = 0.;
                                                   >> 146   if (nbNeutral) meanEkin = energyNeutral/nbNeutral;
                                                   >> 147   
                                                   >> 148   G4cout 
                                                   >> 149     << "\n brems   : eLoss/primary= " 
                                                   >> 150     << G4BestUnit(energyPerEvent, "Energy")
                                                   >> 151     << "\t  nb of gammas= " << nbPerEvent
                                                   >> 152     << "  <Tkin>= " << G4BestUnit(meanEkin, "Energy")
                                                   >> 153     << "  Tmin= "   << G4BestUnit(emin[1],  "Energy")
                                                   >> 154     << "  Tmax= "   << G4BestUnit(emax[1],  "Energy")
                                                   >> 155     << G4endl;
                                                   >> 156     
234                                                   157 
235   G4cout.precision(4);                         << 158   G4EmCalculator emCal;
236                                                << 159   
237   // frequency of processes                    << 160   //local energy deposit
238   //                                           << 
239   G4cout << "\n Process defining step :" << G4 << 
240   G4int index = 0;                             << 
241   for (const auto& procCounter : fProcCounter) << 
242     G4String procName = procCounter.first;     << 
243     G4int count = procCounter.second;          << 
244     G4String space = " ";                      << 
245     if (++index % 4 == 0) space = "\n";        << 
246     G4cout << " " << std::setw(15) << procName << 
247   }                                            << 
248   G4cout << G4endl;                            << 
249                                                << 
250   // track length                              << 
251   //                                           << 
252   G4double trackLPerEvent = fTrackLength / nbE << 
253   G4double nbStepPerEvent = double(fNbSteps) / << 
254   G4double stepSize = fTrackLength / fNbSteps; << 
255                                                << 
256   G4cout << "\n TrackLength = " << G4BestUnit( << 
257          << "  nb of steps = " << nbStepPerEve << 
258          << "  stepSize = " << G4BestUnit(step << 
259          << G4BestUnit(fStepMin, "Length") <<  << 
260          << G4endl;                            << 
261                                                << 
262   // continuous energy deposited by primary tr << 
263   //                                              161   //
264   G4double energyPerEvent = fEdepPrimary / nbE << 162   energyPerEvent = energyDeposit/nbEvents;
265                                                << 
266   G4cout << "\n Energy continuously deposited  << 
267          << " (restricted dE/dx)  dE1 = " << G << 
268          << G4BestUnit(fEdepPrimMin, "Energy") << 
269          << ")" << G4endl;                     << 
270                                                << 
271   // eveluation of dE1 from reading restricted << 
272   //                                              163   //
273   G4EmCalculator emCal;                        << 164   G4double r0  = emCal.GetRangeFromRestricteDEDX(eprimary,particle,material);  
274                                                << 
275   G4double r0 = emCal.GetRangeFromRestricteDED << 
276   G4double r1 = r0 - trackLPerEvent;              165   G4double r1 = r0 - trackLPerEvent;
277   G4double etry = ePrimary - energyPerEvent;   << 166   G4double etry = eprimary - energyPerEvent;  
278   G4double efinal = 0.;                           167   G4double efinal = 0.;
279   if (r1 > 0.) efinal = GetEnergyFromRestricte << 168   if (r1 > 0.) efinal = GetEnergyFromRestrictedRange(r1,particle,material,etry);
280   G4double dEtable = ePrimary - efinal;        << 169   G4double dEtable = eprimary - efinal;
281   G4double ratio = 0.;                            170   G4double ratio = 0.;
282   if (dEtable > 0.) ratio = energyPerEvent / d << 171   if (dEtable > 0.) ratio = energyPerEvent/dEtable;
283                                                << 172     
284   G4cout << "\n Evaluation of dE1 from reading << 173   G4cout 
285          << G4BestUnit(dEtable, "Energy") << " << 174     << "\n deposit : eLoss/primary= " 
286                                                << 175     << G4BestUnit(energyPerEvent, "Energy")
287   // energy transfered to secondary particles  << 176     << "\t <dEcut > table= " 
                                                   >> 177     << G4BestUnit(dEtable, "Energy")
                                                   >> 178     << "   ---> simul/reference= " << ratio           
                                                   >> 179     << G4endl;
                                                   >> 180     
                                                   >> 181   //total energy transferred
288   //                                              182   //
289   G4cout << "\n Energy transfered to secondary << 183   G4double energyTotal = energyDeposit + energyCharged + energyNeutral;
290   std::map<G4String, MinMaxData>::iterator it1 << 184   energyPerEvent = energyTotal/nbEvents;
291   for (it1 = fEtransfByProcess.begin(); it1 != << 
292     G4String name = it1->first;                << 
293     MinMaxData data = it1->second;             << 
294     energyPerEvent = data.fVsum / nbEvents;    << 
295     G4double eMin = data.fVmin;                << 
296     G4double eMax = data.fVmax;                << 
297                                                << 
298     G4cout << "  " << std::setw(17) << "due to << 
299            << G4BestUnit(energyPerEvent, "Ener << 
300            << G4BestUnit(eMax, "Energy") << ") << 
301   }                                            << 
302                                                << 
303   // total energy tranfered : dE3 = sum of dE2 << 
304   //                                           << 
305   energyPerEvent = fEnergyTransfered / nbEvent << 
306                                                << 
307   G4cout << "\n Total energy transfered to sec << 
308          << G4BestUnit(energyPerEvent, "Energy << 
309          << " --> " << G4BestUnit(fEtransfMax, << 
310                                                << 
311   // total energy lost by incident particle :  << 
312   //                                           << 
313   energyPerEvent = fEnergyLost / nbEvents;     << 
314                                                << 
315   G4cout << "\n Total energy lost by incident  << 
316          << G4BestUnit(energyPerEvent, "Energy << 
317          << " --> " << G4BestUnit(fElostMax, " << 
318                                                << 
319   // calcul of energy lost from energy balance << 
320   //                                           << 
321   energyPerEvent = fEnergyBalance / nbEvents;  << 
322                                                << 
323   G4cout << "\n calcul of dE4 from energy bala << 
324          << G4BestUnit(energyPerEvent, "Energy << 
325          << " --> " << G4BestUnit(fEbalMax, "E << 
326                                                << 
327   // eveluation of dE4 from reading full Range << 
328   //                                              185   //
329   r0 = emCal.GetCSDARange(ePrimary, particle,  << 186   r0  = emCal.GetCSDARange(eprimary,particle,material);  
330   r1 = r0 - trackLPerEvent;                       187   r1 = r0 - trackLPerEvent;
331   etry = ePrimary - energyPerEvent;            << 188   etry = eprimary - energyPerEvent;
332   efinal = 0.;                                    189   efinal = 0.;
333   if (r1 > 0.) efinal = GetEnergyFromCSDARange << 190   if (r1 > 0.) efinal = GetEnergyFromCSDARange(r1,particle,material,etry);
334   dEtable = ePrimary - efinal;                 << 191   dEtable = eprimary - efinal;
335   ratio = 0.;                                     192   ratio = 0.;
336   if (dEtable > 0.) ratio = energyPerEvent / d << 193   if (dEtable > 0.) ratio = energyPerEvent/dEtable;
337                                                << 194     
338   G4cout << "\n Evaluation of dE4 from reading << 195   G4cout 
339          << G4BestUnit(dEtable, "Energy") << " << 196     << "\n total   : eLoss/primary= " 
340                                                << 197     << G4BestUnit(energyPerEvent, "Energy")
341   // energy spectrum of secondary particles    << 198     << "\t <dEfull> table= " 
342   //                                           << 199     << G4BestUnit(dEtable, "Energy")
343   G4cout << "\n Energy spectrum of secondary p << 200     << "   ---> simul/reference= " << ratio           
344   std::map<G4String, MinMaxData>::iterator it2 << 201     << G4endl; 
345   for (it2 = fEkinOfSecondaries.begin(); it2 ! << 
346     G4String name = it2->first;                << 
347     MinMaxData data = it2->second;             << 
348     G4int count = data.fCount;                 << 
349     G4double eMean = data.fVsum / count;       << 
350     G4double eMin = data.fVmin;                << 
351     G4double eMax = data.fVmax;                << 
352                                                << 
353     G4cout << "  " << std::setw(13) << name << << 
354            << "  Emean = " << std::setw(6) <<  << 
355            << G4BestUnit(eMin, "Energy") << "  << 
356   }                                            << 
357   G4cout << G4endl;                            << 
358                                                << 
359   // continuous energy deposited by secondary  << 
360   //  (only if secondary particles are tracked << 
361   //                                           << 
362   if (fEdepSecondary > 0.) {                   << 
363     energyPerEvent = fEdepSecondary / nbEvents << 
364                                                << 
365     G4cout << "\n Energy continuously deposite << 
366            << " (restricted dE/dx)  dE5 = " << << 
367            << G4BestUnit(fEdepSecMin, "Energy" << 
368            << ")" << G4endl;                   << 
369                                                << 
370     // total energy deposited : dE6 = dE1 + dE << 
371     //                                         << 
372     energyPerEvent = fEdepTotal / nbEvents;    << 
373                                                << 
374     G4cout << "\n Total energy deposited : dE6 << 
375            << G4BestUnit(energyPerEvent, "Ener << 
376            << " --> " << G4BestUnit(fEdepTotMa << 
377            << G4endl;                          << 
378   }                                            << 
379                                                   202 
380   G4cout.precision(prec);                         203   G4cout.precision(prec);
381                                                   204 
382   // clear maps                                << 205   histoManager->save();
383   //                                           << 
384   fProcCounter.clear();                        << 
385   fEtransfByProcess.clear();                   << 
386   fEkinOfSecondaries.clear();                  << 
387                                                << 
388   // save histograms                           << 
389   G4AnalysisManager* analysisManager = G4Analy << 
390   if (analysisManager->IsActive()) {           << 
391     analysisManager->Write();                  << 
392     analysisManager->CloseFile();              << 
393   }                                            << 
394                                                   206 
395   // show Rndm status                             207   // show Rndm status
396   CLHEP::HepRandom::showEngineStatus();           208   CLHEP::HepRandom::showEngineStatus();
397 }                                                 209 }
398                                                   210 
399 //....oooOO0OOooo........oooOO0OOooo........oo    211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
400                                                   212 
401 G4double RunAction::GetEnergyFromRestrictedRan << 213 G4double RunAction::GetEnergyFromRestrictedRange(G4double range,
402                                                << 214             G4ParticleDefinition* particle, G4Material* material, G4double Etry)
403 {                                                 215 {
404   G4EmCalculator emCal;                           216   G4EmCalculator emCal;
405                                                << 217     
406   G4double Energy = Etry, dE = 0., dEdx;          218   G4double Energy = Etry, dE = 0., dEdx;
407   G4double r, dr;                                 219   G4double r, dr;
408   G4double err = 1., errmax = 0.00001;         << 220   G4double err  = 1., errmax = 0.00001;
409   G4int iter = 0, itermax = 10;                << 221   G4int    iter = 0 , itermax = 10;  
410   while (err > errmax && iter < itermax) {        222   while (err > errmax && iter < itermax) {
411     iter++;                                       223     iter++;
412     Energy -= dE;                                 224     Energy -= dE;
413     r = emCal.GetRangeFromRestricteDEDX(Energy << 225     r = emCal.GetRangeFromRestricteDEDX(Energy,particle,material);
414     dr = r - range;                            << 226     dr = r - range;          
415     dEdx = emCal.GetDEDX(Energy, particle, mat << 227     dEdx = emCal.GetDEDX(Energy,particle,material);
416     dE = dEdx * dr;                            << 228     dE = dEdx*dr;
417     err = std::abs(dE) / Energy;               << 229     err = std::abs(dE)/Energy;    
418   }                                               230   }
419   if (iter == itermax) {                          231   if (iter == itermax) {
420     G4cout << "\n  ---> warning: RunAction::Ge << 232     G4cout 
421            << "   Etry = " << G4BestUnit(Etry, << 233     << "\n  ---> warning: RunAction::GetEnergyFromRestRange() did not converge"
422            << "   Energy = " << G4BestUnit(Ene << 234     << "   Etry = " << G4BestUnit(Etry,"Energy")
423            << "   iter = " << iter << G4endl;  << 235     << "   Energy = " << G4BestUnit(Energy,"Energy")
424   }                                            << 236     << "   err = " << err
425                                                << 237     << "   iter = " << iter << G4endl;
426   return Energy;                               << 238   }  
                                                   >> 239    
                                                   >> 240   return Energy;   
427 }                                                 241 }
428                                                   242 
429 //....oooOO0OOooo........oooOO0OOooo........oo    243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
430                                                   244 
431 G4double RunAction::GetEnergyFromCSDARange(G4d << 245 G4double RunAction::GetEnergyFromCSDARange(G4double range,
432                                            G4M << 246             G4ParticleDefinition* particle, G4Material* material, G4double Etry)
433 {                                                 247 {
434   G4EmCalculator emCal;                           248   G4EmCalculator emCal;
435                                                << 249     
436   G4double Energy = Etry, dE = 0., dEdx;          250   G4double Energy = Etry, dE = 0., dEdx;
437   G4double r, dr;                                 251   G4double r, dr;
438   G4double err = 1., errmax = 0.00001;         << 252   G4double err  = 1., errmax = 0.00001;
439   G4int iter = 0, itermax = 10;                << 253   G4int    iter = 0 , itermax = 10;  
440   while (err > errmax && iter < itermax) {        254   while (err > errmax && iter < itermax) {
441     iter++;                                       255     iter++;
442     Energy -= dE;                                 256     Energy -= dE;
443     r = emCal.GetCSDARange(Energy, particle, m << 257     r = emCal.GetCSDARange(Energy,particle,material);
444     dr = r - range;                            << 258     dr = r - range;          
445     dEdx = emCal.ComputeTotalDEDX(Energy, part << 259     dEdx = emCal.ComputeTotalDEDX(Energy,particle,material);
446     dE = dEdx * dr;                            << 260     dE = dEdx*dr;
447     err = std::abs(dE) / Energy;               << 261     err = std::abs(dE)/Energy;
448   }                                               262   }
449   if (iter == itermax) {                          263   if (iter == itermax) {
450     G4cout << "\n  ---> warning: RunAction::Ge << 264     G4cout 
451            << "   Etry = " << G4BestUnit(Etry, << 265     << "\n  ---> warning: RunAction::GetEnergyFromCSDARange() did not converge"
452            << "   Energy = " << G4BestUnit(Ene << 266     << "   Etry = " << G4BestUnit(Etry,"Energy")
453            << "   iter = " << iter << G4endl;  << 267     << "   Energy = " << G4BestUnit(Energy,"Energy")
454   }                                            << 268     << "   err = " << err
455                                                << 269     << "   iter = " << iter << G4endl;
456   return Energy;                               << 270   }  
                                                   >> 271    
                                                   >> 272   return Energy;   
457 }                                                 273 }
458                                                   274 
459 //....oooOO0OOooo........oooOO0OOooo........oo    275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
460                                                   276