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

Last change on this file since 850 was 807, checked in by garnier, 16 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.