source: trunk/examples/extended/radioactivedecay/exrdm/src/exrdmAnalysisManager.cc@ 810

Last change on this file since 810 was 807, checked in by garnier, 17 years ago

update
File size: 8.6 KB
Line 
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
27
28//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
29//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
30
31#include "exrdmAnalysisManager.hh"
32#include "G4UnitsTable.hh"
33#include "exrdmHisto.hh"
34
35//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
36
37exrdmAnalysisManager* exrdmAnalysisManager::fManager = 0;
38
39//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40
41exrdmAnalysisManager* exrdmAnalysisManager::getInstance()
42{
43 if(!fManager) {
44 fManager = new exrdmAnalysisManager();
45 }
46 return fManager;
47}
48
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50
51exrdmAnalysisManager::exrdmAnalysisManager()
52{
53 verbose = 0;
54 nEvt1 = -1;
55 nEvt2 = -1;
56 targetThresE = 10*keV;
57 detectorThresE = 10*keV;
58 pulseWidth = 1.*microsecond;
59 histo = new exrdmHisto();
60#if defined G4ANALYSIS_USE_AIDA || defined G4ANALYSIS_USE_ROOT
61 bookHisto();
62#endif
63}
64
65//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66
67exrdmAnalysisManager::~exrdmAnalysisManager()
68{
69// delete histo;
70}
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
74void exrdmAnalysisManager::bookHisto()
75{
76 histEMax = 15.0*MeV;
77 histEMin = .0*MeV;
78 histNBin = 100;
79
80 histo->add1D("10",
81 "Energy deposit (MeV) in the traget",histNBin,histEMin,histEMax,MeV);
82 histo->add1D("11",
83 "Energy deposit (MeV) in the detector",histNBin,histEMin,histEMax,MeV);
84 histo->add1D("12",
85 "Total energy spectrum (MeV) of the traget and detector",histNBin,histEMin,histEMax,MeV);
86 histo->add1D("13",
87 "Coincidence spectrum (MeV) between the traget and detector",histNBin,histEMin,histEMax,MeV);
88 histo->add1D("14",
89 "Anti-coincidence spectrum (MeV) in the traget",histNBin,histEMin,histEMax,MeV);
90 histo->add1D("15",
91 "Anti-coincidence spectrum (MeV) in the detector",histNBin,histEMin,histEMax,MeV);
92 histo->add1D("16",
93 "Decay emission spectrum (MeV)",histNBin,histEMin,histEMax,MeV);
94 // in aida these histos are indiced from 0-6
95 //
96 histo->addTuple( "100", "Emitted Particles","float PID, Energy, Time, Weight" );
97 histo->addTuple( "200", "RadioIsotopes","float PID, Time, Weight" );
98 histo->addTuple( "300", "Energy Depositions","float Energy, Time, Weight" );
99
100}
101
102//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
103
104void exrdmAnalysisManager::BeginOfRun()
105{
106#if defined G4ANALYSIS_USE_AIDA || G4ANALYSIS_USE_ROOT
107 histo->book();
108#endif
109 if(verbose > 0) {
110 G4cout << "exrdmAnalysisManager: Histograms are booked and the run has been started"
111 << G4endl;
112 }
113}
114
115//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116
117void exrdmAnalysisManager::EndOfRun()
118{
119#if defined G4ANALYSIS_USE_AIDA || G4ANALYSIS_USE_ROOT
120 histo->save();
121#endif
122}
123
124//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125
126void exrdmAnalysisManager::BeginOfEvent()
127{
128 Edepo.clear();
129}
130
131//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
132
133void exrdmAnalysisManager::EndOfEvent()
134{
135 if (Edepo.size()) {
136 std::sort(Edepo.begin(),Edepo.end());
137 G4double TarE = Edepo[0].GetEnergy();
138 G4double Time = Edepo[0].GetTime();
139 G4double TarW = Edepo[0].GetEnergy()*Edepo[0].GetWeight();
140 G4double DetE = 0.;
141 G4double DetW = 0.;
142 G4double ComW = 0.;
143 if (TarE< 0.) {
144 DetE = -TarE;
145 DetW = -TarW;
146 TarE = 0.;
147 TarW = 0.;
148 }
149 for (size_t i = 1; i < Edepo.size(); i++) {
150 if (std::fabs((Edepo[i].GetTime()- Time)/second) <= pulseWidth) {
151 if ( Edepo[i].GetEnergy() > 0. ) {
152 TarE += Edepo[i].GetEnergy();
153 TarW += Edepo[i].GetEnergy()*Edepo[i].GetWeight();
154 } else {
155 DetE -= Edepo[i].GetEnergy();
156 DetW -= Edepo[i].GetEnergy()*Edepo[i].GetWeight();
157 }
158 } else {
159 // tallying
160 if (TarE || DetE) ComW = (TarW+DetW)/(TarE+DetE);
161 if (TarE) TarW /= TarE;
162 if (DetE) DetW /= DetE;
163 //
164 if (TarE) histo->fillHisto(0,TarE,TarW); // target histogram
165 if (DetE) histo->fillHisto(1,DetE,DetW); // Detector histogram
166 if (TarE+DetE) histo->fillHisto(2,DetE+TarE,ComW); // Total histogram
167 if (DetE >= detectorThresE && TarE >= targetThresE )
168 histo->fillHisto(3,DetE,DetW); // coincidence histogram
169 if (TarE >= targetThresE && DetE < detectorThresE)
170 histo->fillHisto(4,TarE,TarW); // target anti-coincidence histogram
171 if (TarE < targetThresE && DetE >= detectorThresE)
172 histo->fillHisto(5,DetE,DetW); // detector anti-coincidence histogram
173 //
174 //reset
175 TarE = Edepo[i].GetEnergy();
176 Time = Edepo[i].GetTime();
177 TarW = Edepo[i].GetEnergy()*Edepo[i].GetWeight();
178 DetE = 0.;
179 DetW = 0.;
180 if (TarE < 0) {
181 DetE = -TarE;
182 DetW = -TarW;
183 TarE = 0.;
184 TarW = 0.;
185 }
186 }
187 }
188 //tally the last hit
189 if (TarE || DetE) ComW = (TarW+DetW)/(TarE+DetE);
190 if (TarE) TarW /= TarE;
191 if (DetE) DetW /= DetE;
192 //
193 if (TarE) histo->fillHisto(0,TarE,TarW); // target histogram
194 if (DetE) histo->fillHisto(1,DetE,DetW); // Detector histogram
195 if (TarE+DetE) histo->fillHisto(2,DetE+TarE,ComW); // Total histogram
196 if (DetE >= detectorThresE && TarE >= targetThresE )
197 histo->fillHisto(3,DetE,DetW); // coincidence histogram
198 if (TarE >= targetThresE && DetE < detectorThresE)
199 histo->fillHisto(4,TarE,TarW); // target anti-coincidence histogram
200 if (TarE < targetThresE && DetE >= detectorThresE)
201 histo->fillHisto(5,DetE,DetW); // detector anti-coincidence histogram
202 // now add zero energy to separat events
203 AddEnergy(0.,0.,0.);
204 }
205}
206
207//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
208
209void exrdmAnalysisManager::AddEnergy(G4double edep, G4double weight, G4double time)
210{
211 if(1 < verbose) {
212 G4cout << "exrdmAnalysisManager::AddEnergy: e(keV)= " << edep/keV
213 << " weight = " << weight << " time (s) = " << time/second
214 << G4endl;
215 }
216 histo->fillTuple(2, 0, edep/MeV);
217 histo->fillTuple(2,1,time/second);
218 histo->fillTuple(2,2,weight);
219 histo->addRow(2);
220 //
221 exrdmEnergyDeposition A(edep,time,weight);
222 Edepo.push_back(A);
223}
224
225//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
226
227void exrdmAnalysisManager::AddParticle(G4double pid, G4double energy, G4double weight, G4double time )
228{
229 if(1 < verbose) {
230 G4cout << "exrdmAnalysisManager::AddParticle: " << pid
231 << G4endl;
232 }
233 histo->fillTuple(0,0, pid);
234 histo->fillTuple(0,1,energy/MeV);
235 histo->fillTuple(0,2,time/second);
236 histo->fillTuple(0,3,weight);
237 histo->addRow(0);
238 // now fill th emission spectrum
239 if (energy>0.0) histo->fillHisto(6,energy/MeV,weight);
240}
241//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
242
243void exrdmAnalysisManager::AddIsotope(G4double pid,G4double weight, G4double time )
244{
245 if(1 < verbose) {
246 G4cout << "exrdmAnalysisManager::AddIsotope: " << pid
247 << G4endl;
248 }
249 histo->fillTuple(1,0,pid);
250 histo->fillTuple(1,1,time/second);
251 histo->fillTuple(1,2,weight);
252 histo->addRow(1);
253}
254//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
255
Note: See TracBrowser for help on using the repository browser.