ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Run.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Run.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 "Run.hh"
34 #include "DetectorConstruction.hh"
35 #include "PrimaryGeneratorAction.hh"
36 #include "HistoManager.hh"
37 #include "G4Track.hh"
38 #include "G4VPhysicalVolume.hh"
39 
40 #include "G4EmCalculator.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4UnitsTable.hh"
43 
44 #include <iomanip>
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 
48 Run::Run(DetectorConstruction* det, PrimaryGeneratorAction* /*prim*/,
49  HistoManager* histoMgr)
50 : G4Run(),
51  fDetector(det),
52  fHistoMgr(histoMgr)
53 {
54 
55  fSumR=0;
56  fNevt = fNelec = fNposit = fNgam = fNstep = fNgamPh =
57  fNgamTar = fNeTar = fNePh = fNstepTarget = 0;
58 
59  fAnalysisManager = G4AnalysisManager::Instance();
60 
61  fHistoId = fHistoMgr->GetHistoIdentifiers();
62  fNHisto = fHistoId.size();
63 
64  fCheckVolume = fDetector->GetCheckVolume();
65  fGasVolume = fDetector->GetGasVolume();
66  fPhantom = fDetector->GetPhantom();
67  fTarget1 = fDetector->GetTarget1();
68  fTarget2 = fDetector->GetTarget2();
69 
70  fNBinsR = fDetector->GetNumberDivR();
71  fNBinsZ = fDetector->GetNumberDivZ();
72 
73  fScoreZ = fDetector->GetScoreZ();
74  fAbsorberZ = fDetector->GetAbsorberZ();
75  fAbsorberR = fDetector->GetAbsorberR();
76  fMaxEnergy = fDetector->GetMaxEnergy();
77 
78  fNBinsE = fDetector->GetNumberDivE();
79  fMaxEnergy = fDetector->GetMaxEnergy();
80 
81  fStepZ = fAbsorberZ/(G4double)fNBinsZ;
82  fStepR = fAbsorberR/(G4double)fNBinsR;
83  fStepE = fMaxEnergy/(G4double)fNBinsE;
84  fScoreBin = (G4int)(fScoreZ/fStepZ + 0.5);
85 
86  fVerbose = fDetector->GetVerbose();
87 
88  fGamma = G4Gamma::Gamma();
89  fElectron = G4Electron::Electron();
90  fPositron = G4Positron::Positron();
91 
92  G4cout << " "<< fNBinsR << " bins R stepR= " << fStepR/mm << " mm "
93  << G4endl;
94  G4cout << " "<< fNBinsZ << " bins Z stepZ= " << fStepZ/mm << " mm "
95  << G4endl;
96  G4cout << " "<< fNBinsE << " bins E stepE= " << fStepE/MeV << " MeV "
97  << G4endl;
98  G4cout << " "<< fScoreBin << "-th bin in Z is used for R distribution"
99  << G4endl;
100 
101  fVolumeR.clear();
102  fEdep.clear();
103  fGammaE.clear();
104 
105  fVolumeR.resize(fNBinsR,0.0);
106  fEdep.resize(fNBinsR, 0.0);
107  fGammaE.resize(fNBinsE, 0.0);
108 
109  G4double r1 = 0.0;
110  G4double r2 = fStepR;
111  for(G4int i=0; i<fNBinsR; ++i) {
112  fVolumeR[i] = cm*cm/(CLHEP::pi*(r2*r2 - r1*r1));
113  r1 = r2;
114  r2 += fStepR;
115  }
116 
117  if(fAnalysisManager->GetFileName()!="")fHistoMgr->Update(det, true);
118 
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 
123 Run::~Run()
124 {
125 
126 }
127 
128 
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130 
131 void Run::Merge(const G4Run* run)
132 {
133  const Run* localRun = static_cast<const Run*>(run);
134 
135  fNevt += localRun->fNevt;
136 
137  fNelec += localRun->fNelec;
138  fNgam += localRun->fNgam;
139  fNposit += localRun->fNposit;
140 
141  fNgamTar+= localRun->fNgamTar;
142  fNgamPh += localRun->fNgamPh;
143  fNeTar += localRun->fNeTar;
144  fNePh += localRun->fNePh;
145 
146  fNstep += localRun->fNstep;
147  fNstepTarget += localRun->fNstepTarget;
148 
149  fSumR += localRun->fSumR;
150 
151  for(int i=0; i<(int)fEdep.size(); i++) fEdep[i] += localRun->fEdep[i];
152  for(int i=0; i<(int)fGammaE.size(); i++) fGammaE[i] += localRun->fGammaE[i];
153 
154  G4Run::Merge(run);
155 }
156 
157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
158 
159 void Run::EndOfRun()
160 {
161  G4cout << "Histo: End of run actions are started" << G4endl;
162 
163  // average
164  G4cout<<"========================================================"<<G4endl;
165  G4double x = (G4double)fNevt;
166  if(fNevt > 0) { x = 1.0/x; }
167  G4double xe = x*(G4double)fNelec;
168  G4double xg = x*(G4double)fNgam;
169  G4double xp = x*(G4double)fNposit;
170  G4double xs = x*(G4double)fNstep;
171  G4double xph= x*(G4double)fNgamPh;
172  G4double xes= x*(G4double)fNstepTarget;
173  G4double xgt= x*(G4double)fNgamTar;
174  G4double xet= x*(G4double)fNeTar;
175  G4double xphe= x*(G4double)fNePh;
176 
177  G4cout << "Number of events "
178  << std::setprecision(8) << fNevt <<G4endl;
179  G4cout
180  << std::setprecision(4) << "Average number of e- "
181  << xe << G4endl;
182  G4cout
183  << std::setprecision(4) << "Average number of gamma "
184  << xg << G4endl;
185  G4cout
186  << std::setprecision(4) << "Average number of e+ "
187  << xp << G4endl;
188  G4cout
189  << std::setprecision(4) << "Average number of steps in the phantom "
190  << xs << G4endl;
191  G4cout
192  << std::setprecision(4) << "Average number of e- steps in the target "
193  << xes << G4endl;
194  G4cout
195  << std::setprecision(4) << "Average number of g produced in the target "
196  << xgt << G4endl;
197  G4cout
198  << std::setprecision(4) << "Average number of e- produced in the target "
199  << xet << G4endl;
200  G4cout
201  << std::setprecision(4) << "Average number of g produced in the phantom "
202  << xph << G4endl;
203  G4cout
204  << std::setprecision(4) << "Average number of e- produced in the phantom "
205  << xphe << G4endl;
206  G4cout
207  << std::setprecision(4) << "Total fGamma fluence in front of the phantom "
208  << x*fSumR/MeV << " MeV " << G4endl;
209  G4cout<<"========================================================"<<G4endl;
210  G4cout<<G4endl;
211  G4cout<<G4endl;
212 
213  G4double sum = 0.0;
214  for(G4int i=0; i<fNBinsR; i++) {
215  fEdep[i] *= x;
216  sum += fEdep[i];
217  }
218 
219 
220  if(fAnalysisManager) {
221 
222  if(fAnalysisManager->IsActive()) {
223 
224  // normalise histograms
225  for(G4int i=0; i<fNHisto; i++) {
226  fAnalysisManager->GetH1(fHistoId[i])->scale(x);
227  }
228  G4double nr = fEdep[0];
229  if(nr > 0.0) { nr = 1./nr; }
230  fAnalysisManager->GetH1(fHistoId[0])->scale(nr);
231 
232  nr = sum*fStepR;
233  if(nr > 0.0) { nr = 1./nr; }
234  fAnalysisManager->GetH1(fHistoId[1])->scale(nr);
235 
236  fAnalysisManager->GetH1(fHistoId[3])
237  ->scale(1000.0*cm3/(CLHEP::pi*fAbsorberR*fAbsorberR*fStepZ));
238  fAnalysisManager->GetH1(fHistoId[4])
239  ->scale(1000.0*cm3*fVolumeR[0]/fStepZ);
240 
241  // Write histogram file
242  if(!fAnalysisManager->Write()) {
243  G4Exception ("Histo::Save()", "hist01", FatalException,
244  "Cannot write ROOT file.");
245  }
246  G4cout << "### Histo::Save: Histograms are saved" << G4endl;
247  if(fAnalysisManager->CloseFile() && fVerbose) {
248  G4cout << " File is closed" << G4endl;
249  }
250  }
251 
252  delete fAnalysisManager;
253  fAnalysisManager = 0;
254  }
255 
256 
257 }
258 
259 
260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
261 
262 void Run::AddTargetPhoton(const G4DynamicParticle* p)
263 {
264  ++fNgamTar;
265  if(fAnalysisManager) {
266  fAnalysisManager->FillH1(fHistoId[5],p->GetKineticEnergy()/MeV,1.0);
267  }
268 }
269 
270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
271 
272 void Run::AddPhantomPhoton(const G4DynamicParticle*)
273 {
274  ++fNgamPh;
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
278 
279 void Run::AddTargetElectron(const G4DynamicParticle* p)
280 {
281  ++fNeTar;
282  if(fAnalysisManager) {
283  fAnalysisManager->FillH1(fHistoId[8],p->GetKineticEnergy()/MeV,1.0);
284  }
285 }
286 
287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 
289 void Run::AddPhantomElectron(const G4DynamicParticle* p)
290 {
291  ++fNePh;
292  if(fAnalysisManager) {
293  fAnalysisManager->FillH1(fHistoId[7],p->GetKineticEnergy()/MeV,1.0);
294  }
295 }
296 
297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
298 
299 void Run::ScoreNewTrack(const G4Track* aTrack)
300 {
301 
302  //Save primary parameters
303  const G4ParticleDefinition* particle = aTrack->GetParticleDefinition();
304  G4int pid = aTrack->GetParentID();
305  G4VPhysicalVolume* pv = aTrack->GetVolume();
306  const G4DynamicParticle* dp = aTrack->GetDynamicParticle();
307 
308  //primary particle
309  if(0 == pid) {
310 
311  ++fNevt;
312  if(fVerbose)
313  {
314  G4ThreeVector pos = aTrack->GetVertexPosition();
315  G4ThreeVector dir = aTrack->GetMomentumDirection();
316  G4cout << "TrackingAction: Primary "
317  << particle->GetParticleName()
318  << " Ekin(MeV)= "
319  << aTrack->GetKineticEnergy()/MeV
320  << "; pos= " << pos << "; dir= " << dir << G4endl;
321  }
322  // secondary electron
323  }
324  else if (0 < pid && particle == fElectron)
325  {
326  if(fVerbose) {
327  G4cout << "TrackingAction: Secondary electron " << G4endl;
328  }
329  AddElectron();
330  if(pv == fPhantom) { AddPhantomElectron(dp); }
331  else if(pv == fTarget1 || pv == fTarget2) { AddTargetElectron(dp); }
332 
333  // secondary positron
334  }
335  else if (0 < pid && particle == fPositron) {
336  if(fVerbose) {
337  G4cout << "TrackingAction: Secondary positron " << G4endl;
338  }
339  AddPositron();
340 
341  // secondary gamma
342  }
343  else if (0 < pid && particle == fGamma) {
344  if(fVerbose) {
345  G4cout << "TrackingAction: Secondary gamma; parentID= " << pid
346  << " E= " << aTrack->GetKineticEnergy() << G4endl;
347  }
348  AddPhoton();
349  if(pv == fPhantom) { AddPhantomPhoton(dp); }
350  else if(pv == fTarget1 || pv == fTarget2) { AddTargetPhoton(dp); }
351 
352  }
353 }
354 
355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
356 
357 void Run::AddPhantomGamma(G4double e, G4double r)
358 {
359  e /= MeV;
360  fSumR += e;
361  G4int bin = (G4int)(e/fStepE);
362  if(bin >= fNBinsE) { bin = fNBinsE-1; }
363  fGammaE[bin] += e;
364  G4int bin1 = (G4int)(r/fStepR);
365  if(bin1 >= fNBinsR) { bin1 = fNBinsR-1; }
366  if(fAnalysisManager) {
367  G4AnalysisManager::Instance()->FillH1(fHistoId[6],e,1.0);
368  G4AnalysisManager::Instance()->FillH1(fHistoId[9],r,e*fVolumeR[bin1]);
369  }
370 }
371 
372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
373 
374 void Run::AddPhantomStep(G4double edep, G4double r1, G4double z1,
375  G4double r2, G4double z2,
376  G4double r0, G4double z0)
377 {
378  ++fNstep;
379  G4int nzbin = (G4int)(z0/fStepZ);
380  if(fVerbose) {
381  G4cout << "Histo: edep(MeV)= " << edep/MeV << " at binz= " << nzbin
382  << " r1= " << r1 << " z1= " << z1
383  << " r2= " << r2 << " z2= " << z2
384  << " r0= " << r0 << " z0= " << z0
385  << G4endl;
386  }
387  if(nzbin == fScoreBin) {
388  G4int bin = (G4int)(r0/fStepR);
389  if(bin >= fNBinsR) { bin = fNBinsR-1; }
390  double w = edep*fVolumeR[bin];
391  fEdep[bin] += w;
392 
393  if(fAnalysisManager) {
394  G4AnalysisManager::Instance()->FillH1(fHistoId[0],r0,w);
395  G4AnalysisManager::Instance()->FillH1(fHistoId[1],r0,w);
396  G4AnalysisManager::Instance()->FillH1(fHistoId[2],r0,w);
397  }
398  }
399  G4int bin1 = (G4int)(z1/fStepZ);
400  if(bin1 >= fNBinsZ) { bin1 = fNBinsZ-1; }
401  G4int bin2 = (G4int)(z2/fStepZ);
402  if(bin2 >= fNBinsZ) { bin2 = fNBinsZ-1; }
403  if(bin1 == bin2) {
404  if(fAnalysisManager) {
405  G4AnalysisManager::Instance()->FillH1(fHistoId[3],z0,edep);
406  if(r1 < fStepR) {
407  G4double w = edep;
408  if(r2 > fStepR) { w *= (fStepR - r1)/(r2 - r1); }
409  G4AnalysisManager::Instance()->FillH1(fHistoId[4],z0,w);
410  }
411  }
412  } else {
413  G4int bin;
414 
415  if(bin2 < bin1) {
416  bin = bin2;
417  G4double z = z2;
418  bin2 = bin1;
419  z2 = z1;
420  bin1 = bin;
421  z1 = z;
422  }
423  G4double zz1 = z1;
424  G4double zz2 = (bin1+1)*fStepZ;
425  G4double rr1 = r1;
426  G4double dz = z2 - z1;
427  G4double dr = r2 - r1;
428  G4double rr2 = r1 + dr*(zz2-zz1)/dz;
429  for(bin=bin1; bin<=bin2; bin++) {
430  if(fAnalysisManager) {
431  G4double de = edep*(zz2 - zz1)/dz;
432  G4double zf = (zz1+zz2)*0.5;
433  { G4AnalysisManager::Instance()->FillH1(fHistoId[3],zf,de); }
434  if(rr1 < fStepR) {
435  G4double w = de;
436  if(rr2 > fStepR) w *= (fStepR - rr1)/(rr2 - rr1);
437  { G4AnalysisManager::Instance()->FillH1(fHistoId[4],zf,w); }
438  }
439  }
440  zz1 = zz2;
441  zz2 = std::min(z2, zz1+fStepZ);
442  rr1 = rr2;
443  rr2 = rr1 + dr*(zz2 - zz1)/dz;
444  }
445  }
446 }
447 
448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......