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