Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm3/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/TestEm3/src/RunAction.cc (Version 11.3.0) and /examples/extended/electromagnetic/TestEm3/src/RunAction.cc (Version 10.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/TestEm3/src/RunActio     26 /// \file electromagnetic/TestEm3/src/RunAction.cc
 27 /// \brief Implementation of the RunAction cla     27 /// \brief Implementation of the RunAction class
 28 //                                                 28 //
                                                   >>  29 // $Id: RunAction.cc 67268 2013-02-13 11:38:40Z ihrivnac $
 29 //                                                 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 
 35 #include "DetectorConstruction.hh"             << 
 36 #include "HistoManager.hh"                     << 
 37 #include "PrimaryGeneratorAction.hh"               36 #include "PrimaryGeneratorAction.hh"
 38 #include "Run.hh"                              << 
 39 #include "RunActionMessenger.hh"                   37 #include "RunActionMessenger.hh"
                                                   >>  38 #include "HistoManager.hh"
                                                   >>  39 #include "EmAcceptance.hh"
 40                                                    40 
                                                   >>  41 #include "G4Run.hh"
 41 #include "G4RunManager.hh"                         42 #include "G4RunManager.hh"
 42 #include "G4Timer.hh"                          <<  43 
                                                   >>  44 #include "G4ParticleTable.hh"
                                                   >>  45 #include "G4ParticleDefinition.hh"
                                                   >>  46 #include "G4Track.hh"
                                                   >>  47 #include "G4Gamma.hh"
                                                   >>  48 #include "G4Electron.hh"
                                                   >>  49 #include "G4Positron.hh"
                                                   >>  50 #include "G4ProductionCutsTable.hh"
                                                   >>  51 #include "G4LossTableManager.hh"
                                                   >>  52 
                                                   >>  53 #include "G4UnitsTable.hh"
                                                   >>  54 #include "G4SystemOfUnits.hh"
                                                   >>  55 
 43 #include "Randomize.hh"                            56 #include "Randomize.hh"
 44                                                    57 
 45 //....oooOO0OOooo........oooOO0OOooo........oo     58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 46                                                    59 
 47 RunAction::RunAction(DetectorConstruction* det     60 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
 48   : fDetector(det), fPrimary(prim)             <<  61 :G4UserRunAction(),fDetector(det), fPrimary(prim), fRunMessenger(0), 
                                                   >>  62  fHistoManager(0)
 49 {                                                  63 {
 50   fRunMessenger = new RunActionMessenger(this)     64   fRunMessenger = new RunActionMessenger(this);
 51   fHistoManager = new HistoManager();              65   fHistoManager = new HistoManager();
                                                   >>  66   fApplyLimit = false;
                                                   >>  67 
                                                   >>  68   fChargedStep = fNeutralStep = 0.0;
                                                   >>  69 
                                                   >>  70   for (G4int k=0; k<MaxAbsor; k++) { fEdeptrue[k] = fRmstrue[k] = 1.;
                                                   >>  71                                     fLimittrue[k] = DBL_MAX;
                                                   >>  72   }
 52 }                                                  73 }
 53                                                    74 
 54 //....oooOO0OOooo........oooOO0OOooo........oo     75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 55                                                    76 
 56 RunAction::~RunAction()                            77 RunAction::~RunAction()
 57 {                                                  78 {
 58   delete fRunMessenger;                            79   delete fRunMessenger;
 59 }                                                  80 }
 60                                                    81 
 61 //....oooOO0OOooo........oooOO0OOooo........oo     82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 62                                                    83 
 63 G4Run* RunAction::GenerateRun()                <<  84 void RunAction::BeginOfRunAction(const G4Run* aRun)
 64 {                                                  85 {
 65   fRun = new Run(fDetector);                   <<  86   G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
 66   return fRun;                                 << 
 67 }                                              << 
 68                                                    87 
 69 //....oooOO0OOooo........oooOO0OOooo........oo <<  88   // save Rndm status
                                                   >>  89   //
                                                   >>  90   ////G4RunManager::GetRunManager()->SetRandomNumberStore(true);
                                                   >>  91   CLHEP::HepRandom::showEngineStatus();
 70                                                    92 
 71 void RunAction::BeginOfRunAction(const G4Run*) <<  93   //initialize cumulative quantities
 72 {                                              <<  94   //
 73   // keep run condition                        <<  95   for (G4int k=0; k<MaxAbsor; k++) {
 74   if (fPrimary) {                              <<  96     fSumEAbs[k] = fSum2EAbs[k]  = fSumLAbs[k] = fSum2LAbs[k] = 0.;
 75     G4ParticleDefinition* particle = fPrimary- <<  97     fEnergyDeposit[k].clear();  
 76     G4double energy = fPrimary->GetParticleGun << 
 77     fRun->SetPrimary(particle, energy);        << 
 78   }                                                98   }
 79                                                    99 
 80   // histograms                                << 100   fChargedStep = fNeutralStep = 0.0;
                                                   >> 101 
                                                   >> 102   fN_gamma = 0;
                                                   >> 103   fN_elec  = 0;
                                                   >> 104   fN_pos   = 0;
                                                   >> 105 
                                                   >> 106   //initialize Eflow
                                                   >> 107   //
                                                   >> 108   G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2;
                                                   >> 109   fEnergyFlow.resize(nbPlanes);
                                                   >> 110   fLateralEleak.resize(nbPlanes);
                                                   >> 111   for (G4int k=0; k<nbPlanes; k++) {fEnergyFlow[k] = fLateralEleak[k] = 0.; }
                                                   >> 112   
                                                   >> 113   //histograms
 81   //                                              114   //
 82   G4AnalysisManager* analysis = G4AnalysisMana    115   G4AnalysisManager* analysis = G4AnalysisManager::Instance();
 83   if (analysis->IsActive()) analysis->OpenFile    116   if (analysis->IsActive()) analysis->OpenFile();
                                                   >> 117    
                                                   >> 118   //example of print dEdx tables
                                                   >> 119   //
                                                   >> 120   ////PrintDedxTables();
                                                   >> 121 }
 84                                                   122 
 85   // save Rndm status and open the timer       << 123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 86                                                   124 
 87   if (isMaster) {                              << 125 void RunAction::FillPerEvent(G4int kAbs, G4double EAbs, G4double LAbs)
 88     //    G4Random::showEngineStatus();        << 126 {
 89     fTimer = new G4Timer();                    << 127   //accumulate statistic with restriction
 90     fTimer->Start();                           << 128   //
 91   }                                            << 129   if(fApplyLimit) fEnergyDeposit[kAbs].push_back(EAbs);
                                                   >> 130   fSumEAbs[kAbs]  += EAbs;  fSum2EAbs[kAbs]  += EAbs*EAbs;
                                                   >> 131   fSumLAbs[kAbs]  += LAbs;  fSum2LAbs[kAbs]  += LAbs*LAbs;
 92 }                                                 132 }
 93                                                   133 
 94 //....oooOO0OOooo........oooOO0OOooo........oo    134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 95                                                   135 
 96 void RunAction::EndOfRunAction(const G4Run*)   << 136 
                                                   >> 137 void RunAction::EndOfRunAction(const G4Run* aRun)
 97 {                                                 138 {
 98   // compute and print statistic               << 139   G4int nEvt = aRun->GetNumberOfEvent();
 99   if (isMaster) {                              << 140   G4double  norm = G4double(nEvt);
100     fTimer->Stop();                            << 141   if(norm > 0) norm = 1./norm;
101     if (!((G4RunManager::GetRunManager()->GetR << 142   G4double qnorm = std::sqrt(norm);
102       G4cout << "\n"                           << 143 
103              << "Total number of events:  " << << 144   fChargedStep *= norm;
104       G4cout << "Master thread time:  " << *fT << 145   fNeutralStep *= norm;
                                                   >> 146 
                                                   >> 147   //compute and print statistic
                                                   >> 148   //
                                                   >> 149   G4double beamEnergy = fPrimary->GetParticleGun()->GetParticleEnergy();
                                                   >> 150   G4double sqbeam = std::sqrt(beamEnergy/GeV);
                                                   >> 151 
                                                   >> 152   G4double MeanEAbs,MeanEAbs2,rmsEAbs,resolution,rmsres;
                                                   >> 153   G4double MeanLAbs,MeanLAbs2,rmsLAbs;
                                                   >> 154 
                                                   >> 155   std::ios::fmtflags mode = G4cout.flags();
                                                   >> 156   G4int  prec = G4cout.precision(2);
                                                   >> 157   G4cout << "\n------------------------------------------------------------\n";
                                                   >> 158   G4cout << std::setw(14) << "material"
                                                   >> 159          << std::setw(17) << "Edep       RMS"
                                                   >> 160          << std::setw(33) << "sqrt(E0(GeV))*rmsE/Emean"
                                                   >> 161          << std::setw(23) << "total tracklen \n \n";
                                                   >> 162 
                                                   >> 163   for (G4int k=1; k<=fDetector->GetNbOfAbsor(); k++)
                                                   >> 164     {
                                                   >> 165       MeanEAbs  = fSumEAbs[k]*norm;
                                                   >> 166       MeanEAbs2 = fSum2EAbs[k]*norm;
                                                   >> 167       rmsEAbs  = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
                                                   >> 168       //G4cout << "k= " << k << "  RMS= " <<  rmsEAbs 
                                                   >> 169       //     << "  fApplyLimit: " << fApplyLimit << G4endl;
                                                   >> 170       if(fApplyLimit) {
                                                   >> 171         G4int    nn    = 0;
                                                   >> 172         G4double sume  = 0.0;
                                                   >> 173         G4double sume2 = 0.0;
                                                   >> 174         // compute trancated means  
                                                   >> 175         G4double lim   = rmsEAbs * 2.5;
                                                   >> 176         for(G4int i=0; i<nEvt; i++) {
                                                   >> 177           G4double e = (fEnergyDeposit[k])[i];
                                                   >> 178           if(std::abs(e - MeanEAbs) < lim) {
                                                   >> 179             sume  += e;
                                                   >> 180             sume2 += e*e;
                                                   >> 181             nn++;
                                                   >> 182           }
                                                   >> 183         }
                                                   >> 184         G4double norm1 = G4double(nn);
                                                   >> 185         if(norm1 > 0.0) norm1 = 1.0/norm1;
                                                   >> 186         MeanEAbs  = sume*norm1;
                                                   >> 187         MeanEAbs2 = sume2*norm1;
                                                   >> 188         rmsEAbs  = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
                                                   >> 189       }
                                                   >> 190 
                                                   >> 191       resolution= 100.*sqbeam*rmsEAbs/MeanEAbs;
                                                   >> 192       rmsres    = resolution*qnorm;
                                                   >> 193 
                                                   >> 194       // Save mean and RMS
                                                   >> 195       fSumEAbs[k] = MeanEAbs;
                                                   >> 196       fSum2EAbs[k] = rmsEAbs;
                                                   >> 197 
                                                   >> 198       MeanLAbs  = fSumLAbs[k]*norm;
                                                   >> 199       MeanLAbs2 = fSum2LAbs[k]*norm;
                                                   >> 200       rmsLAbs  = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs*MeanLAbs));
                                                   >> 201 
                                                   >> 202       //print
                                                   >> 203       //
                                                   >> 204       G4cout
                                                   >> 205        << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName() << ": "
                                                   >> 206        << std::setprecision(5)
                                                   >> 207        << std::setw(6) << G4BestUnit(MeanEAbs,"Energy") << " :  "
                                                   >> 208        << std::setprecision(4)
                                                   >> 209        << std::setw(5) << G4BestUnit( rmsEAbs,"Energy")  
                                                   >> 210        << std::setw(10) << resolution  << " +- " 
                                                   >> 211        << std::setw(5) << rmsres << " %"
                                                   >> 212        << std::setprecision(3)
                                                   >> 213        << std::setw(10) << G4BestUnit(MeanLAbs,"Length")  << " +- "
                                                   >> 214        << std::setw(4) << G4BestUnit( rmsLAbs,"Length")
                                                   >> 215        << G4endl;
105     }                                             216     }
106     delete fTimer;                             << 217   G4cout << "\n------------------------------------------------------------\n";
107     fRun->EndOfRun();                          << 218 
108   }                                            << 219   G4cout << " Beam particle " 
109   // save histograms                           << 220          << fPrimary->GetParticleGun()->
                                                   >> 221     GetParticleDefinition()->GetParticleName()
                                                   >> 222          << "  E = " << G4BestUnit(beamEnergy,"Energy") << G4endl;
                                                   >> 223   G4cout << " Mean number of gamma          " << (G4double)fN_gamma*norm << G4endl;
                                                   >> 224   G4cout << " Mean number of e-             " << (G4double)fN_elec*norm << G4endl;
                                                   >> 225   G4cout << " Mean number of e+             " << (G4double)fN_pos*norm << G4endl;
                                                   >> 226   G4cout << std::setprecision(6)
                                                   >> 227    << " Mean number of charged steps  " << fChargedStep << G4endl;
                                                   >> 228   G4cout << " Mean number of neutral steps  " << fNeutralStep << G4endl;
                                                   >> 229   G4cout << "------------------------------------------------------------\n";
                                                   >> 230   
                                                   >> 231   //Energy flow
                                                   >> 232   //
110   G4AnalysisManager* analysis = G4AnalysisMana    233   G4AnalysisManager* analysis = G4AnalysisManager::Instance();
111   if (analysis->IsActive()) {                  << 234   G4int Idmax = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor());
                                                   >> 235   for (G4int Id=1; Id<=Idmax+1; Id++) {
                                                   >> 236     analysis->FillH1(2*MaxAbsor+1, (G4double)Id, fEnergyFlow[Id]);
                                                   >> 237     analysis->FillH1(2*MaxAbsor+2, (G4double)Id, fLateralEleak[Id]);
                                                   >> 238   }
                                                   >> 239   
                                                   >> 240   //Energy deposit from energy flow balance
                                                   >> 241   //
                                                   >> 242   G4double EdepTot[MaxAbsor];
                                                   >> 243   for (G4int k=0; k<MaxAbsor; k++) EdepTot[k] = 0.;
                                                   >> 244   
                                                   >> 245   G4int nbOfAbsor = fDetector->GetNbOfAbsor();
                                                   >> 246   for (G4int Id=1; Id<=Idmax; Id++) {
                                                   >> 247     G4int iAbsor = Id%nbOfAbsor; if (iAbsor==0) iAbsor = nbOfAbsor;
                                                   >> 248     EdepTot [iAbsor] += (fEnergyFlow[Id] - fEnergyFlow[Id+1] - fLateralEleak[Id]);
                                                   >> 249   }
                                                   >> 250   
                                                   >> 251   G4cout << std::setprecision(3)
                                                   >> 252    << "\n Energy deposition from Energy flow balance : \n"
                                                   >> 253          << std::setw(10) << "  material \t Total Edep \n \n";
                                                   >> 254   G4cout.precision(6);
                                                   >> 255   
                                                   >> 256   for (G4int k=1; k<=nbOfAbsor; k++) {
                                                   >> 257     EdepTot [k] *= norm;
                                                   >> 258     G4cout << std::setw(10) << fDetector->GetAbsorMaterial(k)->GetName() << ":"
                                                   >> 259            << "\t " << G4BestUnit(EdepTot [k],"Energy") << "\n";
                                                   >> 260   }
                                                   >> 261   
                                                   >> 262   G4cout << "\n------------------------------------------------------------\n" 
                                                   >> 263          << G4endl;
                                                   >> 264     
                                                   >> 265   G4cout.setf(mode,std::ios::floatfield);
                                                   >> 266   G4cout.precision(prec);
                                                   >> 267 
                                                   >> 268   // Acceptance
                                                   >> 269   EmAcceptance acc;
                                                   >> 270   G4bool isStarted = false;
                                                   >> 271   for (G4int j=1; j<=fDetector->GetNbOfAbsor(); j++) {
                                                   >> 272     if (fLimittrue[j] < DBL_MAX) {
                                                   >> 273       if (!isStarted) {
                                                   >> 274         acc.BeginOfAcceptance("Sampling Calorimeter",nEvt);
                                                   >> 275         isStarted = true;
                                                   >> 276       }
                                                   >> 277       MeanEAbs = fSumEAbs[j];
                                                   >> 278       rmsEAbs  = fSum2EAbs[j];
                                                   >> 279       G4String mat = fDetector->GetAbsorMaterial(j)->GetName();
                                                   >> 280       acc.EmAcceptanceGauss("Edep"+mat, nEvt, MeanEAbs,
                                                   >> 281                              fEdeptrue[j], fRmstrue[j], fLimittrue[j]);
                                                   >> 282       acc.EmAcceptanceGauss("Erms"+mat, nEvt, rmsEAbs,
                                                   >> 283                              fRmstrue[j], fRmstrue[j], 2.0*fLimittrue[j]);
                                                   >> 284     }
                                                   >> 285   }
                                                   >> 286   if(isStarted) acc.EndOfAcceptance();
                                                   >> 287 
                                                   >> 288   //normalize histograms
                                                   >> 289   //
                                                   >> 290   for (G4int ih = MaxAbsor+1; ih < MaxHisto; ih++) {
                                                   >> 291     analysis->ScaleH1(ih,norm/MeV);
                                                   >> 292   }
                                                   >> 293   
                                                   >> 294   //save histograms
                                                   >> 295   if (analysis->IsActive()) {   
112     analysis->Write();                            296     analysis->Write();
113     analysis->CloseFile();                        297     analysis->CloseFile();
114   }                                            << 298   }    
115                                                   299 
116   // show Rndm status                             300   // show Rndm status
117   //  if (isMaster)  G4Random::showEngineStatu << 301   CLHEP::HepRandom::showEngineStatus();
118 }                                                 302 }
119                                                   303 
120 //....oooOO0OOooo........oooOO0OOooo........oo    304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
121                                                   305 
122 void RunAction::SetEdepAndRMS(G4int i, G4doubl << 306 void RunAction::PrintDedxTables()
                                                   >> 307 {
                                                   >> 308   //Print dE/dx tables with binning identical to the Geant3 JMATE bank.
                                                   >> 309   //The printout is readable as Geant3 ffread data cards (by the program g4mat).
                                                   >> 310   //
                                                   >> 311   const G4double tkmin=10*keV, tkmax=10*TeV;
                                                   >> 312   const G4int nbin=90;
                                                   >> 313   G4double tk[nbin];
                                                   >> 314 
                                                   >> 315   const G4int ncolumn = 5;
                                                   >> 316 
                                                   >> 317   //compute the kinetic energies
                                                   >> 318   //
                                                   >> 319   const G4double dp = std::log10(tkmax/tkmin)/nbin;
                                                   >> 320   const G4double dt = std::pow(10.,dp);
                                                   >> 321   tk[0] = tkmin;
                                                   >> 322   for (G4int i=1; i<nbin; ++i) tk[i] = tk[i-1]*dt;
                                                   >> 323 
                                                   >> 324   //print the kinetic energies
                                                   >> 325   //
                                                   >> 326   std::ios::fmtflags mode = G4cout.flags();
                                                   >> 327   G4cout.setf(std::ios::fixed,std::ios::floatfield);
                                                   >> 328   G4int  prec = G4cout.precision(3);
                                                   >> 329 
                                                   >> 330   G4cout << "\n kinetic energies \n ";
                                                   >> 331   for (G4int j=0; j<nbin; ++j) {
                                                   >> 332     G4cout << G4BestUnit(tk[j],"Energy") << "\t";
                                                   >> 333     if ((j+1)%ncolumn == 0) G4cout << "\n ";
                                                   >> 334   }
                                                   >> 335   G4cout << G4endl;
                                                   >> 336 
                                                   >> 337   //print the dE/dx tables
                                                   >> 338   //
                                                   >> 339   G4cout.setf(std::ios::scientific,std::ios::floatfield);
                                                   >> 340 
                                                   >> 341   G4ParticleDefinition*
                                                   >> 342   part = fPrimary->GetParticleGun()->GetParticleDefinition();
                                                   >> 343   
                                                   >> 344   G4ProductionCutsTable* theCoupleTable =
                                                   >> 345         G4ProductionCutsTable::GetProductionCutsTable();
                                                   >> 346   size_t numOfCouples = theCoupleTable->GetTableSize();
                                                   >> 347   const G4MaterialCutsCouple* couple = 0;
                                                   >> 348 
                                                   >> 349   for (G4int iab=1;iab <= fDetector->GetNbOfAbsor(); iab++)
                                                   >> 350      {
                                                   >> 351       G4Material* mat = fDetector->GetAbsorMaterial(iab);
                                                   >> 352       G4int index = 0;
                                                   >> 353       for (size_t i=0; i<numOfCouples; i++) {
                                                   >> 354          couple = theCoupleTable->GetMaterialCutsCouple(i);
                                                   >> 355          if (couple->GetMaterial() == mat) {index = i; break;}
                                                   >> 356       }
                                                   >> 357       G4cout << "\nLIST";
                                                   >> 358       G4cout << "\nC \nC  dE/dx (MeV/cm) for " << part->GetParticleName()
                                                   >> 359              << " in " << mat ->GetName() << "\nC";
                                                   >> 360       G4cout << "\nKINE   (" << part->GetParticleName() << ")";
                                                   >> 361       G4cout << "\nMATE   (" << mat ->GetName() << ")";
                                                   >> 362       G4cout.precision(2);
                                                   >> 363       G4cout << "\nERAN  " << tkmin/GeV << " (ekmin)\t"
                                                   >> 364                            << tkmax/GeV << " (ekmax)\t"
                                                   >> 365                            << nbin      << " (nekbin)";
                                                   >> 366       G4double cutgam =
                                                   >> 367          (*(theCoupleTable->GetEnergyCutsVector(idxG4GammaCut)))[index];
                                                   >> 368       if (cutgam < tkmin) cutgam = tkmin;
                                                   >> 369       if (cutgam > tkmax) cutgam = tkmax;
                                                   >> 370       G4double cutele =
                                                   >> 371          (*(theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut)))[index];
                                                   >> 372       if (cutele < tkmin) cutele = tkmin;
                                                   >> 373       if (cutele > tkmax) cutele = tkmax;
                                                   >> 374       G4cout << "\nCUTS  " << cutgam/GeV << " (cutgam)\t"
                                                   >> 375                            << cutele/GeV << " (cutele)";
                                                   >> 376 
                                                   >> 377       G4cout.precision(6);
                                                   >> 378       G4cout << "\nG4VAL \n ";
                                                   >> 379       for (G4int l=0;l<nbin; ++l)
                                                   >> 380          {
                                                   >> 381            G4double dedx = G4LossTableManager::Instance()
                                                   >> 382                                                ->GetDEDX(part,tk[l],couple);
                                                   >> 383            G4cout << dedx/(MeV/cm) << "\t";
                                                   >> 384            if ((l+1)%ncolumn == 0) G4cout << "\n ";
                                                   >> 385          }
                                                   >> 386       G4cout << G4endl;
                                                   >> 387      }
                                                   >> 388 
                                                   >> 389   G4cout.precision(prec);
                                                   >> 390   G4cout.setf(mode,std::ios::floatfield);
                                                   >> 391 }
                                                   >> 392 
                                                   >> 393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 394 
                                                   >> 395 void RunAction::AddSecondaryTrack(const G4Track* track)
123 {                                                 396 {
124   if (fRun) fRun->SetEdepAndRMS(i, edep, rms,  << 397   const G4ParticleDefinition* d = track->GetDefinition();
                                                   >> 398   if(d == G4Gamma::Gamma()) { ++fN_gamma; }
                                                   >> 399   else if (d == G4Electron::Electron()) { ++fN_elec; }
                                                   >> 400   else if (d == G4Positron::Positron()) { ++fN_pos; }
125 }                                                 401 }
126                                                   402 
127 //....oooOO0OOooo........oooOO0OOooo........oo    403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128                                                   404 
129 void RunAction::SetApplyLimit(G4bool val)      << 405 void RunAction::SetEdepAndRMS(G4int i, G4double edep, G4double rms, G4double lim)
130 {                                                 406 {
131   if (fRun) fRun->SetApplyLimit(val);          << 407   if (i>=0 && i<MaxAbsor) {
                                                   >> 408     fEdeptrue [i] = edep;
                                                   >> 409     fRmstrue  [i] = rms;
                                                   >> 410     fLimittrue[i] = lim;
                                                   >> 411   }
132 }                                                 412 }
133                                                   413 
134 //....oooOO0OOooo........oooOO0OOooo........oo    414 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135                                                   415