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

Last change on this file since 1285 was 1230, checked in by garnier, 14 years ago

update to geant4.9.3
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}
48void exrdmAnalysisManager::dispose()
49{
50  delete fManager;
51  fManager = 0;
52}
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();
65  bookHisto();
66}
67
68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69
70exrdmAnalysisManager::~exrdmAnalysisManager()
71{
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
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  //
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" );
110
111}
112
113//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
114
115void exrdmAnalysisManager::BeginOfRun()
116{
117  histo->book();
118  G4cout << "exrdmAnalysisManager: Histograms are booked and the run has been started" << G4endl;
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.