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

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

update to geant4.9.3
File size: 6.8 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.24 2008/08/22 18:30:27 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
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#include "Histo.hh"
45
46//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47
48RunAction::RunAction(DetectorConstruction* det, PhysicsList* phys,
49 PrimaryGeneratorAction* kin)
50:detector(det), physics(phys), kinematic(kin)
51{
52 tallyEdep = new G4double[MaxTally];
53 binLength = offsetX = 0.;
54 histo = new Histo();
55 histo->setFileName("testem7");
56 histo->add1D("1","Edep (MeV/mm) along absorber (mm)", 100, 0, 100);
57 histo->add1D("2","Edep (MeV/mm) along absorber zoomed (mm)", 100, 0, 100);
58 histo->add1D("3","Projectile range (mm)", 100, 0, 100);
59}
60
61//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62
63RunAction::~RunAction()
64{
65 delete [] tallyEdep;
66 delete histo;
67}
68
69//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70
71void RunAction::FillHisto(G4int ih, G4double x, G4double weight)
72{
73 histo->fill(ih, x, weight);
74}
75
76//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
77
78void RunAction::BeginOfRunAction(const G4Run* aRun)
79{
80 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
81
82 // save Rndm status
83 G4RunManager::GetRunManager()->SetRandomNumberStore(true);
84 CLHEP::HepRandom::showEngineStatus();
85
86 //initialize projected range, tallies, Ebeam, and book histograms
87 //
88 nPrimarySteps = 0;
89 nRange = 0;
90 projRange = projRange2 = 0.;
91 edeptot = eniel = 0.;
92 for (G4int j=0; j<MaxTally; j++) tallyEdep[j] = 0.;
93 kinematic->ResetEbeamCumul();
94
95 // define "1" histogram binning
96 length = detector->GetAbsorSizeX();
97 G4double stepMax = physics->GetStepMaxProcess()->GetMaxStep();
98 const G4int nbmin = 100;
99 G4int nbBins = (G4int)(0.5 + length/stepMax);
100 if (nbBins < nbmin) nbBins = nbmin;
101 binLength = length/nbBins;
102 offsetX = 0.5*length;
103
104 // histogram "1" is defined by the length of the target
105 // zoomed histograms are defined by UI command
106 histo->setHisto1D(0, nbBins, 0, length, mm);
107
108 histo->book();
109}
110
111//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
112
113void RunAction::EndOfRunAction(const G4Run* aRun)
114{
115 G4int NbofEvents = aRun->GetNumberOfEvent();
116 if (NbofEvents == 0) return;
117
118 //run conditions
119 //
120 G4Material* material = detector->GetAbsorMaterial();
121 G4double density = material->GetDensity();
122
123 G4String particle = kinematic->GetParticleGun()->GetParticleDefinition()
124 ->GetParticleName();
125 G4double energy = kinematic->GetParticleGun()->GetParticleEnergy();
126 G4cout << "\n The run consists of " << NbofEvents << " "<< particle << " of "
127 << G4BestUnit(energy,"Energy") << " through "
128 << G4BestUnit(detector->GetAbsorSizeX(),"Length") << " of "
129 << material->GetName() << " (density: "
130 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
131
132 //compute projected range and straggling
133 //
134 if(nRange > 0) {
135 projRange /= nRange;
136 projRange2 /= nRange;
137 }
138 G4double rms = projRange2 - projRange*projRange;
139 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
140
141 G4double nstep = G4double(nPrimarySteps)/G4double(NbofEvents);
142
143 G4cout.precision(6);
144 G4cout << "\n Projected Range= "<< G4BestUnit(projRange,"Length")
145 << " rms= " << G4BestUnit( rms,"Length")
146 << G4endl;
147 G4cout << " Mean number of primary steps = "<< nstep << G4endl;
148
149 //compute energy deposition and NIEL
150 //
151 edeptot /= NbofEvents;
152 G4cout << " Total energy deposit= "<< G4BestUnit(edeptot,"Energy")
153 << G4endl;
154 eniel /= NbofEvents;
155 G4cout << " NIEL energy deposit = "<< G4BestUnit(eniel,"Energy")
156 << G4endl;
157
158 //print dose in tallies
159 //
160 G4int tallyNumber = detector->GetTallyNumber();
161 if (tallyNumber > 0) {
162 G4double tallyMass = detector->GetTallyMass();
163 G4double Ebeam = kinematic->GetEbeamCumul();
164 G4cout << "\n---------------------------------------------------------\n";
165 G4cout << " Cumulated Doses : \tEdep \tEdep/Ebeam \tDose" << G4endl;
166 for (G4int j=0; j<tallyNumber; j++) {
167 G4double Edep = tallyEdep[j], ratio = 100*Edep/Ebeam;
168 G4double Dose = Edep/tallyMass;
169 G4cout << " tally " << j << ": \t \t"
170 << G4BestUnit(Edep,"Energy") << "\t"
171 << ratio << " % \t"
172 << G4BestUnit(Dose,"Dose") << G4endl;
173 }
174 G4cout << "\n---------------------------------------------------------\n";
175 G4cout << G4endl;
176 }
177
178 // normalize histogram
179 G4double fac = (mm/MeV)/(NbofEvents * binLength);
180 for (G4int j=0; j<3; j++) {histo->scale(j, fac);}
181
182 // save and clean histo
183 histo->save();
184
185 // show Rndm status
186 CLHEP::HepRandom::showEngineStatus();
187}
188
189//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.