Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm5/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/TestEm5/src/RunAction.cc (Version 11.3.0) and /examples/extended/electromagnetic/TestEm5/src/RunAction.cc (Version 9.1.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/TestEm5/src/RunActio <<  26 // $Id: RunAction.cc,v 1.25 2007/11/21 17:41:19 maire Exp $
 27 /// \brief Implementation of the RunAction cla <<  27 // GEANT4 tag $Name: geant4-09-01-patch-01 $
 28 //                                             << 
 29 //                                                 28 //
 30 //....oooOO0OOooo........oooOO0OOooo........oo     29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oo     30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32                                                    31 
 33 #include "RunAction.hh"                            32 #include "RunAction.hh"
 34                                                << 
 35 #include "DetectorConstruction.hh"                 33 #include "DetectorConstruction.hh"
 36 #include "HistoManager.hh"                     << 
 37 #include "PrimaryGeneratorAction.hh"               34 #include "PrimaryGeneratorAction.hh"
 38 #include "Run.hh"                              <<  35 #include "HistoManager.hh"
 39                                                    36 
 40 #include "G4EmCalculator.hh"                   << 
 41 #include "G4Run.hh"                                37 #include "G4Run.hh"
 42 #include "G4SystemOfUnits.hh"                  <<  38 #include "G4RunManager.hh"
 43 #include "G4UnitsTable.hh"                         39 #include "G4UnitsTable.hh"
 44 #include "Randomize.hh"                        <<  40 #include "G4EmCalculator.hh"
 45                                                    41 
                                                   >>  42 #include "Randomize.hh"
 46 #include <iomanip>                                 43 #include <iomanip>
 47                                                    44 
 48 //....oooOO0OOooo........oooOO0OOooo........oo     45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 49                                                    46 
 50 RunAction::RunAction(DetectorConstruction* det <<  47 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin,
 51   : fDetector(det), fPrimary(kin)              <<  48                      HistoManager* histo)
 52 {                                              <<  49 :detector(det), primary(kin), histoManager(histo)
 53   // Book predefined histograms                <<  50 { }
 54   fHistoManager = new HistoManager();          << 
 55 }                                              << 
 56                                                    51 
 57 //....oooOO0OOooo........oooOO0OOooo........oo     52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 58                                                    53 
 59 RunAction::~RunAction()                            54 RunAction::~RunAction()
 60 {                                              <<  55 { }
 61   delete fHistoManager;                        << 
 62 }                                              << 
 63                                                    56 
 64 //....oooOO0OOooo........oooOO0OOooo........oo     57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 65                                                    58 
 66 G4Run* RunAction::GenerateRun()                <<  59 void RunAction::BeginOfRunAction(const G4Run* aRun)
 67 {                                                  60 {
 68   fRun = new Run(fDetector);                   <<  61   G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
 69   return fRun;                                 <<  62 
                                                   >>  63   //initialisation
                                                   >>  64   EnergyDeposit  = EnergyDeposit2  = 0.;
                                                   >>  65   TrakLenCharged = TrakLenCharged2 = 0.;
                                                   >>  66   TrakLenNeutral = TrakLenNeutral2 = 0.;
                                                   >>  67   nbStepsCharged = nbStepsCharged2 = 0.;
                                                   >>  68   nbStepsNeutral = nbStepsNeutral2 = 0.;
                                                   >>  69   MscProjecTheta = MscProjecTheta2 = 0.;
                                                   >>  70   MscThetaCentral = 3*ComputeMscHighland();
                                                   >>  71 
                                                   >>  72   nbGamma = nbElect = nbPosit = 0;
                                                   >>  73 
                                                   >>  74   Transmit[0] = Transmit[1] = Reflect[0] = Reflect[1] = 0;
                                                   >>  75   
                                                   >>  76   MscEntryCentral = 0;
                                                   >>  77   
                                                   >>  78   EnergyLeak[0] = EnergyLeak[1] = EnergyLeak2[0] = EnergyLeak2[1] = 0.;
                                                   >>  79 
                                                   >>  80   histoManager->book();
                                                   >>  81 
                                                   >>  82   // save Rndm status
                                                   >>  83   G4RunManager::GetRunManager()->SetRandomNumberStore(true);
                                                   >>  84   CLHEP::HepRandom::showEngineStatus();
 70 }                                                  85 }
 71                                                    86 
 72 //....oooOO0OOooo........oooOO0OOooo........oo     87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 73                                                    88 
 74 void RunAction::BeginOfRunAction(const G4Run*) <<  89 void RunAction::EndOfRunAction(const G4Run* aRun)
 75 {                                                  90 {
 76   // show Rndm status                          <<  91   // compute mean and rms
 77   if (isMaster) G4Random::showEngineStatus();  <<  92   //
 78                                                <<  93   G4int TotNbofEvents = aRun->GetNumberOfEvent();
 79   // keep run condition                        <<  94   if (TotNbofEvents == 0) return;
 80   if (fPrimary) {                              <<  95   
 81     G4ParticleDefinition* particle = fPrimary- <<  96   G4double EnergyBalance = EnergyDeposit + EnergyLeak[0] + EnergyLeak[1];
 82     G4double energy = fPrimary->GetParticleGun <<  97   EnergyBalance /= TotNbofEvents;
 83     fRun->SetPrimary(particle, energy);        <<  98 
                                                   >>  99   EnergyDeposit /= TotNbofEvents; EnergyDeposit2 /= TotNbofEvents;
                                                   >> 100   G4double rmsEdep = EnergyDeposit2 - EnergyDeposit*EnergyDeposit;
                                                   >> 101   if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep/TotNbofEvents);
                                                   >> 102   else            rmsEdep = 0.;
                                                   >> 103 
                                                   >> 104   TrakLenCharged /= TotNbofEvents; TrakLenCharged2 /= TotNbofEvents;
                                                   >> 105   G4double rmsTLCh = TrakLenCharged2 - TrakLenCharged*TrakLenCharged;
                                                   >> 106   if (rmsTLCh>0.) rmsTLCh = std::sqrt(rmsTLCh/TotNbofEvents);
                                                   >> 107   else            rmsTLCh = 0.;
                                                   >> 108 
                                                   >> 109   TrakLenNeutral /= TotNbofEvents; TrakLenNeutral2 /= TotNbofEvents;
                                                   >> 110   G4double rmsTLNe = TrakLenNeutral2 - TrakLenNeutral*TrakLenNeutral;
                                                   >> 111   if (rmsTLNe>0.) rmsTLNe = std::sqrt(rmsTLNe/TotNbofEvents);
                                                   >> 112   else            rmsTLNe = 0.;
                                                   >> 113 
                                                   >> 114   nbStepsCharged /= TotNbofEvents; nbStepsCharged2 /= TotNbofEvents;
                                                   >> 115   G4double rmsStCh = nbStepsCharged2 - nbStepsCharged*nbStepsCharged;
                                                   >> 116   if (rmsStCh>0.) rmsStCh = std::sqrt(rmsTLCh/TotNbofEvents);
                                                   >> 117   else            rmsStCh = 0.;
                                                   >> 118 
                                                   >> 119   nbStepsNeutral /= TotNbofEvents; nbStepsNeutral2 /= TotNbofEvents;
                                                   >> 120   G4double rmsStNe = nbStepsNeutral2 - nbStepsNeutral*nbStepsNeutral;
                                                   >> 121   if (rmsStNe>0.) rmsStNe = std::sqrt(rmsTLCh/TotNbofEvents);
                                                   >> 122   else            rmsStNe = 0.;
                                                   >> 123 
                                                   >> 124   G4double Gamma = (double)nbGamma/TotNbofEvents;
                                                   >> 125   G4double Elect = (double)nbElect/TotNbofEvents;
                                                   >> 126   G4double Posit = (double)nbPosit/TotNbofEvents;
                                                   >> 127 
                                                   >> 128   G4double transmit[2];
                                                   >> 129   transmit[0] = 100.*Transmit[0]/TotNbofEvents;
                                                   >> 130   transmit[1] = 100.*Transmit[1]/TotNbofEvents;
                                                   >> 131 
                                                   >> 132   G4double reflect[2];
                                                   >> 133   reflect[0] = 100.*Reflect[0]/TotNbofEvents;
                                                   >> 134   reflect[1] = 100.*Reflect[1]/TotNbofEvents;
                                                   >> 135 
                                                   >> 136   G4double rmsMsc = 0., tailMsc = 0.;
                                                   >> 137   if (MscEntryCentral > 0) {
                                                   >> 138     MscProjecTheta /= MscEntryCentral; MscProjecTheta2 /= MscEntryCentral;
                                                   >> 139     rmsMsc = MscProjecTheta2 - MscProjecTheta*MscProjecTheta;
                                                   >> 140     if (rmsMsc > 0.) rmsMsc = std::sqrt(rmsMsc);
                                                   >> 141     tailMsc = 100.- (100.*MscEntryCentral)/(2*Transmit[1]);    
 84   }                                               142   }
 85                                                << 143   
 86   // histograms                                << 144   EnergyLeak[0] /= TotNbofEvents; EnergyLeak2[0] /= TotNbofEvents;
                                                   >> 145   G4double rmsEl0 = EnergyLeak2[0] - EnergyLeak[0]*EnergyLeak[0];
                                                   >> 146   if (rmsEl0>0.) rmsEl0 = std::sqrt(rmsEl0/TotNbofEvents);
                                                   >> 147   else           rmsEl0 = 0.;
                                                   >> 148   
                                                   >> 149   EnergyLeak[1] /= TotNbofEvents; EnergyLeak2[1] /= TotNbofEvents;
                                                   >> 150   G4double rmsEl1 = EnergyLeak2[1] - EnergyLeak[1]*EnergyLeak[1];
                                                   >> 151   if (rmsEl1>0.) rmsEl1 = std::sqrt(rmsEl1/TotNbofEvents);
                                                   >> 152   else           rmsEl1 = 0.;    
                                                   >> 153   
                                                   >> 154       
                                                   >> 155   //Stopping Power from input Table.
 87   //                                              156   //
 88   G4AnalysisManager* analysisManager = G4Analy << 157   G4Material* material = detector->GetAbsorberMaterial();
 89   if (analysisManager->IsActive()) {           << 158   G4double length  = detector->GetAbsorberThickness();
 90     analysisManager->OpenFile();               << 159   G4double density = material->GetDensity();
                                                   >> 160    
                                                   >> 161   G4ParticleDefinition* particle = primary->GetParticleGun()
                                                   >> 162                                           ->GetParticleDefinition();
                                                   >> 163   G4String partName = particle->GetParticleName();
                                                   >> 164   G4double energy = primary->GetParticleGun()->GetParticleEnergy();
                                                   >> 165 
                                                   >> 166   G4EmCalculator emCalculator;
                                                   >> 167   G4double dEdxTable = 0., dEdxFull = 0.;
                                                   >> 168   if (particle->GetPDGCharge()!= 0.) { 
                                                   >> 169     dEdxTable = emCalculator.GetDEDX(energy,particle,material);
                                                   >> 170     dEdxFull  = emCalculator.ComputeTotalDEDX(energy,particle,material);    
 91   }                                               171   }
 92 }                                              << 172   G4double stopTable = dEdxTable/density;
                                                   >> 173   G4double stopFull  = dEdxFull /density; 
                                                   >> 174    
                                                   >> 175   //Stopping Power from simulation.
                                                   >> 176   //    
                                                   >> 177   G4double meandEdx  = EnergyDeposit/length;
                                                   >> 178   G4double stopPower = meandEdx/density;  
                                                   >> 179 
                                                   >> 180   G4cout << "\n ======================== run summary ======================\n";
                                                   >> 181 
                                                   >> 182   G4int prec = G4cout.precision(3);
                                                   >> 183   
                                                   >> 184   G4cout << "\n The run was " << TotNbofEvents << " " << partName << " of "
                                                   >> 185          << G4BestUnit(energy,"Energy") << " through " 
                                                   >> 186    << G4BestUnit(length,"Length") << " of "
                                                   >> 187    << material->GetName() << " (density: " 
                                                   >> 188    << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
                                                   >> 189   
                                                   >> 190   G4cout.precision(4);
                                                   >> 191   
                                                   >> 192   G4cout << "\n Total energy deposit in absorber per event = "
                                                   >> 193          << G4BestUnit(EnergyDeposit,"Energy") << " +- "
                                                   >> 194          << G4BestUnit(rmsEdep,      "Energy") 
                                                   >> 195          << G4endl;
                                                   >> 196    
                                                   >> 197   G4cout << "\n -----> Mean dE/dx = " << meandEdx/(MeV/cm) << " MeV/cm"
                                                   >> 198          << "\t(" << stopPower/(MeV*cm2/g) << " MeV*cm2/g)"
                                                   >> 199    << G4endl;
                                                   >> 200    
                                                   >> 201   G4cout << "\n From formulas :" << G4endl; 
                                                   >> 202   G4cout << "   restricted dEdx = " << dEdxTable/(MeV/cm) << " MeV/cm"
                                                   >> 203          << "\t(" << stopTable/(MeV*cm2/g) << " MeV*cm2/g)"
                                                   >> 204    << G4endl;
                                                   >> 205    
                                                   >> 206   G4cout << "   full dEdx       = " << dEdxFull/(MeV/cm) << " MeV/cm"
                                                   >> 207          << "\t(" << stopFull/(MeV*cm2/g) << " MeV*cm2/g)"
                                                   >> 208    << G4endl;
                                                   >> 209    
                                                   >> 210   G4cout << "\n Leakage :  primary = "
                                                   >> 211          << G4BestUnit(EnergyLeak[0],"Energy") << " +- "
                                                   >> 212          << G4BestUnit(rmsEl0,       "Energy")
                                                   >> 213    << "   secondaries = "
                                                   >> 214          << G4BestUnit(EnergyLeak[1],"Energy") << " +- "
                                                   >> 215          << G4BestUnit(rmsEl1,       "Energy")    
                                                   >> 216          << G4endl;
                                                   >> 217    
                                                   >> 218   G4cout << " Energy balance :  edep + eleak = "
                                                   >> 219          << G4BestUnit(EnergyBalance,"Energy")
                                                   >> 220          << G4endl;  
                                                   >> 221        
                                                   >> 222   G4cout << "\n Total track length (charged) in absorber per event = "
                                                   >> 223          << G4BestUnit(TrakLenCharged,"Length") << " +- "
                                                   >> 224          << G4BestUnit(rmsTLCh,       "Length") << G4endl;
                                                   >> 225 
                                                   >> 226   G4cout << " Total track length (neutral) in absorber per event = "
                                                   >> 227          << G4BestUnit(TrakLenNeutral,"Length") << " +- "
                                                   >> 228          << G4BestUnit(rmsTLNe,       "Length") << G4endl;
                                                   >> 229 
                                                   >> 230   G4cout << "\n Number of steps (charged) in absorber per event = "
                                                   >> 231          << nbStepsCharged << " +- " << rmsStCh << G4endl;
                                                   >> 232 
                                                   >> 233   G4cout << " Number of steps (neutral) in absorber per event = "
                                                   >> 234          << nbStepsNeutral << " +- " << rmsStNe << G4endl;
                                                   >> 235 
                                                   >> 236   G4cout << "\n Number of secondaries per event : Gammas = " << Gamma
                                                   >> 237          << ";   electrons = " << Elect
                                                   >> 238      << ";   positrons = " << Posit << G4endl;
                                                   >> 239 
                                                   >> 240   G4cout << "\n Number of events with the primary particle transmitted = "
                                                   >> 241          << transmit[1] << " %" << G4endl;
                                                   >> 242 
                                                   >> 243   G4cout << " Number of events with at least  1 particle transmitted "
                                                   >> 244          << "(same charge as primary) = " << transmit[0] << " %" << G4endl;
 93                                                   245 
 94 //....oooOO0OOooo........oooOO0OOooo........oo << 246   G4cout << "\n Number of events with the primary particle reflected = "
                                                   >> 247          << reflect[1] << " %" << G4endl;
 95                                                   248 
 96 void RunAction::EndOfRunAction(const G4Run*)   << 249   G4cout << " Number of events with at least  1 particle reflected "
 97 {                                              << 250          << "(same charge as primary) = " << reflect[0] << " %" << G4endl;
 98   // print Run summary                         << 
 99   //                                           << 
100   if (isMaster) fRun->EndOfRun();              << 
101                                                   251 
102   // save histograms                           << 252   // compute width of the Gaussian central part of the MultipleScattering
103   G4AnalysisManager* analysisManager = G4Analy << 253   //
104   if (analysisManager->IsActive()) {           << 254   if (histoManager->HistoExist(13)) {
105     analysisManager->Write();                  << 255     G4cout << "\n MultipleScattering:" 
106     analysisManager->CloseFile();              << 256            << "\n  rms proj angle of transmit primary particle = "
                                                   >> 257            << rmsMsc/mrad << " mrad (central part only)" << G4endl;
                                                   >> 258 
                                                   >> 259     G4cout << "  computed theta0 (Highland formula)          = "
                                                   >> 260      << ComputeMscHighland()/mrad << " mrad" << G4endl;
                                                   >> 261      
                                                   >> 262     G4cout << "  central part defined as +- "
                                                   >> 263      << MscThetaCentral/mrad << " mrad; " 
                                                   >> 264      << "  Tail ratio = " << tailMsc << " %" << G4endl;    
107   }                                               265   }
108                                                   266 
                                                   >> 267   G4cout.precision(prec);
                                                   >> 268 
                                                   >> 269   histoManager->save();
                                                   >> 270 
109   // show Rndm status                             271   // show Rndm status
110   if (isMaster) G4Random::showEngineStatus();  << 272   CLHEP::HepRandom::showEngineStatus();
                                                   >> 273 }
                                                   >> 274 
                                                   >> 275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 276 
                                                   >> 277 G4double RunAction::ComputeMscHighland()
                                                   >> 278 {
                                                   >> 279   //compute the width of the Gaussian central part of the MultipleScattering
                                                   >> 280   //projected angular distribution.
                                                   >> 281   //Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
                                                   >> 282 
                                                   >> 283   G4double t = (detector->GetAbsorberThickness())
                                                   >> 284               /(detector->GetAbsorberMaterial()->GetRadlen());
                                                   >> 285   if (t < DBL_MIN) return 0.;
                                                   >> 286 
                                                   >> 287   G4ParticleGun* particle = primary->GetParticleGun();
                                                   >> 288   G4double T = particle->GetParticleEnergy();
                                                   >> 289   G4double M = particle->GetParticleDefinition()->GetPDGMass();
                                                   >> 290   G4double z = std::abs(particle->GetParticleDefinition()->GetPDGCharge()/eplus);
                                                   >> 291 
                                                   >> 292   G4double bpc = T*(T+2*M)/(T+M);
                                                   >> 293   G4double teta0 = 13.6*MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc;
                                                   >> 294   return teta0;
111 }                                                 295 }
112                                                   296 
113 //....oooOO0OOooo........oooOO0OOooo........oo    297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114                                                   298