Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm11/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/TestEm11/src/RunAction.cc (Version 11.3.0) and /examples/extended/electromagnetic/TestEm11/src/RunAction.cc (Version 9.6.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/TestEm11/src/RunActi     26 /// \file electromagnetic/TestEm11/src/RunAction.cc
 27 /// \brief Implementation of the RunAction cla     27 /// \brief Implementation of the RunAction class
 28 //                                                 28 //
 29 //                                             <<  29 // $Id$
                                                   >>  30 // 
 30 //....oooOO0OOooo........oooOO0OOooo........oo     31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oo     32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32                                                    33 
 33 #include "RunAction.hh"                            34 #include "RunAction.hh"
 34                                                << 
 35 #include "HistoManager.hh"                     << 
 36 #include "PhysicsList.hh"                          35 #include "PhysicsList.hh"
 37 #include "PrimaryGeneratorAction.hh"           << 
 38 #include "Run.hh"                              << 
 39 #include "StepMax.hh"                              36 #include "StepMax.hh"
                                                   >>  37 #include "PrimaryGeneratorAction.hh"
                                                   >>  38 #include "HistoManager.hh"
 40                                                    39 
 41 #include "G4EmCalculator.hh"                   << 
 42 #include "G4EmParameters.hh"                   << 
 43 #include "G4Run.hh"                                40 #include "G4Run.hh"
 44 #include "G4SystemOfUnits.hh"                  <<  41 #include "G4RunManager.hh"
 45 #include "G4UnitsTable.hh"                         42 #include "G4UnitsTable.hh"
                                                   >>  43 #include "G4EmCalculator.hh"
                                                   >>  44 
 46 #include "Randomize.hh"                            45 #include "Randomize.hh"
                                                   >>  46 #include "G4SystemOfUnits.hh"
                                                   >>  47 #include <iomanip>
 47                                                    48 
 48 //....oooOO0OOooo........oooOO0OOooo........oo     49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 49                                                    50 
 50 RunAction::RunAction(DetectorConstruction* det <<  51 RunAction::RunAction(DetectorConstruction* det, PhysicsList* phys,
 51   : fDetector(det), fPhysics(phys), fPrimary(k <<  52                      PrimaryGeneratorAction* kin)
                                                   >>  53 :fDetector(det),fPhysics(phys),fKinematic(kin)
 52 {                                                  54 {
 53   // Book predefined histograms                    55   // Book predefined histograms
 54   fHistoManager = new HistoManager();              56   fHistoManager = new HistoManager();
 55 }                                                  57 }
 56                                                    58 
 57 //....oooOO0OOooo........oooOO0OOooo........oo     59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 58                                                    60 
 59 RunAction::~RunAction()                            61 RunAction::~RunAction()
 60 {                                                  62 {
 61   delete fHistoManager;                            63   delete fHistoManager;
 62 }                                                  64 }
 63                                                    65 
 64 //....oooOO0OOooo........oooOO0OOooo........oo     66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 65                                                    67 
 66 G4Run* RunAction::GenerateRun()                <<  68 void RunAction::BeginOfRunAction(const G4Run* aRun)
 67 {                                              <<  69 {  
 68   fRun = new Run(fDetector);                   <<  70   G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
 69   return fRun;                                 <<  71   
 70 }                                              <<  72   // save Rndm status
 71                                                <<  73   ////G4RunManager::GetRunManager()->SetRandomNumberStore(true);
 72 //....oooOO0OOooo........oooOO0OOooo........oo <<  74   CLHEP::HepRandom::showEngineStatus();
 73                                                <<  75   
 74 void RunAction::BeginOfRunAction(const G4Run*) <<  76   //initialize total energy deposit
 75 {                                              <<  77   //
 76   // show Rndm status                          <<  78   fEdeposit = fEdeposit2 = 0.; 
 77   if (isMaster) {                              <<  79    
 78     G4Random::showEngineStatus();              <<  80   //initialize track legth of primary
 79     G4EmParameters::Instance()->Dump();        <<  81   //
 80   }                                            <<  82   fTrackLen = fTrackLen2 = 0.;
 81                                                <<  83     
 82   // keep run condition                        <<  84   //initialize projected range
 83   if (fPrimary) {                              <<  85   //
 84     G4ParticleDefinition* particle = fPrimary- <<  86   fProjRange = fProjRange2 = 0.;
 85     G4double energy = fPrimary->GetParticleGun <<  87     
 86     fRun->SetPrimary(particle, energy);        <<  88   //initialize mean step size
                                                   >>  89   //
                                                   >>  90   fNbOfSteps = fNbOfSteps2 = 0;  fStepSize = fStepSize2 = 0.;
                                                   >>  91   
                                                   >>  92   //initialize track status
                                                   >>  93   //
                                                   >>  94   fStatus[0] = fStatus[1] = fStatus[2] = 0;
                                                   >>  95   
                                                   >>  96   //normalized binwidth
                                                   >>  97   //
                                                   >>  98   for (G4int i=0; i< MaxAbsor; i++) {
                                                   >>  99   fCsdaRange[i] = fXfrontNorm[i] = 0.;
 87   }                                               100   }
 88                                                << 101     
 89   // histograms                                << 102   //histograms
 90   //                                              103   //
 91   G4AnalysisManager* analysisManager = G4Analy    104   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
 92   if (analysisManager->IsActive()) {           << 105   if ( analysisManager->IsActive() ) {
 93     analysisManager->OpenFile();                  106     analysisManager->OpenFile();
 94   }                                            << 107   } 
 95                                                << 108     
 96   if (!fPrimary) return;                       << 109   //set StepMax from histos 1 and 8
 97                                                << 110   //e
 98   // set StepMax from histos 1 and 8           << 
 99   //                                           << 
100   G4double stepMax = DBL_MAX;                     111   G4double stepMax = DBL_MAX;
101   G4int ih = 1;                                   112   G4int ih = 1;
102   if (analysisManager->GetH1Activation(ih)) {  << 113   if (analysisManager->GetH1Activation(ih))
103     stepMax = analysisManager->GetH1Width(ih)  << 114                                     stepMax = analysisManager->GetH1Width(ih);
104   }                                            << 115   //
105                                                << 
106   ih = 8;                                         116   ih = 8;
107   G4ParticleDefinition* particle = fPrimary->G << 117   G4ParticleDefinition* particle = fKinematic->GetParticleGun()
                                                   >> 118                                           ->GetParticleDefinition();
108   if (particle->GetPDGCharge() != 0.) {           119   if (particle->GetPDGCharge() != 0.) {
109     G4double width = analysisManager->GetH1Wid    120     G4double width = analysisManager->GetH1Width(ih);
110     if (width == 0.) width = 1.;               << 121     if (width == 0.) width = 1.;                                               
111     G4EmCalculator emCalculator;                  122     G4EmCalculator emCalculator;
112     G4double energy = fPrimary->GetParticleGun << 123     G4double energy = fKinematic->GetParticleGun()->GetParticleEnergy();
113     G4int nbOfAbsor = fDetector->GetNbOfAbsor( << 124     G4int NbOfAbsor = fDetector->GetNbOfAbsor();
114     for (G4int i = 1; i <= nbOfAbsor; i++) {   << 125     for (G4int i=1; i<= NbOfAbsor; i++) {
115       G4Material* material = fDetector->GetAbs << 126       G4Material* material = fDetector->GetAbsorMaterial(i);  
116       G4double newCsdaRange = emCalculator.Get << 127       fCsdaRange[i] = emCalculator.GetCSDARange(energy,particle,material);
117       fRun->SetCsdaRange(i, newCsdaRange);     << 128       if (analysisManager->GetH1Activation(ih))
118       if (analysisManager->GetH1Activation(ih) << 129         stepMax = std::min(stepMax, width*fCsdaRange[i]);
119       if (i > 1) {                             << 130       if (i>1) {
120         G4double thickness = fDetector->GetAbs << 131         G4double thickness = fDetector->GetAbsorThickness(i-1);
121         G4double xfrontNorm = fRun->GetXfrontN << 132         fXfrontNorm[i] = fXfrontNorm[i-1] + thickness/fCsdaRange[i-1];
122         G4double csdaRange = fRun->GetCsdaRang << 133       }                              
123         G4double newXfrontNorm = xfrontNorm +  << 134     }        
124         fRun->SetXfrontNorm(i, newXfrontNorm); << 135   }     
125       }                                        << 136   fPhysics->GetStepMaxProcess()->SetMaxStep2(stepMax);           
126     }                                          << 
127   }                                            << 
128   fPhysics->GetStepMaxProcess()->SetMaxStep2(s << 
129 }                                                 137 }
130                                                   138 
131 //....oooOO0OOooo........oooOO0OOooo........oo    139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132                                                   140 
133 void RunAction::EndOfRunAction(const G4Run*)   << 141 void RunAction::EndOfRunAction(const G4Run* aRun)
134 {                                                 142 {
135   if (isMaster) fRun->EndOfRun();              << 143   std::ios::fmtflags mode = G4cout.flags();
136                                                << 144   G4cout.setf(std::ios::fixed,std::ios::floatfield);
137   // save histograms                           << 145   
                                                   >> 146   G4int NbofEvents = aRun->GetNumberOfEvent();
                                                   >> 147   if (NbofEvents == 0) return;
                                                   >> 148   G4double fNbofEvents = double(NbofEvents);
                                                   >> 149 
                                                   >> 150   //run conditions
                                                   >> 151   //     
                                                   >> 152   G4ParticleDefinition* particle = fKinematic->GetParticleGun()
                                                   >> 153                                           ->GetParticleDefinition();
                                                   >> 154   G4String partName = particle->GetParticleName();                             
                                                   >> 155   G4double energy = fKinematic->GetParticleGun()->GetParticleEnergy();
                                                   >> 156   
                                                   >> 157   G4int NbOfAbsor = fDetector->GetNbOfAbsor();
                                                   >> 158   
                                                   >> 159   G4cout << "\n ======================== run summary ======================\n";
                                                   >> 160   
                                                   >> 161   G4int prec = G4cout.precision(2);
                                                   >> 162   
                                                   >> 163   G4cout 
                                                   >> 164     << "\n The run consists of " << NbofEvents << " "<< partName << " of "
                                                   >> 165     << G4BestUnit(energy,"Energy") 
                                                   >> 166     << " through "  << NbOfAbsor << " absorbers: \n";
                                                   >> 167   for (G4int i=1; i<= NbOfAbsor; i++) {
                                                   >> 168      G4Material* material = fDetector->GetAbsorMaterial(i);
                                                   >> 169      G4double thickness = fDetector->GetAbsorThickness(i);
                                                   >> 170      G4double density = material->GetDensity();    
                                                   >> 171      G4cout << std::setw(20) << G4BestUnit(thickness,"Length") << " of "
                                                   >> 172             << material->GetName() << " (density: " 
                                                   >> 173             << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
                                                   >> 174   }         
                                                   >> 175   G4cout << "\n ============================================================\n";
                                                   >> 176   
                                                   >> 177   //compute total energy deposit
                                                   >> 178   //
                                                   >> 179   fEdeposit /= NbofEvents; fEdeposit2 /= NbofEvents;
                                                   >> 180   G4double rms = fEdeposit2 - fEdeposit*fEdeposit;        
                                                   >> 181   if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
                                                   >> 182 
                                                   >> 183   G4cout.precision(3);       
                                                   >> 184   G4cout 
                                                   >> 185     << "\n Total Energy deposited        = " << G4BestUnit(fEdeposit,"Energy")
                                                   >> 186     << " +- "                                << G4BestUnit( rms,"Energy")
                                                   >> 187     << G4endl;
                                                   >> 188               
                                                   >> 189   //compute track length of primary track
                                                   >> 190   //
                                                   >> 191   fTrackLen /= NbofEvents; fTrackLen2 /= NbofEvents;
                                                   >> 192   rms = fTrackLen2 - fTrackLen*fTrackLen;        
                                                   >> 193   if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
                                                   >> 194 
                                                   >> 195   G4cout.precision(3);       
                                                   >> 196   G4cout 
                                                   >> 197     << "\n Track length of primary track = " << G4BestUnit(fTrackLen,"Length")
                                                   >> 198     << " +- "                                << G4BestUnit( rms,"Length");
                                                   >> 199     
                                                   >> 200   //compare with csda range
                                                   >> 201   //
                                                   >> 202   if (NbOfAbsor == 1) {
                                                   >> 203     G4cout 
                                                   >> 204      << "\n Range from EmCalculator = " << G4BestUnit(fCsdaRange[1],"Length")
                                                   >> 205      << " (from full dE/dx)" << G4endl;
                                                   >> 206   }
                                                   >> 207                      
                                                   >> 208   //compute projected range of primary track
                                                   >> 209   //
                                                   >> 210   fProjRange /= NbofEvents; fProjRange2 /= NbofEvents;
                                                   >> 211   rms = fProjRange2 - fProjRange*fProjRange;        
                                                   >> 212   if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
                                                   >> 213    
                                                   >> 214   G4cout 
                                                   >> 215     << "\n Projected range               = " << G4BestUnit(fProjRange,"Length")
                                                   >> 216     << " +- "                                << G4BestUnit( rms,"Length")    
                                                   >> 217     << G4endl;
                                                   >> 218     
                                                   >> 219   //nb of steps and step size of primary track
                                                   >> 220   //
                                                   >> 221   G4double fNbSteps = fNbOfSteps/fNbofEvents, fNbSteps2 = fNbOfSteps2/fNbofEvents;
                                                   >> 222   rms = fNbSteps2 - fNbSteps*fNbSteps;       
                                                   >> 223   if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
                                                   >> 224 
                                                   >> 225   G4cout.precision(2);       
                                                   >> 226   G4cout << "\n Nb of steps of primary track  = " << fNbSteps << " +- " << rms;
                                                   >> 227     
                                                   >> 228   fStepSize /= NbofEvents; fStepSize2 /= NbofEvents;
                                                   >> 229   rms = fStepSize2 - fStepSize*fStepSize;        
                                                   >> 230   if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
                                                   >> 231 
                                                   >> 232   G4cout.precision(3);       
                                                   >> 233   G4cout 
                                                   >> 234     << "\t Step size= " << G4BestUnit(fStepSize,"Length")
                                                   >> 235     << " +- "           << G4BestUnit( rms,"Length")
                                                   >> 236     << G4endl;
                                                   >> 237     
                                                   >> 238   //transmission coefficients
                                                   >> 239   //
                                                   >> 240   G4double absorbed  = 100.*fStatus[0]/fNbofEvents;
                                                   >> 241   G4double transmit  = 100.*fStatus[1]/fNbofEvents;
                                                   >> 242   G4double reflected = 100.*fStatus[2]/fNbofEvents;  
                                                   >> 243 
                                                   >> 244   G4cout.precision(2);       
                                                   >> 245   G4cout 
                                                   >> 246     << "\n absorbed = "  << absorbed  << " %"
                                                   >> 247     << "   transmit = "  << transmit  << " %"
                                                   >> 248     << "   reflected = " << reflected << " %" << G4endl;
                                                   >> 249      
                                                   >> 250   // normalize histograms of longitudinal energy profile
                                                   >> 251   //
138   G4AnalysisManager* analysisManager = G4Analy    252   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
139   if (analysisManager->IsActive()) {           << 253   G4int ih = 1;
                                                   >> 254   G4double binWidth = analysisManager->GetH1Width(ih);
                                                   >> 255   G4double fac = (1./(NbofEvents*binWidth))*(mm/MeV);
                                                   >> 256   analysisManager->ScaleH1(ih,fac);
                                                   >> 257     
                                                   >> 258   ih = 8;
                                                   >> 259   binWidth = analysisManager->GetH1Width(ih);
                                                   >> 260   fac = (1./(NbofEvents*binWidth))*(g/(MeV*cm2));
                                                   >> 261   analysisManager->ScaleH1(ih,fac);
                                                   >> 262       
                                                   >> 263    // reset default formats
                                                   >> 264   G4cout.setf(mode,std::ios::floatfield);
                                                   >> 265   G4cout.precision(prec);
                                                   >> 266   
                                                   >> 267   // save histograms
                                                   >> 268   if ( analysisManager->IsActive() ) {  
140     analysisManager->Write();                     269     analysisManager->Write();
141     analysisManager->CloseFile();                 270     analysisManager->CloseFile();
142   }                                            << 271   }    
143                                                << 272  
144   // show Rndm status                             273   // show Rndm status
145   if (isMaster) G4Random::showEngineStatus();  << 274   CLHEP::HepRandom::showEngineStatus();
146 }                                                 275 }
147                                                   276 
148 //....oooOO0OOooo........oooOO0OOooo........oo    277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
149                                                   278