ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RunAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RunAction.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 
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"
43 
44 #include "RunAction.hh"
45 #include "G4Run.hh"
46 #include "G4RunManager.hh"
47 #include "G4MuonMinus.hh"
48 
49 #include "Randomize.hh"
50 
51 #include <sstream>
52 
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 
55 RunAction::RunAction(DetectorConstruction* det)
56  : G4UserRunAction(),fDetector(det),fProcCounter(0),fAnalysis(0),fMat(0)
57 {
58  fMinE = 40*GeV;
59  fMaxE = 10000*GeV;
60  fnBin = 10000;
61  fNtColId[0] = fNtColId[1] = fNtColId[2] = 0;
62 }
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65 
66 RunAction::~RunAction()
67 {}
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 
71 void RunAction::BeginOfRunAction(const G4Run* aRun)
72 {
73  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
74 
75  //get material
76  //
77  fMat = fDetector->GetMaterial();
78  G4cout << "###RunAction::BeginOfRunAction Material:"
79  << fMat->GetName() << G4endl;
80 
81  fProcCounter = new ProcessesCount;
82 
83  fAnalysis = G4AnalysisManager::Instance();
84 
85  // Open an output file
86  //
87  std::stringstream tmp;
88  tmp << "testem6_" << aRun->GetRunID();
89  G4String fileName = tmp.str();
90  fAnalysis->OpenFile(fileName);
91  fAnalysis->SetVerboseLevel(2);
92  G4String extension = fAnalysis->GetFileType();
93  fileName = fileName + "." + extension;
94 
95  // Creating histograms 1,2,3,4,5,6
96  //
97  fAnalysis->SetFirstHistoId(1);
98  fAnalysis->CreateH1("h1","1/(1+(theta+[g]+)**2)",100, 0 ,1.);
99  fAnalysis->CreateH1("h2","log10(theta+ [g]+)", 100,-3.,1.);
100  fAnalysis->CreateH1("h3","log10(theta- [g]-)", 100,-3.,1.);
101  fAnalysis->CreateH1("h4","log10(theta+ [g]+ -theta- [g]-)", 100,-3.,1.);
102  fAnalysis->CreateH1("h5","xPlus" ,100,0.,1.);
103  fAnalysis->CreateH1("h6","xMinus",100,0.,1.);
104 
105  //creating histogram 7,8,9,10,11 (CrossSectionPerAtom)
106  //
107  G4double minBin = std::log10(fMinE/GeV);
108  G4double maxBin = std::log10(fMaxE/GeV);
109  fAnalysis->CreateH1("h7","CrossSectionPerAtom of AnnihiToMuMu (microbarn)",
110  fnBin,minBin,maxBin);
111  fAnalysis->CreateH1("h8",
112  "CrossSectionPerAtom of AnnihiToTwoGamma (microbarn)",fnBin,minBin,maxBin);
113  fAnalysis->CreateH1("h9","CrossSectionPerAtom of AnnihiToHadrons (microbarn)",
114  fnBin,minBin,maxBin);
115  fAnalysis->CreateH1("h10",
116  "Theoretical CrossSectionPerAtom of AnnihiToTwoGamma (microbarn)",
117  fnBin,minBin,maxBin);
118  fAnalysis->CreateH1("h11",
119  "Theoretical CrossSectionPerAtom of AnnihiToMuMu (microbarn)",
120  fnBin,minBin,maxBin);
121 
122  //creating histogram 12,13,14,15,16(CrossSectionPerVolume)
123  //
124  fAnalysis->CreateH1("h12","CrossSectionPerVol of Bremsstraulung (1/mm) ",
125  fnBin,minBin,maxBin);
126  fAnalysis->CreateH1("h13","CrossSectionPerVol of Ionization (1/mm)",
127  fnBin,minBin,maxBin);
128  fAnalysis->CreateH1("h14","CrossSectionPerVol of AnnihiToMuMu (1/mm)",
129  fnBin,minBin,maxBin);
130  fAnalysis->CreateH1("h15","CrossSectionPerVol of AnnihiToTwoGamma (1/mm)",
131  fnBin,minBin,maxBin);
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;
139 }
140 
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
142 
143 void RunAction::CountProcesses(G4String procName)
144 {
145  //does the process already encounted ?
146  //
147  size_t nbProc = fProcCounter->size();
148  size_t i = 0;
149  while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++;
150  if (i == nbProc) fProcCounter->push_back( new OneProcessCount(procName));
151 
152  (*fProcCounter)[i]->Count();
153 }
154 
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
156 
157 void RunAction::EndOfRunAction(const G4Run*)
158 {
159  // show Rndm status : not needed
160  //
161  //CLHEP::HepRandom::showEngineStatus();
162 
163  //total number of process calls
164  //
165  G4double countTot = 0.;
166  G4cout << "\n Number of process calls --->";
167  for (size_t i=0; i< fProcCounter->size();i++) {
168  G4String procName = (*fProcCounter)[i]->GetName();
169  if (procName != "Transportation") {
170  G4int count = (*fProcCounter)[i]->GetCounter();
171  G4cout << "\t" << procName << " : " << count;
172  countTot += count;
173  }
174  }
175 
176  //instance EmCalculator
177  //
178  G4EmCalculator emCal;
179  emCal.SetVerbose(0);
180 
181  //get positron
182  //
183  G4String positronName = "e+";
184  G4ParticleDefinition* positron =
185  G4ParticleTable::GetParticleTable()->FindParticle(positronName);
186 
187  //process name
188  //
189  G4String annihilName = "annihil";
190  G4String annihiToMuName = "AnnihiToMuPair";
191  G4String annihiToHadrName = "ee2hadr";
192  G4String BremName = "eBrem";
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.;
225 
226  //set parameters for theory
227  //
228  const G4ParticleDefinition* muon = G4MuonMinus::MuonMinus();
229  Mu = muon->GetPDGMass();
230  Me = electron_mass_c2;
231  Re = classic_electr_radius;
232  Ru = Re*Me/Mu;
233  Eth = 2*Mu*Mu/Me-Me;
234 
235  //Compute CrossSections and Fill histgrams
236  //
237  for(int i=0;i<fnBin;i++){
238  x += de;
239  energy=std::pow(10,x)*GeV;
240 
241  //CrossSectionPerAtom
242  //
243  crs_annihiToMu = annihiToMu->ComputeCrossSectionPerAtom(energy,atomicZ);
244  fAnalysis->FillH1(7,x,crs_annihiToMu/microbarn);
245 
246  crs_annihil =
247  emCal.ComputeCrossSectionPerAtom(energy,positron,annihilName,
248  atomicZ,atomicA);
249  fAnalysis->FillH1(8,x,crs_annihil/microbarn);
250 
251  crs_annihiToHadr =
252  emCal.ComputeCrossSectionPerAtom(energy,positron,annihiToHadrName,
253  atomicZ,atomicA);
254  fAnalysis->FillH1(9,x,crs_annihiToHadr/microbarn);
255 
256  //CrossSectionPerVolume
257  //
258  crsVol_Brem =
259  emCal.ComputeCrossSectionPerVolume(energy,positron,BremName,fMat,100*keV);
260  fAnalysis->FillH1(12,x,crsVol_Brem*mm);
261 
262  crsVol_Ioni =
263  emCal.ComputeCrossSectionPerVolume(energy,positron,IoniName,fMat,100*keV);
264  fAnalysis->FillH1(13,x,crsVol_Ioni*mm);
265 
266  crsVol_annihiToMu = annihiToMu->CrossSectionPerVolume(energy,fMat);
267  fAnalysis->FillH1(14,x,crsVol_annihiToMu*mm);
268 
269  crsVol_annihil =
270  emCal.ComputeCrossSectionPerVolume(energy,positron,annihilName,fMat);
271  fAnalysis->FillH1(15,x,crsVol_annihil*mm);
272 
273  crsVol_annihiToHadr =
274  emCal.ComputeCrossSectionPerVolume(energy,positron,annihiToHadrName,fMat);
275  fAnalysis->FillH1(16,x,crsVol_annihiToHadr*mm);
276 
277  //R ratio
278  //
279  RR = 0.0;
280  if(crsVol_annihiToMu != 0) RR = crsVol_annihiToHadr/crsVol_annihiToMu;
281  fAnalysis->FillH1(17,x,RR);
282 
283  //Theoretical calculation
284  //
285  X1 = energy/Me;
286  if(X1>1 && i%1000==0){
287  crs_annihil_theory = atomicZ*pi*Re*Re*
288  ( (X1*X1+4*X1+1)*G4Log(X1+std::sqrt(X1*X1-1))/(X1*X1-1)
289  -(X1+3)/std::sqrt(X1*X1-1) )/(X1+1);
290  fAnalysis->FillH1(10,x,crs_annihil_theory/microbarn);
291  }
292 
293  X2 = Eth/energy;
294  if(X2<1 && i%1000==0){
295  crs_annihiToMu_theory = atomicZ*pi*Ru*Ru/3*X2*(1+X2/2)*std::sqrt(1-X2);
296  fAnalysis->FillH1(11,x,crs_annihiToMu_theory/microbarn);
297  }
298 
299  //if(i%1000==0)G4cout <<"###energy:" << energy << "/crs_ToMuMu:"
300  // << crs_annihiToMu << "/crs_ToTwoGamma:"<< crs_annihil
301  // <<"/crs_ToToHadr:"<<crs_annihiToHadr<< G4endl;
302  }
303 
304  fAnalysis->Write();
305  fAnalysis->CloseFile();
306 
307  delete fAnalysis;
308  G4cout << G4endl;
309 }
310 
311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......