source: trunk/examples/extended/electromagnetic/TestEm7/src/RunAction.cc @ 812

Last change on this file since 812 was 807, checked in by garnier, 16 years ago

update
File size: 8.1 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// $Id: RunAction.cc,v 1.21 2008/01/14 12:11:39 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31
32#include "RunAction.hh"
33#include "DetectorConstruction.hh"
34#include "PhysicsList.hh"
35#include "StepMax.hh"
36#include "PrimaryGeneratorAction.hh"
37
38#include "G4Run.hh"
39#include "G4RunManager.hh"
40#include "G4UnitsTable.hh"
41#include "G4ios.hh"
42
43#include "Randomize.hh"
44
45#ifdef G4ANALYSIS_USE
46#include "AIDA/AIDA.h"
47#endif
48
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50
51RunAction::RunAction(DetectorConstruction* det, PhysicsList* phys,
52                     PrimaryGeneratorAction* kin)
53:detector(det), physics(phys), kinematic(kin), af(0), tree(0)
54{ 
55  tallyEdep = new G4double[MaxTally];
56  binLength = offsetX = 0.;
57  histo[0] = 0;
58 
59#ifdef G4ANALYSIS_USE
60  // Creating the analysis factory
61  af = AIDA_createAnalysisFactory();
62  if(!af) {
63    G4cout << "RunAction::RunAction() :" 
64           << " problem creating the AIDA analysis factory."
65           << G4endl;
66  }       
67#endif 
68}
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71
72RunAction::~RunAction()
73{
74  delete [] tallyEdep;
75 
76#ifdef G4ANALYSIS_USE
77  delete af; 
78#endif     
79}
80
81//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82
83void RunAction::bookHisto()
84{
85  length  = detector->GetAbsorSizeX();
86  G4double stepMax = physics->GetStepMaxProcess()->GetMaxStep();
87  const G4int nbmin = 100;
88  G4int nbBins = (int)(0.5 + length/stepMax);
89  if (nbBins < nbmin) nbBins = nbmin;
90  binLength = length/nbBins;
91  offsetX   = 0.5*length;
92 
93#ifdef G4ANALYSIS_USE
94  if (!af) return;
95
96  // Create a tree mapped to an hbook file.
97  G4bool readOnly  = false;
98  G4bool createNew = true;
99  G4String options = "--noErrors uncompress";
100  AIDA::ITreeFactory* tf  = af->createTreeFactory(); 
101  tree = tf->create("testem7.hbook","hbook", readOnly, createNew, options);
102  //tree = tf->create("testem7.root", "root",readOnly, createNew, options);
103  //tree = tf->create("testem7.XML" , "XML" ,readOnly, createNew, options);
104  delete tf;
105  if (!tree) {
106    G4cout << "RunAction::bookHisto()" << G4endl;
107    return;
108  }
109
110  // Create a histogram factory, whose histograms will be handled by the tree
111  AIDA::IHistogramFactory* hf = af->createHistogramFactory(*tree);
112 
113  // Create histogram
114  histo[0] = hf->createHistogram1D("1","Edep (MeV/mm) along absorber (mm)",
115             nbBins, 0, length/mm);
116             
117  delete hf;         
118  G4cout << "\n----> Histogram Tree opened" << G4endl;
119#endif
120}
121
122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123
124void RunAction::cleanHisto()
125{
126#ifdef G4ANALYSIS_USE
127  tree->commit();       // Writing the histograms to the file
128  tree->close();        // and closing the tree (and the file)
129  delete tree;
130  tree = 0;
131 
132  G4cout << "\n----> Histogram Tree saved" << G4endl; 
133#endif
134}
135
136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137
138void RunAction::FillHisto(G4int ih, G4double x, G4double weight)
139{
140#ifdef G4ANALYSIS_USE
141  if(histo[ih]) histo[ih]->fill(x, weight);
142#endif
143}
144
145//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146
147void RunAction::BeginOfRunAction(const G4Run* aRun)
148{ 
149  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
150 
151  // save Rndm status
152  G4RunManager::GetRunManager()->SetRandomNumberStore(true);
153  CLHEP::HepRandom::showEngineStatus();
154     
155  //initialize projected range, tallies, Ebeam, and book histograms
156  //
157  nPrimarySteps = 0;
158  projRange = projRange2 = 0.;
159  edeptot = eniel = 0.;
160  for (G4int j=0; j<MaxTally; j++) tallyEdep[j] = 0.;
161  kinematic->ResetEbeamCumul();
162  bookHisto();     
163}
164
165//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
166
167void RunAction::EndOfRunAction(const G4Run* aRun)
168{
169  G4int NbofEvents = aRun->GetNumberOfEvent();
170  if (NbofEvents == 0) return;
171
172  //run conditions
173  // 
174  G4Material* material = detector->GetAbsorMaterial();
175  G4double density = material->GetDensity();
176   
177  G4String particle = kinematic->GetParticleGun()->GetParticleDefinition()
178                      ->GetParticleName();   
179  G4double energy = kinematic->GetParticleGun()->GetParticleEnergy();
180  G4cout << "\n The run consists of " << NbofEvents << " "<< particle << " of "
181         << G4BestUnit(energy,"Energy") << " through " 
182         << G4BestUnit(detector->GetAbsorSizeX(),"Length") << " of "
183         << material->GetName() << " (density: " 
184         << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
185         
186  //compute projected range and straggling
187  //
188  projRange /= NbofEvents; projRange2 /= NbofEvents;
189  G4double rms = projRange2 - projRange*projRange;       
190  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
191
192  G4double nstep = G4double(nPrimarySteps)/G4double(NbofEvents);
193
194  G4cout.precision(6);       
195  G4cout << "\n Projected Range= "<< G4BestUnit(projRange,"Length")
196         << "   rms= "            << G4BestUnit( rms,"Length")
197         << G4endl;
198  G4cout << " Mean number of primary steps = "<< nstep << G4endl;
199
200  //compute energy deposition and NIEL
201  //
202  edeptot /= NbofEvents; 
203  G4cout << " Total energy deposit= "<< G4BestUnit(edeptot,"Energy")
204         << G4endl;
205  eniel /= NbofEvents; 
206  G4cout << " NIEL energy deposit = "<< G4BestUnit(eniel,"Energy")
207         << G4endl;
208     
209  //print dose in tallies
210  //
211  G4int tallyNumber = detector->GetTallyNumber();
212  if (tallyNumber > 0) {
213    G4double tallyMass = detector->GetTallyMass();
214    G4double Ebeam = kinematic->GetEbeamCumul();
215    G4cout << "\n---------------------------------------------------------\n";
216    G4cout << " Cumulated Doses : \tEdep      \tEdep/Ebeam \tDose" << G4endl;
217    for (G4int j=0; j<tallyNumber; j++) {
218      G4double Edep = tallyEdep[j], ratio = 100*Edep/Ebeam;
219      G4double Dose = Edep/tallyMass;
220      G4cout << " tally " << j << ": \t \t"
221             << G4BestUnit(Edep,"Energy") << "\t"
222             << ratio << " % \t"
223             << G4BestUnit(Dose,"Dose")   << G4endl;
224    }
225    G4cout << "\n---------------------------------------------------------\n";
226    G4cout << G4endl; 
227  }
228
229#ifdef G4ANALYSIS_USE
230  // normalize histogram
231  G4double fac = (mm/MeV)/(NbofEvents *  binLength);
232  histo[0]->scale(fac);
233#endif
234 
235 
236  // save and clean histo
237  cleanHisto();
238 
239  // show Rndm status
240  CLHEP::HepRandom::showEngineStatus();
241}
242
243//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.