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.7)


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