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

Last change on this file since 1330 was 1230, checked in by garnier, 16 years ago

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