Geant4 Cross Reference

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


  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/TestEm17/src/RunActi     26 /// \file electromagnetic/TestEm17/src/RunAction.cc
 27 /// \brief Implementation of the RunAction cla     27 /// \brief Implementation of the RunAction class
 28 //                                                 28 //
 29 //                                             <<  29 // 
 30 //....oooOO0OOooo........oooOO0OOooo........oo     30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oo     31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32                                                    32 
 33 #include "RunAction.hh"                            33 #include "RunAction.hh"
 34                                                    34 
 35 #include "DetectorConstruction.hh"                 35 #include "DetectorConstruction.hh"
                                                   >>  36 #include "PrimaryGeneratorAction.hh"
 36 #include "HistoManager.hh"                         37 #include "HistoManager.hh"
 37 #include "MuCrossSections.hh"                      38 #include "MuCrossSections.hh"
 38 #include "PrimaryGeneratorAction.hh"           << 
 39                                                << 
 40 #include "G4EmCalculator.hh"                   << 
 41 #include "G4PhysicalConstants.hh"              << 
 42 #include "G4ProductionCutsTable.hh"                39 #include "G4ProductionCutsTable.hh"
                                                   >>  40 
 43 #include "G4Run.hh"                                41 #include "G4Run.hh"
 44 #include "G4RunManager.hh"                         42 #include "G4RunManager.hh"
 45 #include "G4SystemOfUnits.hh"                  << 
 46 #include "G4UnitsTable.hh"                         43 #include "G4UnitsTable.hh"
                                                   >>  44 #include "G4EmCalculator.hh"
                                                   >>  45 
                                                   >>  46 #include "G4PhysicalConstants.hh"
                                                   >>  47 #include "G4SystemOfUnits.hh"
 47 #include "Randomize.hh"                            48 #include "Randomize.hh"
 48                                                    49 
 49 //....oooOO0OOooo........oooOO0OOooo........oo     50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 50                                                    51 
 51 RunAction::RunAction(DetectorConstruction* det <<  52 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim,
 52   : G4UserRunAction(), fDetector(det), fPrimar <<  53                      HistoManager* HistM)
                                                   >>  54   : G4UserRunAction(),
                                                   >>  55     fDetector(det), fPrimary(prim), fProcCounter(0), fHistoManager(HistM)
 53 {                                                  56 {
 54   fMucs = new MuCrossSections();                   57   fMucs = new MuCrossSections();
 55 }                                                  58 }
 56                                                    59 
 57 //....oooOO0OOooo........oooOO0OOooo........oo     60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 58                                                    61 
 59 RunAction::~RunAction()                            62 RunAction::~RunAction()
 60 {                                                  63 {
 61   delete fMucs;                                    64   delete fMucs;
 62 }                                                  65 }
 63                                                    66 
 64 //....oooOO0OOooo........oooOO0OOooo........oo     67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 65                                                    68 
 66 void RunAction::BeginOfRunAction(const G4Run*      69 void RunAction::BeginOfRunAction(const G4Run* aRun)
 67 {                                              <<  70 {  
 68   G4cout << "### Run " << aRun->GetRunID() <<      71   G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
 69                                                <<  72   
 70   // save Rndm status                              73   // save Rndm status
 71   CLHEP::HepRandom::showEngineStatus();            74   CLHEP::HepRandom::showEngineStatus();
 72                                                    75 
 73   fProcCounter = new ProcessesCount();             76   fProcCounter = new ProcessesCount();
 74   fHistoManager->Book();                           77   fHistoManager->Book();
 75 }                                                  78 }
 76                                                    79 
 77 //....oooOO0OOooo........oooOO0OOooo........oo     80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 78                                                    81 
 79 void RunAction::CountProcesses(const G4String&     82 void RunAction::CountProcesses(const G4String& procName)
 80 {                                                  83 {
 81   // does the process  already encounted ?     <<  84    //does the process  already encounted ?
 82   size_t n = fProcCounter->size();             <<  85    size_t n = fProcCounter->size();
 83   for (size_t i = 0; i < n; ++i) {             <<  86    for(size_t i = 0; i<n; ++i) {
 84     if ((*fProcCounter)[i]->GetName() == procN <<  87      if((*fProcCounter)[i]->GetName()==procName) {
 85       (*fProcCounter)[i]->Count();             <<  88        (*fProcCounter)[i]->Count();
 86       return;                                  <<  89        return;
 87     }                                          <<  90      }
 88   }                                            <<  91    }
 89   OneProcessCount* count = new OneProcessCount <<  92    OneProcessCount* count = new OneProcessCount(procName);
 90   count->Count();                              <<  93    count->Count();
 91   fProcCounter->push_back(count);              <<  94    fProcCounter->push_back(count);
 92 }                                                  95 }
 93                                                    96 
 94 //....oooOO0OOooo........oooOO0OOooo........oo     97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 95                                                    98 
 96 void RunAction::EndOfRunAction(const G4Run* aR     99 void RunAction::EndOfRunAction(const G4Run* aRun)
 97 {                                                 100 {
 98   G4int NbOfEvents = aRun->GetNumberOfEvent();    101   G4int NbOfEvents = aRun->GetNumberOfEvent();
 99   if (NbOfEvents == 0) return;                    102   if (NbOfEvents == 0) return;
100                                                << 103   
101   //  std::ios::fmtflags mode = G4cout.flags()    104   //  std::ios::fmtflags mode = G4cout.flags();
102   G4int prec = G4cout.precision(2);            << 105   G4int  prec = G4cout.precision(2);
103                                                << 106     
104   const G4Material* material = fDetector->GetM    107   const G4Material* material = fDetector->GetMaterial();
105   G4double length = fDetector->GetSize();      << 108   G4double length  = fDetector->GetSize();
106   G4double density = material->GetDensity();      109   G4double density = material->GetDensity();
107                                                << 110    
108   G4String particle = fPrimary->GetParticleGun << 111   G4String particle = fPrimary->GetParticleGun()->GetParticleDefinition()
                                                   >> 112                       ->GetParticleName();    
109   G4double energy = fPrimary->GetParticleGun()    113   G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
110                                                << 114   
111   G4cout << "\n The run consists of " << NbOfE << 115   G4cout << "\n The run consists of " << NbOfEvents << " "<< particle << " of "
112          << G4BestUnit(energy, "Energy") << "  << 116          << G4BestUnit(energy,"Energy") << " through " 
113          << material->GetName() << " (density: << 117          << G4BestUnit(length,"Length") << " of "
114          << G4endl;                            << 118          << material->GetName() << " (density: " 
115                                                << 119          << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
116   // total number of process calls             << 120   
                                                   >> 121   //total number of process calls
117   G4double countTot = 0.;                         122   G4double countTot = 0.;
118   G4cout << "\n Number of process calls --->";    123   G4cout << "\n Number of process calls --->";
119   for (size_t i = 0; i < fProcCounter->size(); << 124   for (size_t i=0; i< fProcCounter->size();++i) {
120     G4String procName = (*fProcCounter)[i]->Ge << 125      G4String procName = (*fProcCounter)[i]->GetName();
121     if (procName != "Transportation") {        << 126      if (procName != "Transportation") {
122       G4int count = (*fProcCounter)[i]->GetCou << 127        G4int count = (*fProcCounter)[i]->GetCounter(); 
123       G4cout << "\t" << procName << " : " << c << 128        G4cout << "\t" << procName << " : " << count;
124       countTot += count;                       << 129        countTot += count;
125     }                                          << 130      }
126   }                                               131   }
127                                                << 132   
128   // compute totalCrossSection, meanFreePath a << 133   //compute totalCrossSection, meanFreePath and massicCrossSection
129   //                                              134   //
130   G4double totalCrossSection = countTot / (NbO << 135   G4double totalCrossSection = countTot/(NbOfEvents*length);
131   G4double MeanFreePath = 1. / totalCrossSecti << 136   G4double MeanFreePath      = 1./totalCrossSection;        
132   G4double massCrossSection = totalCrossSectio << 137   G4double massCrossSection  = totalCrossSection/density;     
133                                                << 138    
134   G4cout.precision(5);                            139   G4cout.precision(5);
135   G4cout << "\n Simulation: "                     140   G4cout << "\n Simulation: "
136          << "total CrossSection = " << totalCr << 141          <<    "total CrossSection = " << totalCrossSection*cm << " /cm"
137          << "\t MeanFreePath = " << G4BestUnit << 142          << "\t MeanFreePath = "       << G4BestUnit(MeanFreePath,"Length")
138          << "\t massicCrossSection = " << mass << 143          << "\t massicCrossSection = " << massCrossSection*g/cm2 << " cm2/g"
139                                                << 144          << G4endl;
140   // compute theoretical predictions           << 145   
                                                   >> 146   //compute theoretical predictions
141   //                                              147   //
142   if (particle == "mu+" || particle == "mu-")  << 148   if(particle == "mu+" || particle == "mu-") { 
143     totalCrossSection = 0.;                       149     totalCrossSection = 0.;
144     for (size_t i = 0; i < fProcCounter->size( << 150     for (size_t i=0; i< fProcCounter->size();++i) {
145       G4String procName = (*fProcCounter)[i]->    151       G4String procName = (*fProcCounter)[i]->GetName();
146       if (procName != "Transportation") {         152       if (procName != "Transportation") {
147         totalCrossSection += ComputeTheory(pro    153         totalCrossSection += ComputeTheory(procName, NbOfEvents);
148         FillCrossSectionHisto(procName, NbOfEv    154         FillCrossSectionHisto(procName, NbOfEvents);
149       }                                           155       }
150     }                                             156     }
151                                                << 157   
152     MeanFreePath = 1. / totalCrossSection;     << 158     MeanFreePath     = 1./totalCrossSection;
153     massCrossSection = totalCrossSection / den << 159     massCrossSection = totalCrossSection/density;
154                                                << 160   
155     G4cout << " Theory:     "                     161     G4cout << " Theory:     "
156            << "total CrossSection = " << total << 162            <<    "total CrossSection = " << totalCrossSection*cm << " /cm"
157            << "\t MeanFreePath = " << G4BestUn << 163            << "\t MeanFreePath = "       << G4BestUnit(MeanFreePath,"Length")
158            << "\t massicCrossSection = " << ma << 164            << "\t massicCrossSection = " << massCrossSection*g/cm2 << " cm2/g"
                                                   >> 165            << G4endl;
159   }                                               166   }
160                                                << 167                                                                             
161   //  G4cout.setf(mode,std::ios::floatfield);     168   //  G4cout.setf(mode,std::ios::floatfield);
162   G4cout.precision(prec);                      << 169   G4cout.precision(prec);         
163                                                   170 
164   // delete and remove all contents in fProcCo << 171   // delete and remove all contents in fProcCounter 
165   size_t n = fProcCounter->size();                172   size_t n = fProcCounter->size();
166   for (size_t i = 0; i < n; ++i) {             << 173   for(size_t i = 0; i<n; ++i) { delete (*fProcCounter)[i]; }
167     delete (*fProcCounter)[i];                 << 
168   }                                            << 
169   delete fProcCounter;                            174   delete fProcCounter;
170                                                << 175   
171   fHistoManager->Save();                          176   fHistoManager->Save();
172                                                << 177   
173   // show Rndm status                             178   // show Rndm status
174   // CLHEP::HepRandom::showEngineStatus();     << 179   //CLHEP::HepRandom::showEngineStatus();
175 }                                                 180 }
176                                                   181 
177 //....oooOO0OOooo........oooOO0OOooo........oo    182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
178                                                   183 
179 G4double RunAction::ComputeTheory(const G4Stri << 184 G4double RunAction::ComputeTheory(const G4String& process, G4int NbOfMu)    
180 {                                              << 185 {   
181   const G4Material* material = fDetector->GetM    186   const G4Material* material = fDetector->GetMaterial();
182   G4double ekin = fPrimary->GetParticleGun()->    187   G4double ekin = fPrimary->GetParticleGun()->GetParticleEnergy();
183   G4double particleMass = fPrimary->GetParticl    188   G4double particleMass = fPrimary->GetParticleGun()->GetParticleDefinition()->GetPDGMass();
184                                                   189 
185   G4int id = 0;                                << 190   G4int id = 0; G4double cut = 1.e-10*ekin;
186   G4double cut = 1.e-10 * ekin;                << 191   if (process == "muIoni")          {id = 11; cut =  GetEnergyCut(material,1);}
187   if (process == "muIoni") {                   << 192   else if (process == "muPairProd") {id = 12; cut = 2*(GetEnergyCut(material,1)
188     id = 11;                                   << 193                                                       + electron_mass_c2); }
189     cut = GetEnergyCut(material, 1);           << 194   else if (process == "muBrems")    {id = 13; cut =  GetEnergyCut(material,0);}
190   }                                            << 195   else if (process == "muonNuclear"){id = 14; cut = 100*MeV;}
191   else if (process == "muPairProd") {          << 196   else if (process == "muToMuonPairProd"){id = 18; cut = 2*particleMass;}
192     id = 12;                                   << 197   if (id == 0) { return 0.; }
193     cut = 2 * (GetEnergyCut(material, 1) + ele << 198   
194   }                                            << 
195   else if (process == "muBrems") {             << 
196     id = 13;                                   << 
197     cut = GetEnergyCut(material, 0);           << 
198   }                                            << 
199   else if (process == "muonNuclear") {         << 
200     id = 14;                                   << 
201     cut = 100 * MeV;                           << 
202   }                                            << 
203   else if (process == "muToMuonPairProd") {    << 
204     id = 18;                                   << 
205     cut = 2 * particleMass;                    << 
206   }                                            << 
207   if (id == 0) {                               << 
208     return 0.;                                 << 
209   }                                            << 
210                                                << 
211   G4int nbOfBins = 100;                           199   G4int nbOfBins = 100;
212   // G4double binMin = -10.;                   << 200   //G4double binMin = -10.;
213   G4double binMin = std::log10(cut / ekin);    << 201   G4double binMin = std::log10(cut/ekin);
214   G4double binMax = 0.;                           202   G4double binMax = 0.;
215   G4double binWidth = (binMax - binMin) / G4do << 203   G4double binWidth = (binMax-binMin)/G4double(nbOfBins);
216                                                   204 
217   // create histo for theoretical crossSection << 205   //create histo for theoretical crossSections, with same bining as simulation
218   //                                              206   //
219   G4AnalysisManager* analysisManager = G4Analy    207   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
220                                                << 208     
221   G4H1* histoTh = 0;                              209   G4H1* histoTh = 0;
222   if (fHistoManager->HistoExist(id)) {            210   if (fHistoManager->HistoExist(id)) {
223     histoTh = analysisManager->GetH1(fHistoMan << 211     histoTh  = analysisManager->GetH1(fHistoManager->GetHistoID(id));  
224     nbOfBins = fHistoManager->GetNbins(id);       212     nbOfBins = fHistoManager->GetNbins(id);
225     binMin = fHistoManager->GetVmin(id);       << 213     binMin   = fHistoManager->GetVmin (id);
226     binMax = fHistoManager->GetVmax(id);       << 214     binMax   = fHistoManager->GetVmax (id);
227     binWidth = fHistoManager->GetBinWidth(id); << 215     binWidth = fHistoManager->GetBinWidth(id);    
228   }                                            << 216   }
229                                                << 217   
230   // compute and plot differential crossSectio << 218   //compute and plot differential crossSection, as function of energy transfert.
231   // compute and return integrated crossSectio << 219   //compute and return integrated crossSection for a given process.
232   //(note: to compare with simulation, the int    220   //(note: to compare with simulation, the integrated crossSection is function
233   //        of the energy cut.)                << 221   //       of the energy cut.) 
234   //                                           << 222   // 
235   G4double lgeps, etransf, sigmaE, dsigma;        223   G4double lgeps, etransf, sigmaE, dsigma;
236   G4double sigmaTot = 0.;                         224   G4double sigmaTot = 0.;
237   const G4double ln10 = std::log(10.);         << 225   const G4double ln10 = std::log(10.);  
238   G4double length = fDetector->GetSize();         226   G4double length = fDetector->GetSize();
239                                                   227 
240   // G4cout << "MU: " << process << " E= " <<  << 228   //G4cout << "MU: " << process << " E= " << ekin 
241   //        <<"  binMin= " << binMin << " binW << 229   //       <<"  binMin= " << binMin << " binW= " << binWidth << G4endl;
242                                                   230 
243   for (G4int ibin = 0; ibin < nbOfBins; ibin++ << 231   for (G4int ibin=0; ibin<nbOfBins; ibin++) {
244     lgeps = binMin + (ibin + 0.5) * binWidth;  << 232     lgeps = binMin + (ibin+0.5)*binWidth;
245     etransf = ekin * std::pow(10., lgeps);     << 233     etransf = ekin*std::pow(10.,lgeps);
246     sigmaE = fMucs->CR_Macroscopic(process, ma << 234     sigmaE = fMucs->CR_Macroscopic(process,material,ekin,etransf);
247     dsigma = sigmaE * etransf * binWidth * ln1 << 235     dsigma = sigmaE*etransf*binWidth*ln10;
248     if (etransf > cut) sigmaTot += dsigma;     << 236     if (etransf > cut) sigmaTot += dsigma;    
249     if (histoTh) {                                237     if (histoTh) {
250       G4double NbProcess = NbOfMu * length * d << 238       G4double NbProcess = NbOfMu*length*dsigma;
251       histoTh->fill(lgeps, NbProcess);            239       histoTh->fill(lgeps, NbProcess);
252     }                                             240     }
253   }                                               241   }
254                                                << 242      
255   // return integrated crossSection            << 243   //return integrated crossSection
256   //                                              244   //
257   return sigmaTot;                             << 245   return sigmaTot;   
258 }                                                 246 }
259                                                   247 
260 //....oooOO0OOooo........oooOO0OOooo........oo    248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
261                                                   249 
262 void RunAction::FillCrossSectionHisto(const G4    250 void RunAction::FillCrossSectionHisto(const G4String& process, G4int)
263 {                                                 251 {
264   const G4Material* material = fDetector->GetM    252   const G4Material* material = fDetector->GetMaterial();
265   G4double ekin = fPrimary->GetParticleGun()->    253   G4double ekin = fPrimary->GetParticleGun()->GetParticleEnergy();
266   G4ParticleDefinition* particle = fPrimary->G << 254   G4ParticleDefinition *particle = fPrimary->GetParticleGun()->GetParticleDefinition();
267   G4double particleMass = particle->GetPDGMass    255   G4double particleMass = particle->GetPDGMass();
268                                                << 256   
269   G4EmCalculator emCal;                           257   G4EmCalculator emCal;
270                                                   258 
271   G4int id = 0;                                << 259   G4int id = 0; G4double cut = 1.e-10*ekin;
272   G4double cut = 1.e-10 * ekin;                << 260   if (process == "muIoni")          {id = 21; cut = GetEnergyCut(material,1);}
273   if (process == "muIoni") {                   << 261   else if (process == "muPairProd") {id = 22; cut = 2*(GetEnergyCut(material,1)
274     id = 21;                                   << 262                                                       + electron_mass_c2); }
275     cut = GetEnergyCut(material, 1);           << 263   else if (process == "muBrems")    {id = 23; cut = GetEnergyCut(material,0);}
276   }                                            << 264   else if (process == "muonNuclear"){id = 24; cut = 100*MeV;}
277   else if (process == "muPairProd") {          << 265   else if (process == "muToMuonPairProd"){id = 28; cut = 2*particleMass;}
278     id = 22;                                   << 266   if (id == 0) { return; }
279     cut = 2 * (GetEnergyCut(material, 1) + ele << 
280   }                                            << 
281   else if (process == "muBrems") {             << 
282     id = 23;                                   << 
283     cut = GetEnergyCut(material, 0);           << 
284   }                                            << 
285   else if (process == "muonNuclear") {         << 
286     id = 24;                                   << 
287     cut = 100 * MeV;                           << 
288   }                                            << 
289   else if (process == "muToMuonPairProd") {    << 
290     id = 28;                                   << 
291     cut = 2 * particleMass;                    << 
292   }                                            << 
293   if (id == 0) {                               << 
294     return;                                    << 
295   }                                            << 
296                                                   267 
297   G4int nbOfBins = 100;                           268   G4int nbOfBins = 100;
298   G4double binMin = cut;                          269   G4double binMin = cut;
299   G4double binMax = ekin;                         270   G4double binMax = ekin;
300   G4double binWidth = (binMax - binMin) / G4do << 271   G4double binWidth = (binMax-binMin)/G4double(nbOfBins);
301                                                   272 
302   G4AnalysisManager* analysisManager = G4Analy    273   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
303                                                << 274     
304   G4H1* histoTh = 0;                              275   G4H1* histoTh = 0;
305   if (fHistoManager->HistoExist(id)) {            276   if (fHistoManager->HistoExist(id)) {
306     histoTh = analysisManager->GetH1(fHistoMan << 277     histoTh  = analysisManager->GetH1(fHistoManager->GetHistoID(id));  
307     nbOfBins = fHistoManager->GetNbins(id);       278     nbOfBins = fHistoManager->GetNbins(id);
308     binMin = fHistoManager->GetVmin(id);       << 279     binMin   = fHistoManager->GetVmin (id);
309     binMax = fHistoManager->GetVmax(id);       << 280     binMax   = fHistoManager->GetVmax (id);
310     binWidth = fHistoManager->GetBinWidth(id); << 281     binWidth = fHistoManager->GetBinWidth(id);    
311   }                                               282   }
312                                                   283 
313   G4double sigma, primaryEnergy;                  284   G4double sigma, primaryEnergy;
314                                                   285 
315   for (G4int ibin = 0; ibin < nbOfBins; ibin++ << 286   for(G4int ibin=0; ibin<nbOfBins; ibin++){
316     primaryEnergy = binMin + (ibin + 0.5) * bi << 287     primaryEnergy = binMin + (ibin+0.5)*binWidth;
317     sigma = emCal.GetCrossSectionPerVolume(pri    288     sigma = emCal.GetCrossSectionPerVolume(primaryEnergy, particle, process, material);
318     if (histoTh) {                                289     if (histoTh) {
319       histoTh->fill(primaryEnergy, sigma);        290       histoTh->fill(primaryEnergy, sigma);
320     }                                             291     }
321   }                                               292   }
322 }                                                 293 }
323                                                   294 
324 //....oooOO0OOooo........oooOO0OOooo........oo    295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
325                                                   296 
326 G4double RunAction::GetEnergyCut(const G4Mater    297 G4double RunAction::GetEnergyCut(const G4Material* material, G4int idParticle)
327 {                                              << 298 { 
328   G4ProductionCutsTable* table = G4ProductionC << 299  G4ProductionCutsTable* table = G4ProductionCutsTable::GetProductionCutsTable();
329                                                << 300  
330   size_t index = 0;                            << 301  size_t index = 0;
331   while ((table->GetMaterialCutsCouple(index)- << 302  while ( (table->GetMaterialCutsCouple(index)->GetMaterial() != material) &&
332          && (index < table->GetTableSize()))   << 303         (index < table->GetTableSize())) index++;
333     index++;                                   << 
334                                                   304 
335   return (*(table->GetEnergyCutsVector(idParti << 305  return (*(table->GetEnergyCutsVector(idParticle)))[index];
336 }                                              << 306 } 
337                                                   307 
338 //....oooOO0OOooo........oooOO0OOooo........oo    308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 309                    
339                                                   310