Geant4  10.00.p01
exrdmAnalysisManager.cc
Go to the documentation of this file.
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 //
26 // $Id: exrdmAnalysisManager.cc 79266 2014-02-20 16:41:56Z gcosmo $
27 //
30 //
31 
32 
33 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35 
36 #include "exrdmAnalysisManager.hh"
37 #include "G4UnitsTable.hh"
38 #include "exrdmHisto.hh"
39 #include "G4ProcessTable.hh"
40 #include "G4RadioactiveDecay.hh"
41 #include "G4TwoVector.hh"
42 #include "G4SystemOfUnits.hh"
43 
44 #include <fstream>
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 
49 
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 
53 {
54  if(!fManager) {
56  }
57  return fManager;
58 }
60 {
61  delete fManager;
62  fManager = 0;
63 }
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
68 : fVerbose(0), fNEvt1(-1), fNEvt2(-2),
69  fHistEMax (15.0*MeV), fHistEMin (0.),fHistNBin(100),
70  fTargetThresE(10*keV), fDetectorThresE(10*keV),fPulseWidth(1.*microsecond)
71 {
72  fEdepo.clear();
73  fHisto = new exrdmHisto();
74  BookHisto();
75 }
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 
80 {
81 #ifdef G4ANALYSIS_USE
82 #define HISTDELETE
83 #endif
84 #ifdef G4ANALYSIS_USE_ROOT
85 #define HISTDELETE
86 #endif
87 #ifdef HISTDELETE
88  delete fHisto;
89 #endif
90 }
91 
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
93 
95 {
96  fHistEMax = 15.0*MeV;
97  fHistEMin = .0*MeV;
98  fHistNBin = 100;
99 
100  fHisto->Add1D("H10",
101  "Energy deposit (MeV) in the traget",fHistNBin,fHistEMin,fHistEMax,MeV);
102  fHisto->Add1D("H11",
103  "Energy deposit (MeV) in the detector",fHistNBin,fHistEMin,fHistEMax,MeV);
104  fHisto->Add1D("H12",
105  "Total energy spectrum (MeV) of the traget and detector",fHistNBin,
107  fHisto->Add1D("H13",
108  "Coincidence spectrum (MeV) between the traget and detector",fHistNBin,
110  fHisto->Add1D("H14",
111  "Anti-coincidence spectrum (MeV) in the traget",fHistNBin,
113  fHisto->Add1D("H15",
114  "Anti-coincidence spectrum (MeV) in the detector",fHistNBin,
116  fHisto->Add1D("H16",
117  "Decay emission spectrum (MeV)",fHistNBin,fHistEMin,fHistEMax,MeV);
118  // in aida these histos are indiced from 0-6
119  //
120  fHisto->AddTuple( "T1", "Emitted Particles",
121  "double PID, Energy, Time, Weight" );
122  fHisto->AddTuple( "T2", "RadioIsotopes","double PID, Time, Weight" );
123  fHisto->AddTuple( "T3", "Energy Depositions","double Energy, Time, Weight" );
124  fHisto->AddTuple( "RDecayProducts", "All Products of RDecay",
125  "double PID, Z, A, Energy, Time, Weight" );
126 }
127 
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129 
131 {
132  fHisto->Book();
133  G4cout
134  << "exrdmAnalysisManager: Histograms are booked and the run has been started"
135  << G4endl;
138  pTable->FindProcess("RadioactiveDecay", "GenericIon");
139  if (rDecay != NULL) {
140  if (!rDecay->IsAnalogueMonteCarlo()) {
141  std::vector<G4RadioactivityTable*> theTables =
142  rDecay->GetTheRadioactivityTables();
143  for (size_t i = 0 ; i < theTables.size(); i++) {
144  theTables[i]->GetTheMap()->clear();
145  G4cout << " Clear the radioactivity map: 0, new size = "
146  << theTables[i]->GetTheMap()->size() << G4endl;
147  }
148  }
149  }
150 
151 }
152 
153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
154 
156 {
157  fHisto->Save();
158  G4cout << "exrdmAnalysisManager: Histograms have been saved!" << G4endl;
159 
160  // Output the induced radioactivities
161  // in VR mode only
162  //
165  pTable->FindProcess("RadioactiveDecay", "GenericIon");
166  if (rDecay != NULL) {
167  if (!rDecay->IsAnalogueMonteCarlo()) {
168  G4String fileName = fHisto->GetFileName() + ".activity" ;
169  std::ofstream outfile (fileName, std::ios::out );
170  std::vector<G4RadioactivityTable*> theTables =
171  rDecay->GetTheRadioactivityTables();
172  for (size_t i = 0 ; i < theTables.size(); i++) {
173  G4double rate, error;
174  outfile << "Radioactivities in decay window no. " << i << G4endl;
175  outfile
176  << "Z \tA \tE \tActivity (decays/window) \tError (decays/window) "
177  << G4endl;
178  map<G4ThreeVector,G4TwoVector> *aMap = theTables[i]->GetTheMap();
179  map<G4ThreeVector,G4TwoVector>::iterator iter;
180  for(iter=aMap->begin(); iter != aMap->end(); iter++) {
181  rate = iter->second.x()/nevent;
182  error = std::sqrt(iter->second.y())/nevent;
183  if ( rate < 0.) rate = 0.; // statically it can be < 0.
184  outfile << iter->first.x() <<"\t"<< iter->first.y() <<"\t"
185  << iter->first.z() << "\t" << rate <<"\t" << error
186  << G4endl;
187  }
188  outfile << G4endl;
189  }
190  outfile.close();
191  }
192  }
193 }
194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
195 
197 {
198  fEdepo.clear();
199 }
200 
201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
202 
204 {
205  if (fEdepo.size()) {
206  std::sort(fEdepo.begin(),fEdepo.end());
207  G4double TarE = fEdepo[0].GetEnergy();
208  G4double Time = fEdepo[0].GetTime();
209  G4double TarW = fEdepo[0].GetEnergy()*fEdepo[0].GetWeight();
210  G4double DetE = 0.;
211  G4double DetW = 0.;
212  G4double ComW = 0.;
213  if (TarE< 0.) {
214  DetE = -TarE;
215  DetW = -TarW;
216  TarE = 0.;
217  TarW = 0.;
218  }
219  for (size_t i = 1; i < fEdepo.size(); i++) {
220  if (std::fabs(fEdepo[i].GetTime() - Time) <= fPulseWidth) {
221  if ( fEdepo[i].GetEnergy() > 0. ) {
222  TarE += fEdepo[i].GetEnergy();
223  TarW += fEdepo[i].GetEnergy()*fEdepo[i].GetWeight();
224  } else {
225  DetE -= fEdepo[i].GetEnergy();
226  DetW -= fEdepo[i].GetEnergy()*fEdepo[i].GetWeight();
227  }
228  } else {
229  // tallying
230  if (TarE || DetE) ComW = (TarW+DetW)/(TarE+DetE);
231  if (TarE) TarW /= TarE;
232  if (DetE) DetW /= DetE;
233  //
234  if (TarE) fHisto->FillHisto(0,TarE,TarW); // target histogram
235  if (DetE) fHisto->FillHisto(1,DetE,DetW); // Detector histogram
236  if (TarE+DetE) fHisto->FillHisto(2,DetE+TarE,ComW); // Total histogram
237  if (DetE >= fDetectorThresE && TarE >= fTargetThresE )
238  fHisto->FillHisto(3,DetE,DetW); // coincidence histogram
239  if (TarE >= fTargetThresE && DetE < fDetectorThresE)
240  fHisto->FillHisto(4,TarE,TarW); // target anti-coincidence histogram
241  if (TarE < fTargetThresE && DetE >= fDetectorThresE)
242  fHisto->FillHisto(5,DetE,DetW); // detector anti-coincidence histogram
243  //
244  //reset
245  TarE = fEdepo[i].GetEnergy();
246  Time = fEdepo[i].GetTime();
247  TarW = fEdepo[i].GetEnergy()*fEdepo[i].GetWeight();
248  DetE = 0.;
249  DetW = 0.;
250  if (TarE < 0) {
251  DetE = -TarE;
252  DetW = -TarW;
253  TarE = 0.;
254  TarW = 0.;
255  }
256  }
257  }
258  //tally the last hit
259  if (TarE || DetE) ComW = (TarW+DetW)/(TarE+DetE);
260  if (TarE) TarW /= TarE;
261  if (DetE) DetW /= DetE;
262  //
263  if (TarE) fHisto->FillHisto(0,TarE,TarW); // target histogram
264  if (DetE) fHisto->FillHisto(1,DetE,DetW); // Detector histogram
265  if (TarE+DetE) fHisto->FillHisto(2,DetE+TarE,ComW); // Total histogram
266  if (DetE >= fDetectorThresE && TarE >= fTargetThresE )
267  fHisto->FillHisto(3,DetE,DetW); // coincidence histogram
268  if (TarE >= fTargetThresE && DetE < fDetectorThresE)
269  fHisto->FillHisto(4,TarE,TarW); // target anti-coincidence histogram
270  if (TarE < fTargetThresE && DetE >= fDetectorThresE)
271  fHisto->FillHisto(5,DetE,DetW); // detector anti-coincidence histogram
272  // now add zero energy to separat events
273  AddEnergy(0.,0.,0.);
274  }
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
278 
280  G4double time)
281 {
282  if(1 < fVerbose) {
283  G4cout << "exrdmAnalysisManager::AddEnergy: e(keV)= " << edep/keV
284  << " weight = " << weight << " time (s) = " << time/second
285  << G4endl;
286  }
287  fHisto->FillTuple(2, 0, edep/MeV);
288  fHisto->FillTuple(2,1,time/second);
289  fHisto->FillTuple(2,2,weight);
290  fHisto->AddRow(2);
291  //
292  exrdmEnergyDeposition A(edep,time,weight);
293  fEdepo.push_back(A);
294 }
295 
296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
297 
299  G4double weight, G4double time )
300 {
301  if(1 < fVerbose) {
302  G4cout << "exrdmAnalysisManager::AddParticle: " << pid
303  << G4endl;
304  }
305  fHisto->FillTuple(0,0, pid);
306  fHisto->FillTuple(0,1,energy/MeV);
307  fHisto->FillTuple(0,2,time/second);
308  fHisto->FillTuple(0,3,weight);
309  fHisto->AddRow(0);
310  // now fill th emission spectrum
311  if (energy>0.0) fHisto->FillHisto(6,energy/MeV,weight);
312 }
313 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
314 
316  G4double time )
317 {
318  if(1 < fVerbose) {
319  G4cout << "exrdmAnalysisManager::AddIsotope: " << pid
320  << G4endl;
321  }
322  fHisto->FillTuple(1,0,pid);
323  fHisto->FillTuple(1,1,time/second);
324  fHisto->FillTuple(1,2,weight);
325  fHisto->AddRow(1);
326 }
327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
328 
330  G4double energy, G4double time,
331  G4double weight)
332 {
333  fHisto->FillTuple(3,0,pid);
334  fHisto->FillTuple(3,1,double(Z));
335  fHisto->FillTuple(3,2,double(A));
336  fHisto->FillTuple(3,3,energy);
337  fHisto->FillTuple(3,4,time/s);
338  fHisto->FillTuple(3,5,weight);
339  fHisto->AddRow(3);
340 }
341 
342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
const G4String & GetFileName() const
Definition: exrdmHisto.cc:447
static const double MeV
Definition: G4SIunits.hh:193
void AddEnergy(G4double, G4double, G4double)
void FillTuple(G4int, const G4String &, G4double)
Definition: exrdmHisto.cc:375
Definition of the exrdmHisto class.
void FillHisto(G4int, G4double, G4double)
Definition: exrdmHisto.cc:274
void AddTuple(const G4String &, const G4String &, const G4String &)
Definition: exrdmHisto.cc:333
int G4int
Definition: G4Types.hh:78
static const double s
Definition: G4SIunits.hh:150
void AddRow(G4int)
Definition: exrdmHisto.cc:419
static const double microsecond
Definition: G4SIunits.hh:141
G4GLOB_DLL std::ostream G4cout
Definition of the exrdmAnalysisManager class.
std::vector< G4RadioactivityTable * > GetTheRadioactivityTables()
void AddDecayProduct(G4double pid, G4int Z, G4int A, G4double energy, G4double time, G4double weight)
static const double second
Definition: G4SIunits.hh:138
void AddParticle(G4double, G4double, G4double, G4double)
static const G4double A[nN]
std::vector< exrdmEnergyDeposition > fEdepo
static exrdmAnalysisManager * fManager
void Book()
Definition: exrdmHisto.cc:113
G4double energy(const ThreeVector &p, const G4double m)
static exrdmAnalysisManager * GetInstance()
#define G4endl
Definition: G4ios.hh:61
static PROLOG_HANDLER error
Definition: xmlrole.cc:108
static const double keV
Definition: G4SIunits.hh:195
double G4double
Definition: G4Types.hh:76
static G4ProcessTable * GetProcessTable()
G4VProcess * FindProcess(const G4String &processName, const G4String &particleName) const
void Add1D(const G4String &, const G4String &, G4int nb=100, G4double x1=0., G4double x2=1., G4double u=1.)
Definition: exrdmHisto.cc:226
void Save()
Definition: exrdmHisto.cc:196
void AddIsotope(G4double, G4double, G4double)