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 
31 #include "Run.hh"
32 #include "DetectorConstruction.hh"
33 #include "PrimaryGeneratorAction.hh"
34 #include "HistoManager.hh"
35 
36 #include "G4ProcessTable.hh"
37 #include "G4Radioactivation.hh"
38 #include "G4TwoVector.hh"
39 #include "G4UnitsTable.hh"
40 #include "G4SystemOfUnits.hh"
41 
42 #include <fstream>
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 
46 Run::Run(DetectorConstruction* det)
47 : G4Run(),
48  fDetector(det), fParticle(0), fEkin(0.)
49 {
50  fEdepTarget = fEdepTarget2 = 0.;
51  fEdepDetect = fEdepDetect2 = 0.;
52 }
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 
56 Run::~Run()
57 { }
58 
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 
61 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
62 {
63  fParticle = particle;
64  fEkin = energy;
65 }
66 
67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68 
69 void Run::CountProcesses(const G4VProcess* process, G4int iVol)
70 {
71  G4String procName = process->GetProcessName();
72 
73  if (iVol == 1) {
74  std::map<G4String,G4int>::iterator it1 = fProcCounter1.find(procName);
75  if ( it1 == fProcCounter1.end()) {
76  fProcCounter1[procName] = 1;
77  }
78  else {
79  fProcCounter1[procName]++;
80  }
81  }
82 
83  if (iVol == 2) {
84  std::map<G4String,G4int>::iterator it2 = fProcCounter2.find(procName);
85  if ( it2 == fProcCounter2.end()) {
86  fProcCounter2[procName] = 1;
87  }
88  else {
89  fProcCounter2[procName]++;
90  }
91  }
92 }
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95 
96 void Run::ParticleCount(G4String name, G4double Ekin, G4int iVol)
97 {
98  if (iVol == 1) {
99  std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
100  if ( it == fParticleDataMap1.end()) {
101  fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin);
102  }
103  else {
104  ParticleData& data = it->second;
105  data.fCount++;
106  data.fEmean += Ekin;
107  //update min max
108  G4double emin = data.fEmin;
109  if (Ekin < emin) data.fEmin = Ekin;
110  G4double emax = data.fEmax;
111  if (Ekin > emax) data.fEmax = Ekin;
112  }
113  }
114 
115  if (iVol == 2) {
116  std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
117  if ( it == fParticleDataMap2.end()) {
118  fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin);
119  }
120  else {
121  ParticleData& data = it->second;
122  data.fCount++;
123  data.fEmean += Ekin;
124  //update min max
125  G4double emin = data.fEmin;
126  if (Ekin < emin) data.fEmin = Ekin;
127  G4double emax = data.fEmax;
128  if (Ekin > emax) data.fEmax = Ekin;
129  }
130  }
131 }
132 
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134 
135 void Run::AddEdep(G4double edep1, G4double edep2)
136 {
137  fEdepTarget += edep1;
138  fEdepTarget2 += edep1*edep1;
139  fEdepDetect += edep2;
140  fEdepDetect2 += edep2*edep2;
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 
145 void Run::Merge(const G4Run* run)
146 {
147  const Run* localRun = static_cast<const Run*>(run);
148 
149  //primary particle info
150  //
151  fParticle = localRun->fParticle;
152  fEkin = localRun->fEkin;
153 
154  // accumulate sums
155  //
156  fEdepTarget += localRun->fEdepTarget;
157  fEdepTarget2 += localRun->fEdepTarget2;
158  fEdepDetect += localRun->fEdepDetect;
159  fEdepDetect2 += localRun->fEdepDetect2;
160 
161  //map: processes count in target
162 
163  std::map<G4String,G4int>::const_iterator itp1;
164  for ( itp1 = localRun->fProcCounter1.begin();
165  itp1 != localRun->fProcCounter1.end(); ++itp1 ) {
166 
167  G4String procName = itp1->first;
168  G4int localCount = itp1->second;
169  if ( fProcCounter1.find(procName) == fProcCounter1.end()) {
170  fProcCounter1[procName] = localCount;
171  }
172  else {
173  fProcCounter1[procName] += localCount;
174  }
175  }
176 
177  //map: processes count in detector
178 
179  std::map<G4String,G4int>::const_iterator itp2;
180  for ( itp2 = localRun->fProcCounter2.begin();
181  itp2 != localRun->fProcCounter2.end(); ++itp2 ) {
182 
183  G4String procName = itp2->first;
184  G4int localCount = itp2->second;
185  if ( fProcCounter2.find(procName) == fProcCounter2.end()) {
186  fProcCounter2[procName] = localCount;
187  }
188  else {
189  fProcCounter2[procName] += localCount;
190  }
191  }
192 
193  //map: created particles in target
194  std::map<G4String,ParticleData>::const_iterator itc;
195  for (itc = localRun->fParticleDataMap1.begin();
196  itc != localRun->fParticleDataMap1.end(); ++itc) {
197 
198  G4String name = itc->first;
199  const ParticleData& localData = itc->second;
200  if ( fParticleDataMap1.find(name) == fParticleDataMap1.end()) {
201  fParticleDataMap1[name]
202  = ParticleData(localData.fCount,
203  localData.fEmean,
204  localData.fEmin,
205  localData.fEmax);
206  }
207  else {
208  ParticleData& data = fParticleDataMap1[name];
209  data.fCount += localData.fCount;
210  data.fEmean += localData.fEmean;
211  G4double emin = localData.fEmin;
212  if (emin < data.fEmin) data.fEmin = emin;
213  G4double emax = localData.fEmax;
214  if (emax > data.fEmax) data.fEmax = emax;
215  }
216  }
217 
218  //map: created particle in detector
219  std::map<G4String,ParticleData>::const_iterator itn;
220  for (itn = localRun->fParticleDataMap2.begin();
221  itn != localRun->fParticleDataMap2.end(); ++itn) {
222 
223  G4String name = itn->first;
224  const ParticleData& localData = itn->second;
225  if ( fParticleDataMap2.find(name) == fParticleDataMap2.end()) {
226  fParticleDataMap2[name]
227  = ParticleData(localData.fCount,
228  localData.fEmean,
229  localData.fEmin,
230  localData.fEmax);
231  }
232  else {
233  ParticleData& data = fParticleDataMap2[name];
234  data.fCount += localData.fCount;
235  data.fEmean += localData.fEmean;
236  G4double emin = localData.fEmin;
237  if (emin < data.fEmin) data.fEmin = emin;
238  G4double emax = localData.fEmax;
239  if (emax > data.fEmax) data.fEmax = emax;
240  }
241  }
242 
243  G4Run::Merge(run);
244 }
245 
246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
247 
248 void Run::EndOfRun()
249 {
250  G4int prec = 5, wid = prec + 2;
251  G4int dfprec = G4cout.precision(prec);
252 
253  // run condition
254  //
255  G4String Particle = fParticle->GetParticleName();
256  G4cout << "\n The run is " << numberOfEvent << " "<< Particle << " of "
257  << G4BestUnit(fEkin,"Energy") << " through : ";
258 
259  G4cout << "\n Target : Length = "
260  << G4BestUnit(fDetector->GetTargetLength(),"Length")
261  << " Radius = "
262  << G4BestUnit(fDetector->GetTargetRadius(),"Length")
263  << " Material = "
264  << fDetector->GetTargetMaterial()->GetName();
265  G4cout << "\n Detector : Length = "
266  << G4BestUnit(fDetector->GetDetectorLength(),"Length")
267  << " Thickness = "
268  << G4BestUnit(fDetector->GetDetectorThickness(),"Length")
269  << " Material = "
270  << fDetector->GetDetectorMaterial()->GetName() << G4endl;
271 
272  if (numberOfEvent == 0) { G4cout.precision(dfprec); return;}
273 
274  // compute mean Energy deposited and rms in target
275  //
276  G4int TotNbofEvents = numberOfEvent;
277  fEdepTarget /= TotNbofEvents; fEdepTarget2 /= TotNbofEvents;
278  G4double rmsEdep = fEdepTarget2 - fEdepTarget*fEdepTarget;
279  if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep);
280  else rmsEdep = 0.;
281 
282  G4cout << "\n Mean energy deposit in target, in time window = "
283  << G4BestUnit(fEdepTarget,"Energy") << "; rms = "
284  << G4BestUnit(rmsEdep, "Energy")
285  << G4endl;
286 
287  // compute mean Energy deposited and rms in detector
288  //
289  fEdepDetect /= TotNbofEvents; fEdepDetect2 /= TotNbofEvents;
290  rmsEdep = fEdepDetect2 - fEdepDetect*fEdepDetect;
291  if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep);
292  else rmsEdep = 0.;
293 
294  G4cout << " Mean energy deposit in detector, in time window = "
295  << G4BestUnit(fEdepDetect,"Energy") << "; rms = "
296  << G4BestUnit(rmsEdep, "Energy")
297  << G4endl;
298 
299  // frequency of processes in target
300  //
301  G4cout << "\n Process calls frequency in target :" << G4endl;
302  G4int index = 0;
303  std::map<G4String,G4int>::iterator it1;
304  for (it1 = fProcCounter1.begin(); it1 != fProcCounter1.end(); it1++) {
305  G4String procName = it1->first;
306  G4int count = it1->second;
307  G4String space = " "; if (++index%3 == 0) space = "\n";
308  G4cout << " " << std::setw(20) << procName << "="<< std::setw(7) << count
309  << space;
310  }
311  G4cout << G4endl;
312 
313  // frequency of processes in detector
314  //
315  G4cout << "\n Process calls frequency in detector:" << G4endl;
316  index = 0;
317  std::map<G4String,G4int>::iterator it2;
318  for (it2 = fProcCounter2.begin(); it2 != fProcCounter2.end(); it2++) {
319  G4String procName = it2->first;
320  G4int count = it2->second;
321  G4String space = " "; if (++index%3 == 0) space = "\n";
322  G4cout << " " << std::setw(20) << procName << "="<< std::setw(7) << count
323  << space;
324  }
325  G4cout << G4endl;
326 
327  // particles count in target
328  //
329  G4cout << "\n List of generated particles in target:" << G4endl;
330 
331  std::map<G4String,ParticleData>::iterator itc;
332  for (itc = fParticleDataMap1.begin(); itc != fParticleDataMap1.end(); itc++) {
333  G4String name = itc->first;
334  ParticleData data = itc->second;
335  G4int count = data.fCount;
336  G4double eMean = data.fEmean/count;
337  G4double eMin = data.fEmin;
338  G4double eMax = data.fEmax;
339 
340  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
341  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
342  << "\t( " << G4BestUnit(eMin, "Energy")
343  << " --> " << G4BestUnit(eMax, "Energy")
344  << ")" << G4endl;
345  }
346 
347  // particles count in detector
348  //
349  G4cout << "\n List of generated particles in detector:" << G4endl;
350 
351  std::map<G4String,ParticleData>::iterator itn;
352  for (itn = fParticleDataMap2.begin(); itn != fParticleDataMap2.end(); itn++) {
353  G4String name = itn->first;
354  ParticleData data = itn->second;
355  G4int count = data.fCount;
356  G4double eMean = data.fEmean/count;
357  G4double eMin = data.fEmin;
358  G4double eMax = data.fEmax;
359 
360  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
361  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
362  << "\t( " << G4BestUnit(eMin, "Energy")
363  << " --> " << G4BestUnit(eMax, "Energy") << ")" << G4endl;
364  }
365  G4cout << G4endl;
366 
367  // activities in VR mode
368  //
369  WriteActivity(numberOfEvent);
370 
371  //remove all contents in fProcCounter, fCount
372  fProcCounter1.clear();
373  fProcCounter2.clear();
374  fParticleDataMap1.clear();
375  fParticleDataMap2.clear();
376 
377  //restore default format
378  G4cout.precision(dfprec);
379 }
380 
381 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
382 
383 void Run::WriteActivity(G4int nevent)
384 {
385  G4ProcessTable *pTable = G4ProcessTable::GetProcessTable();
386  G4Radioactivation* rDecay = (G4Radioactivation *)
387  pTable->FindProcess("Radioactivation", "GenericIon");
388 
389  // output the induced radioactivities (in VR mode only)
390  //
391  if ((rDecay == 0) || (rDecay->IsAnalogueMonteCarlo())) return;
392 
393  G4String fileName = G4AnalysisManager::Instance()->GetFileName() + ".activity";
394  std::ofstream outfile (fileName, std::ios::out );
395 
396  std::vector<G4RadioactivityTable*> theTables =
397  rDecay->GetTheRadioactivityTables();
398 
399  for (size_t i = 0 ; i < theTables.size(); i++) {
400  G4double rate, error;
401  outfile << "Radioactivities in decay window no. " << i << G4endl;
402  outfile << "Z \tA \tE \tActivity (decays/window) \tError (decays/window) "
403  << G4endl;
404 
405  map<G4ThreeVector,G4TwoVector> *aMap = theTables[i]->GetTheMap();
406  map<G4ThreeVector,G4TwoVector>::iterator iter;
407  for (iter=aMap->begin(); iter != aMap->end(); iter++) {
408  rate = iter->second.x()/nevent;
409  error = std::sqrt(iter->second.y())/nevent;
410  if (rate < 0.) rate = 0.; // statically it can be < 0.
411  outfile << iter->first.x() <<"\t"<< iter->first.y() <<"\t"
412  << iter->first.z() << "\t" << rate <<"\t" << error << G4endl;
413  }
414  outfile << G4endl;
415  }
416  outfile.close();
417 }
418 
419 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....