ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RunAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RunAction.cc
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 //
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 
33 #include "RunAction.hh"
34 #include "DetectorConstruction.hh"
35 #include "PhysicsList.hh"
36 #include "StepMax.hh"
37 #include "PrimaryGeneratorAction.hh"
38 
39 #include "G4Run.hh"
40 #include "G4RunManager.hh"
41 #include "G4UnitsTable.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4ios.hh"
44 
45 #include "Randomize.hh"
46 
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 
51  : G4UserRunAction(),
52  fAnalysisManager(0), fDetector(det), fPhysics(phys), fKinematic(kin),
53  fTallyEdep(new G4double[kMaxTally]), fProjRange(0.), fProjRange2(0.),
54  fEdeptot(0.), fEniel(0.), fNbPrimarySteps(0), fRange(0)
55 {
56  // Book predefined histograms
57  BookHisto();
58 }
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
63 {
64  delete [] fTallyEdep;
65 }
66 
67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68 
69 void RunAction::BeginOfRunAction(const G4Run* aRun)
70 {
71  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
72 
73  if(!fAnalysisManager) { BookHisto(); }
74 
76 
77  //initialize projected range, tallies, Ebeam, and book histograms
78  //
79  fNbPrimarySteps = 0;
80  fRange = 0;
81  fProjRange = fProjRange2 = 0.;
82  fEdeptot = fEniel = 0.;
83  for (G4int j=0; j<kMaxTally; ++j) { fTallyEdep[j] = 0.; }
85 
86  if (fAnalysisManager->IsActive()) {
88 
89  // histogram "1" is defined by the length of the target
90  // zoomed histograms are defined by UI command
92  G4double stepMax = fPhysics->GetStepMaxProcess()->GetMaxStep();
93  G4int nbmin = 100;
94  G4int nbBins = (G4int)(0.5 + length/stepMax);
95  if (nbBins < nbmin) nbBins = nbmin;
96  fAnalysisManager->SetH1(1, nbBins, 0., length, "mm");
97  }
98 }
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 
102 void RunAction::EndOfRunAction(const G4Run* aRun)
103 {
104  G4int nbofEvents = aRun->GetNumberOfEvent();
105  if (nbofEvents == 0) return;
106 
107  //run conditions
108  //
110  G4double density = material->GetDensity();
111 
113  ->GetParticleName();
115  G4cout << "\n The run consists of " << nbofEvents << " "<< particle << " of "
116  << G4BestUnit(energy,"Energy") << " through "
117  << G4BestUnit(fDetector->GetAbsorSizeX(),"Length") << " of "
118  << material->GetName() << " (density: "
119  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
120 
121  //compute projected range and straggling
122  //
123  if(fRange > 0) {
124  fProjRange /= fRange;
125  fProjRange2 /= fRange;
126  }
128  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
129 
130  G4double nstep = G4double(fNbPrimarySteps)/G4double(nbofEvents);
131 
132  G4cout.precision(6);
133  G4cout << "\n Projected Range= "<< G4BestUnit(fProjRange,"Length")
134  << " rms= " << G4BestUnit( rms,"Length")
135  << G4endl;
136  G4cout << " Mean number of primary steps = "<< nstep << G4endl;
137 
138  //compute energy deposition and niel
139  //
140  fEdeptot /= nbofEvents;
141  G4cout << " Total energy deposit= "<< G4BestUnit(fEdeptot,"Energy")
142  << G4endl;
143  fEniel /= nbofEvents;
144  G4cout << " niel energy deposit = "<< G4BestUnit(fEniel,"Energy")
145  << G4endl;
146 
147  //print dose in tallies
148  //
149  G4int tallyNumber = fDetector->GetTallyNumber();
150  if (tallyNumber > 0) {
151  G4double Ebeam = fKinematic->GetEbeamCumul();
152  G4cout << "\n---------------------------------------------------------\n";
153  G4cout << " Cumulated Doses : \tEdep \tEdep/Ebeam \tDose" << G4endl;
154  for (G4int j=0; j < tallyNumber; ++j) {
155  G4double Edep = fTallyEdep[j], ratio = 100*Edep/Ebeam;
156  G4double tallyMass = fDetector->GetTallyMass(j);
157  G4double Dose = Edep/tallyMass;
158  G4cout << " tally " << j << ": \t \t"
159  << G4BestUnit(Edep,"Energy") << "\t"
160  << ratio << " % \t"
161  << G4BestUnit(Dose,"Dose") << G4endl;
162  }
163  G4cout << "\n---------------------------------------------------------\n";
164  G4cout << G4endl;
165  }
166 
167  if (fAnalysisManager->IsActive() ) {
168  // normalize histograms
169  //
170  for (G4int j=1; j<3; ++j) {
171  G4double binWidth = fAnalysisManager->GetH1Width(j);
172  G4double fac = (mm/MeV)/(nbofEvents * binWidth);
173  fAnalysisManager->ScaleH1(j, fac);
174  }
175  fAnalysisManager->ScaleH1(3, 1./nbofEvents);
176 
177  // save histograms
180  delete fAnalysisManager;
181  fAnalysisManager = 0;
182  }
183 
184  // show Rndm status
185  //
187 }
188 
189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
190 
192 {
193  // Create or get analysis manager
194  // The choice of analysis technology is done via selection of a namespace
195  // in HistoManager.hh
197  fAnalysisManager->SetFileName("testem7");
199  fAnalysisManager->SetActivation(true); // enable inactivation of histograms
200 
201  // Define histograms start values
202  const G4int kMaxHisto = 4;
203  const G4String id[] = { "h0", "h1", "h2", "h3" };
204  const G4String title[] =
205  { "dummy", //0
206  "Edep (MeV/mm) along absorber ", //1
207  "Edep (MeV/mm) along absorber zoomed", //2
208  "projectile range" //3
209  };
210 
211  // Default values (to be reset via /analysis/h1/set command)
212  G4int nbins = 100;
213  G4double vmin = 0.;
214  G4double vmax = 100.;
215 
216  // Create all histograms as inactivated
217  // as we have not yet set nbins, vmin, vmax
218  for (G4int k=0; k<kMaxHisto; ++k) {
219  G4int ih = fAnalysisManager->CreateH1(id[k], title[k], nbins, vmin, vmax);
220  G4bool activ = false;
221  if (k == 1) activ = true;
222  fAnalysisManager->SetH1Activation(ih, activ);
223  }
224 }
225 
226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......