Geant4 Cross Reference

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


  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/TestEm6/src/RunActio     26 /// \file electromagnetic/TestEm6/src/RunAction.cc
 27 /// \brief Implementation of the RunAction cla     27 /// \brief Implementation of the RunAction class
 28 //                                                 28 //
 29 //                                             <<  29 // $Id$
                                                   >>  30 // 
 30 //....oooOO0OOooo........oooOO0OOooo........oo     31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oo     32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32                                                    33 
 33 #include "RunAction.hh"                        <<  34 #include "G4SystemOfUnits.hh"  
                                                   >>  35 #include "G4PhysicalConstants.hh"  
                                                   >>  36 #include "G4EmCalculator.hh" 
                                                   >>  37 #include "G4ParticleTable.hh" 
                                                   >>  38 #include "G4ParticleDefinition.hh" 
                                                   >>  39 #include "G4Positron.hh" 
                                                   >>  40 #include "G4AnnihiToMuPair.hh"  
                                                   >>  41 #include "G4eeToHadrons.hh"  
                                                   >>  42 #include "G4eeToHadronsModel.hh"  
                                                   >>  43 #include "G4eBremsstrahlung.hh"  
 34                                                    44 
 35 #include "G4AnnihiToMuPair.hh"                 <<  45 #include "RunAction.hh"
 36 #include "G4EmCalculator.hh"                   << 
 37 #include "G4MuonMinus.hh"                      << 
 38 #include "G4ParticleDefinition.hh"             << 
 39 #include "G4ParticleTable.hh"                  << 
 40 #include "G4PhysicalConstants.hh"              << 
 41 #include "G4Positron.hh"                       << 
 42 #include "G4Run.hh"                                46 #include "G4Run.hh"
 43 #include "G4RunManager.hh"                         47 #include "G4RunManager.hh"
 44 #include "G4SystemOfUnits.hh"                  <<  48 #include "G4MuonMinus.hh"
 45 #include "G4eBremsstrahlung.hh"                <<  49 
 46 #include "G4eeToHadrons.hh"                    << 
 47 #include "G4eeToHadronsModel.hh"               << 
 48 #include "Randomize.hh"                            50 #include "Randomize.hh"
 49                                                    51 
 50 #include <sstream>                                 52 #include <sstream>
 51                                                    53 
 52 //....oooOO0OOooo........oooOO0OOooo........oo     54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 53                                                    55 
 54 RunAction::RunAction(DetectorConstruction* det     56 RunAction::RunAction(DetectorConstruction* det)
 55   : G4UserRunAction(), fDetector(det), fProcCo <<  57  : G4UserRunAction(),fDetector(det),fProcCounter(0),fAnalysis(0),fMat(0)
 56 {                                                  58 {
 57   fMinE = 40 * GeV;                            <<  59   fMinE = 40*GeV;
 58   fMaxE = 10000 * GeV;                         <<  60   fMaxE = 10000*GeV;
 59   fnBin = 10000;                                   61   fnBin = 10000;
                                                   >>  62   fNtColId[0] = fNtColId[1] = fNtColId[2] = 0;
 60 }                                                  63 }
 61                                                    64 
 62 //....oooOO0OOooo........oooOO0OOooo........oo     65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 63                                                    66 
 64 RunAction::~RunAction() {}                     <<  67 RunAction::~RunAction()
                                                   >>  68 {}
 65                                                    69 
 66 //....oooOO0OOooo........oooOO0OOooo........oo     70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 67                                                    71 
 68 void RunAction::BeginOfRunAction(const G4Run*      72 void RunAction::BeginOfRunAction(const G4Run* aRun)
 69 {                                              <<  73 {  
 70   G4cout << "### Run " << aRun->GetRunID() <<      74   G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
 71                                                    75 
 72   // get material                              <<  76   //get material
 73   //                                               77   //
 74   fMat = fDetector->GetMaterial();                 78   fMat = fDetector->GetMaterial();
 75   G4cout << "###RunAction::BeginOfRunAction  M <<  79   G4cout << "###RunAction::BeginOfRunAction  Material:" 
 76                                                <<  80          << fMat->GetName() << G4endl;
                                                   >>  81  
 77   fProcCounter = new ProcessesCount;               82   fProcCounter = new ProcessesCount;
 78                                                    83 
 79   fAnalysis = G4AnalysisManager::Instance();       84   fAnalysis = G4AnalysisManager::Instance();
 80   fAnalysis->SetDefaultFileType("root");       <<  85   
 81                                                << 
 82   // Open an output file                           86   // Open an output file
 83   //                                               87   //
 84   std::stringstream tmp;                           88   std::stringstream tmp;
 85   tmp << "testem6_" << aRun->GetRunID();           89   tmp << "testem6_" << aRun->GetRunID();
 86   G4String fileName = tmp.str();                   90   G4String fileName = tmp.str();
 87   fAnalysis->OpenFile(fileName);               <<  91   fAnalysis->OpenFile(fileName);    
 88   fAnalysis->SetVerboseLevel(2);                   92   fAnalysis->SetVerboseLevel(2);
 89   fAnalysis->SetActivation(true);              <<  93   G4String extension = fAnalysis->GetFileType();
 90                                                <<  94   fileName = fileName + "." + extension;
                                                   >>  95     
 91   // Creating histograms 1,2,3,4,5,6               96   // Creating histograms 1,2,3,4,5,6
 92   //                                           <<  97   // 
 93   fAnalysis->SetFirstHistoId(1);               <<  98   fAnalysis->SetFirstHistoId(1);  
 94   fAnalysis->CreateH1("h1", "1/(1+(theta+[g]+) <<  99   fAnalysis->CreateH1("h1","1/(1+(theta+[g]+)**2)",100, 0 ,1.);
 95   fAnalysis->CreateH1("h2", "log10(theta+ [g]+ << 100   fAnalysis->CreateH1("h2","log10(theta+ [g]+)",   100,-3.,1.);
 96   fAnalysis->CreateH1("h3", "log10(theta- [g]- << 101   fAnalysis->CreateH1("h3","log10(theta- [g]-)",   100,-3.,1.);
 97   fAnalysis->CreateH1("h4", "log10(theta+ [g]+ << 102   fAnalysis->CreateH1("h4","log10(theta+ [g]+ -theta- [g]-)", 100,-3.,1.);
 98   fAnalysis->CreateH1("h5", "xPlus", 100, 0.,  << 103   fAnalysis->CreateH1("h5","xPlus" ,100,0.,1.);
 99   fAnalysis->CreateH1("h6", "xMinus", 100, 0., << 104   fAnalysis->CreateH1("h6","xMinus",100,0.,1.);
100                                                << 105   
101   // creating histogram 7,8,9,10,11 (CrossSect << 106   //creating histogram 7,8,9,10,11 (CrossSectionPerAtom)
102   //                                           << 107   //
103   G4double minBin = std::log10(fMinE / GeV);   << 108   G4double minBin = std::log10(fMinE/GeV);
104   G4double maxBin = std::log10(fMaxE / GeV);   << 109   G4double maxBin = std::log10(fMaxE/GeV);
105   fAnalysis->CreateH1("h7", "CrossSectionPerAt << 110   fAnalysis->CreateH1("h7","CrossSectionPerAtom of AnnihiToMuMu (microbarn)",
106                       maxBin);                 << 111                       fnBin,minBin,maxBin);
107   fAnalysis->CreateH1("h8", "CrossSectionPerAt << 112   fAnalysis->CreateH1("h8",
108                       maxBin);                 << 113     "CrossSectionPerAtom of AnnihiToTwoGamma (microbarn)",fnBin,minBin,maxBin);
109   fAnalysis->CreateH1("h9", "CrossSectionPerAt << 114   fAnalysis->CreateH1("h9","CrossSectionPerAtom of AnnihiToHadrons (microbarn)",
110                       maxBin);                 << 115                       fnBin,minBin,maxBin);
111   fAnalysis->CreateH1("h10", "Theoretical Cros << 116   fAnalysis->CreateH1("h10",
112                       fnBin, minBin, maxBin);  << 117     "Theoretical CrossSectionPerAtom of AnnihiToTwoGamma (microbarn)",
113   fAnalysis->CreateH1("h11", "Theoretical Cros << 118                       fnBin,minBin,maxBin);
114                       minBin, maxBin);         << 119   fAnalysis->CreateH1("h11",
115                                                << 120     "Theoretical CrossSectionPerAtom of AnnihiToMuMu (microbarn)",
116   // creating histogram 12,13,14,15,16(CrossSe << 121                       fnBin,minBin,maxBin);
117   //                                           << 122 
118   fAnalysis->CreateH1("h12", "CrossSectionPerV << 123   //creating histogram 12,13,14,15,16(CrossSectionPerVolume)
119   fAnalysis->CreateH1("h13", "CrossSectionPerV << 124   //
120   fAnalysis->CreateH1("h14", "CrossSectionPerV << 125   fAnalysis->CreateH1("h12","CrossSectionPerVol of Bremsstraulung (1/mm) ",
121   fAnalysis->CreateH1("h15", "CrossSectionPerV << 126                       fnBin,minBin,maxBin);
122                       maxBin);                 << 127   fAnalysis->CreateH1("h13","CrossSectionPerVol of Ionization (1/mm)",
123   fAnalysis->CreateH1("h16", "CrossSectionPerV << 128                       fnBin,minBin,maxBin);
124                                                << 129   fAnalysis->CreateH1("h14","CrossSectionPerVol of AnnihiToMuMu (1/mm)",
125   // creating histogram 17 (R ratio)           << 130                       fnBin,minBin,maxBin);
126   fAnalysis->CreateH1("h17", "R : eeToHadr/eeT << 131   fAnalysis->CreateH1("h15","CrossSectionPerVol of AnnihiToTwoGamma (1/mm)",
127                                                << 132                       fnBin,minBin,maxBin);
128   G4cout << "\n----> Histogram file is opened  << 133   fAnalysis->CreateH1("h16","CrossSectionPerVol of AnnihiToHadrons (1/mm)",
                                                   >> 134                       fnBin,minBin,maxBin);
                                                   >> 135   
                                                   >> 136   //creating histogram 17 (R ratio)
                                                   >> 137   fAnalysis->CreateH1("h17","R : eeToHadr/eeToMu",fnBin,minBin,maxBin);
                                                   >> 138     
                                                   >> 139   G4cout << "\n----> Histogram file is opened in " << fileName << G4endl;   
129 }                                                 140 }
130                                                   141 
131 //....oooOO0OOooo........oooOO0OOooo........oo    142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132                                                   143 
133 void RunAction::CountProcesses(G4String procNa    144 void RunAction::CountProcesses(G4String procName)
134 {                                                 145 {
135   // does the process  already encounted ?     << 146   //does the process  already encounted ?
136   //                                              147   //
137   size_t nbProc = fProcCounter->size();           148   size_t nbProc = fProcCounter->size();
138   size_t i = 0;                                   149   size_t i = 0;
139   while ((i < nbProc) && ((*fProcCounter)[i]-> << 150   while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++;
140     i++;                                       << 151   if (i == nbProc) fProcCounter->push_back( new OneProcessCount(procName));
141   if (i == nbProc) fProcCounter->push_back(new << 152   
142                                                << 
143   (*fProcCounter)[i]->Count();                    153   (*fProcCounter)[i]->Count();
144 }                                                 154 }
145                                                   155 
146 //....oooOO0OOooo........oooOO0OOooo........oo    156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147                                                   157 
148 void RunAction::EndOfRunAction(const G4Run*)      158 void RunAction::EndOfRunAction(const G4Run*)
149 {                                                 159 {
150   G4cout << "### RunAction::EndOfRunAction" << << 160   // show Rndm status : not needed
151   // total number of process calls             << 161   //
                                                   >> 162   //CLHEP::HepRandom::showEngineStatus();
                                                   >> 163 
                                                   >> 164   //total number of process calls
152   //                                              165   //
                                                   >> 166   G4double countTot = 0.;
153   G4cout << "\n Number of process calls --->";    167   G4cout << "\n Number of process calls --->";
154   for (size_t i = 0; i < fProcCounter->size(); << 168   for (size_t i=0; i< fProcCounter->size();i++) {
155     G4String procName = (*fProcCounter)[i]->Ge << 169         G4String procName = (*fProcCounter)[i]->GetName();
156     if (procName != "Transportation") {        << 170         if (procName != "Transportation") {
157       G4int count = (*fProcCounter)[i]->GetCou << 171           G4int count    = (*fProcCounter)[i]->GetCounter(); 
158       G4cout << "\t" << procName << " : " << c << 172           G4cout << "\t" << procName << " : " << count;
159     }                                          << 173           countTot += count;
                                                   >> 174         }
160   }                                               175   }
161   G4cout << G4endl;                            << 176   
162                                                << 177   //instance EmCalculator
163   // instance EmCalculator                     << 
164   //                                              178   //
165   G4EmCalculator emCal;                           179   G4EmCalculator emCal;
166   emCal.SetVerbose(0);                            180   emCal.SetVerbose(0);
167                                                   181 
168   // get positron                              << 182   //get positron
169   //                                              183   //
170   G4String positronName = "e+";                   184   G4String positronName = "e+";
171   G4ParticleDefinition* positron = G4ParticleT << 185   G4ParticleDefinition* positron = 
                                                   >> 186     G4ParticleTable::GetParticleTable()->FindParticle(positronName);
172                                                   187 
173   // process name                              << 188   //process name
174   //                                              189   //
175   G4String annihilName = "annihil";            << 190   G4String annihilName      = "annihil";
176   G4String annihiToMuName = "AnnihiToMuPair";  << 191   G4String annihiToMuName   = "AnnihiToMuPair";
177   G4String annihiToHadrName = "ee2hadr";          192   G4String annihiToHadrName = "ee2hadr";
178   G4String BremName = "eBrem";                 << 193   G4String BremName         = "eBrem";
179   G4String IoniName = "eIoni";                 << 194   G4String IoniName         = "eIoni";
                                                   >> 195   
                                                   >> 196 
                                                   >> 197   //get AnnihiToMuPair
                                                   >> 198   //
                                                   >> 199   G4AnnihiToMuPair* annihiToMu = 
                                                   >> 200     reinterpret_cast<G4AnnihiToMuPair*>(emCal.FindProcess(positron,
                                                   >> 201                                                           annihiToMuName));
                                                   >> 202 
                                                   >> 203   G4double de, x, energy;
                                                   >> 204   G4double crs_annihil, crs_annihiToMu, crs_annihiToHadr;
                                                   >> 205   G4double crsVol_annihil, crsVol_annihiToMu, crsVol_annihiToHadr,
                                                   >> 206     crsVol_Brem, crsVol_Ioni;
                                                   >> 207   G4double crs_annihil_theory, crs_annihiToMu_theory, RR;
                                                   >> 208   G4double X1, X2;
                                                   >> 209 
                                                   >> 210   G4double atomicZ, atomicA;
                                                   >> 211 
                                                   >> 212   G4double Mu;  //for muon mass
                                                   >> 213   G4double Me;  //for electron mass
                                                   >> 214   G4double Re;  //for classical electron radius
                                                   >> 215   G4double Ru;  //for classical muon radius
                                                   >> 216   G4double Eth; //for energy threshould to annihiToMu
                                                   >> 217 
                                                   >> 218   //parameters for ComputeCrossSection
                                                   >> 219   //
                                                   >> 220   G4double minBin = std::log10(fMinE/GeV);
                                                   >> 221   G4double maxBin = std::log10(fMaxE/GeV);
                                                   >> 222   de      = (maxBin - minBin)/G4double(fnBin);
                                                   >> 223   x       = minBin - de*0.5;
                                                   >> 224   atomicZ = 1.;
                                                   >> 225   atomicA = 2.;  
180                                                   226 
181   // get AnnihiToMuPair                        << 227   //set parameters for theory
182   //                                           << 
183   G4AnnihiToMuPair* annihiToMu =               << 
184     reinterpret_cast<G4AnnihiToMuPair*>(emCal. << 
185                                                << 
186   // parameters for ComputeCrossSection        << 
187   //                                           << 
188   G4double atomicZ = 1.;                       << 
189   G4double atomicA = 2.;                       << 
190                                                << 
191   // set parameters for theory                 << 
192   //                                              228   //
193   const G4ParticleDefinition* muon = G4MuonMin    229   const G4ParticleDefinition* muon = G4MuonMinus::MuonMinus();
194   G4double Mu = muon->GetPDGMass();            << 230   Mu   = muon->GetPDGMass();
195   G4double Me = electron_mass_c2;              << 231   Me   = electron_mass_c2;
196   G4double Re = classic_electr_radius;         << 232   Re   = classic_electr_radius;
197   G4double Ru = Re * Me / Mu;                  << 233   Ru   = Re*Me/Mu;
198   G4double Eth = 2 * Mu * Mu / Me - Me;        << 234   Eth  = 2*Mu*Mu/Me-Me;
199   G4PhysicsLogVector v(fMinE, fMaxE, fnBin, fa << 235   
200                                                << 236   //Compute CrossSections and Fill histgrams
201   // Compute CrossSections and Fill histgrams  << 237   //
202   //                                           << 238   for(int i=0;i<fnBin;i++){
203   for (G4int i = 0; i <= fnBin; ++i) {         << 239     x += de;
204     G4double energy = v.Energy(i);             << 240     energy=std::pow(10,x)*GeV;
205     G4double x = std::log10(energy / GeV);     << 
206                                                   241 
207     // CrossSectionPerAtom                     << 242     //CrossSectionPerAtom
208     //                                            243     //
209     G4double crs_annihiToMu = annihiToMu->Comp << 244     crs_annihiToMu = annihiToMu->ComputeCrossSectionPerAtom(energy,atomicZ);
210     // G4cout << "crs_annihiToMu(mkb)=" << crs << 245     fAnalysis->FillH1(7,x,crs_annihiToMu/microbarn);
211     fAnalysis->FillH1(7, x, crs_annihiToMu / m << 246   
212                                                << 247     crs_annihil = 
213     G4double crs_annihil =                     << 248       emCal.ComputeCrossSectionPerAtom(energy,positron,annihilName,
214       emCal.ComputeCrossSectionPerAtom(energy, << 249                                        atomicZ,atomicA);
215     fAnalysis->FillH1(8, x, crs_annihil / micr << 250     fAnalysis->FillH1(8,x,crs_annihil/microbarn);
216                                                << 251 
217     G4double crs_annihiToHadr =                << 252     crs_annihiToHadr = 
218       emCal.ComputeCrossSectionPerAtom(energy, << 253       emCal.ComputeCrossSectionPerAtom(energy,positron,annihiToHadrName,
219     fAnalysis->FillH1(9, x, crs_annihiToHadr / << 254                                        atomicZ,atomicA);
                                                   >> 255     fAnalysis->FillH1(9,x,crs_annihiToHadr/microbarn);
220                                                   256 
221     // CrossSectionPerVolume                   << 257     //CrossSectionPerVolume
222     //                                            258     //
223     G4double crsVol_Brem =                     << 259     crsVol_Brem = 
224       emCal.ComputeCrossSectionPerVolume(energ << 260       emCal.ComputeCrossSectionPerVolume(energy,positron,BremName,fMat,100*keV);
225     fAnalysis->FillH1(12, x, crsVol_Brem * mm) << 261     fAnalysis->FillH1(12,x,crsVol_Brem*mm);
226                                                << 262 
227     G4double crsVol_Ioni =                     << 263     crsVol_Ioni = 
228       emCal.ComputeCrossSectionPerVolume(energ << 264       emCal.ComputeCrossSectionPerVolume(energy,positron,IoniName,fMat,100*keV);
229     fAnalysis->FillH1(13, x, crsVol_Ioni * mm) << 265     fAnalysis->FillH1(13,x,crsVol_Ioni*mm);
230                                                << 266                 
231     G4double crsVol_annihiToMu = annihiToMu->C << 267     crsVol_annihiToMu = annihiToMu->CrossSectionPerVolume(energy,fMat);
232     fAnalysis->FillH1(14, x, crsVol_annihiToMu << 268     fAnalysis->FillH1(14,x,crsVol_annihiToMu*mm);
233                                                << 269   
234     G4double crsVol_annihil =                  << 270     crsVol_annihil = 
235       emCal.ComputeCrossSectionPerVolume(energ << 271       emCal.ComputeCrossSectionPerVolume(energy,positron,annihilName,fMat);
236     fAnalysis->FillH1(15, x, crsVol_annihil *  << 272     fAnalysis->FillH1(15,x,crsVol_annihil*mm);
237                                                << 273 
238     G4double crsVol_annihiToHadr =             << 274     crsVol_annihiToHadr = 
239       emCal.ComputeCrossSectionPerVolume(energ << 275       emCal.ComputeCrossSectionPerVolume(energy,positron,annihiToHadrName,fMat);
240     fAnalysis->FillH1(16, x, crsVol_annihiToHa << 276     fAnalysis->FillH1(16,x,crsVol_annihiToHadr*mm);
241                                                   277 
242     // R ratio                                 << 278     //R ratio
243     //                                            279     //
244     G4double RR = 0.0;                         << 280     RR = 0.0;
245     if (crsVol_annihiToMu > 0.) RR = crsVol_an << 281     if(crsVol_annihiToMu != 0) RR = crsVol_annihiToHadr/crsVol_annihiToMu;
246     fAnalysis->FillH1(17, x, RR);              << 282     fAnalysis->FillH1(17,x,RR);
247                                                   283 
248     // Theoretical calculation                 << 284     //Theoretical calculation
249     //                                            285     //
250     G4double X1 = energy / Me;                 << 286     X1 = energy/Me;
251     if (X1 > 1 && i % 1000 == 0) {             << 287     if(X1>1 && i%1000==0){
252       G4double crs_annihil_theory =            << 288       crs_annihil_theory = atomicZ*pi*Re*Re*
253         atomicZ * pi * Re * Re                 << 289         ( (X1*X1+4*X1+1)*G4Log(X1+std::sqrt(X1*X1-1))/(X1*X1-1)
254         * ((X1 * X1 + 4 * X1 + 1) * G4Log(X1 + << 290          -(X1+3)/std::sqrt(X1*X1-1) )/(X1+1);
255            - (X1 + 3) / std::sqrt(X1 * X1 - 1) << 291       fAnalysis->FillH1(10,x,crs_annihil_theory/microbarn);
256         / (X1 + 1);                            << 
257       fAnalysis->FillH1(10, x, crs_annihil_the << 
258     }                                             292     }
259                                                   293 
260     G4double X2 = Eth / energy;                << 294     X2 = Eth/energy;
261     if (X2 < 1. && i % 1000 == 0) {            << 295     if(X2<1 && i%1000==0){
262       G4double crs_annihiToMu_theory =         << 296       crs_annihiToMu_theory = atomicZ*pi*Ru*Ru/3*X2*(1+X2/2)*std::sqrt(1-X2);
263         atomicZ * pi * Ru * Ru / 3 * X2 * (1 + << 297       fAnalysis->FillH1(11,x,crs_annihiToMu_theory/microbarn);
264       fAnalysis->FillH1(11, x, crs_annihiToMu_ << 
265     }                                             298     }
266                                                   299 
267     // if(i%1000==0)G4cout <<"###energy:" << e << 300     //if(i%1000==0)G4cout <<"###energy:" << energy << "/crs_ToMuMu:"  
268     //         << crs_annihiToMu << "/crs_ToTw << 301     //        << crs_annihiToMu << "/crs_ToTwoGamma:"<< crs_annihil 
269     //         <<"/crs_ToToHadr:"<<crs_annihiT << 302     //        <<"/crs_ToToHadr:"<<crs_annihiToHadr<< G4endl;
270   }                                               303   }
271                                                << 304  
272   fAnalysis->Write();                             305   fAnalysis->Write();
273   fAnalysis->CloseFile();                         306   fAnalysis->CloseFile();
274   fAnalysis->Clear();                          << 
275                                                   307 
                                                   >> 308   delete fAnalysis;
276   G4cout << G4endl;                               309   G4cout << G4endl;
277 }                                                 310 }
278                                                   311 
279 //....oooOO0OOooo........oooOO0OOooo........oo    312 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
280                                                   313