| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /// \file electromagnetic/TestEm12/src/RunActi 26 /// \file electromagnetic/TestEm12/src/RunAction.cc
27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction class
28 // 28 //
29 // << 29 // $Id: RunAction.cc 69108 2013-04-18 13:10:10Z gcosmo $
>> 30 //
30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 33
33 #include "RunAction.hh" 34 #include "RunAction.hh"
34 <<
35 #include "DetectorConstruction.hh" 35 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh" <<
37 #include "PhysicsList.hh" 36 #include "PhysicsList.hh"
38 #include "PrimaryGeneratorAction.hh" <<
39 #include "Run.hh" <<
40 #include "StepMax.hh" 37 #include "StepMax.hh"
>> 38 #include "PrimaryGeneratorAction.hh"
>> 39 #include "HistoManager.hh"
41 40
>> 41 #include "G4Run.hh"
>> 42 #include "G4RunManager.hh"
>> 43 #include "G4UnitsTable.hh"
42 #include "G4EmCalculator.hh" 44 #include "G4EmCalculator.hh"
43 #include "G4EmParameters.hh" << 45
44 #include "Randomize.hh" 46 #include "Randomize.hh"
>> 47 #include "G4SystemOfUnits.hh"
>> 48 #include <iomanip>
45 49
46 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 51
48 RunAction::RunAction(DetectorConstruction* det << 52 RunAction::RunAction(DetectorConstruction* det, PhysicsList* phys,
49 : fDetector(det), fPhysics(phys), fPrimary(k << 53 PrimaryGeneratorAction* kin)
>> 54 :G4UserRunAction(),
>> 55 fDetector(det),fPhysics(phys),fKinematic(kin),fHistoManager(0)
50 { 56 {
51 // Book predefined histograms 57 // Book predefined histograms
52 fHistoManager = new HistoManager(); 58 fHistoManager = new HistoManager();
53 } << 59 }
54 60
55 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 62
57 RunAction::~RunAction() 63 RunAction::~RunAction()
58 { 64 {
59 delete fHistoManager; << 65 delete fHistoManager;
60 } <<
61 <<
62 //....oooOO0OOooo........oooOO0OOooo........oo <<
63 <<
64 G4Run* RunAction::GenerateRun() <<
65 { <<
66 fRun = new Run(fDetector); <<
67 return fRun; <<
68 } 66 }
69 67
70 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 69
72 void RunAction::BeginOfRunAction(const G4Run*) << 70 void RunAction::BeginOfRunAction(const G4Run* aRun)
73 { << 71 {
74 // show Rndm status << 72 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
75 if (isMaster) { << 73
76 G4Random::showEngineStatus(); << 74 // save Rndm status
77 G4EmParameters::Instance()->Dump(); << 75 ////G4RunManager::GetRunManager()->SetRandomNumberStore(true);
78 } << 76 CLHEP::HepRandom::showEngineStatus();
79 // histograms << 77
>> 78 //initialize total energy deposit
80 // 79 //
81 G4AnalysisManager* analysisManager = G4Analy << 80 fEdeposit = fEdeposit2 = 0.;
82 if (analysisManager->IsActive()) { << 81
83 analysisManager->OpenFile(); << 82 //initialize track legth of primary
84 } << 83 //
85 << 84 fTrackLen = fTrackLen2 = 0.;
86 if (!fPrimary) return; << 85
87 << 86 //initialize projected range
88 // keep run condition << 87 //
89 G4ParticleDefinition* particle = fPrimary->G << 88 fProjRange = fProjRange2 = 0.;
90 G4double energy = fPrimary->GetParticleGun() << 89
91 fRun->SetPrimary(particle, energy); << 90 //initialize mean step size
92 << 91 //
93 // get CsdaRange from EmCalculator << 92 fNbOfSteps = fNbOfSteps2 = 0; fStepSize = fStepSize2 = 0.;
>> 93
>> 94 //get fCsdaRange from EmCalculator
94 // 95 //
95 G4EmCalculator emCalculator; 96 G4EmCalculator emCalculator;
96 G4Material* material = fDetector->GetAbsorMa 97 G4Material* material = fDetector->GetAbsorMaterial();
97 << 98 G4ParticleDefinition* particle = fKinematic->GetParticleGun()
98 G4double csdaRange = DBL_MAX; << 99 ->GetParticleDefinition();
99 if (particle->GetPDGCharge() != 0.) { << 100 G4double energy = fKinematic->GetParticleGun()->GetParticleEnergy();
100 csdaRange = emCalculator.GetCSDARange(ener << 101 fCsdaRange = DBL_MAX;
101 fRun->SetCsdaRange(csdaRange); << 102 if (particle->GetPDGCharge() != 0.)
>> 103 fCsdaRange = emCalculator.GetCSDARange(energy,particle,material);
>> 104
>> 105 //histograms
>> 106 //
>> 107 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
>> 108 if ( analysisManager->IsActive() ) {
>> 109 analysisManager->OpenFile();
102 } 110 }
103 << 111
104 // set StepMax from histos 1 and 8 << 112 //set StepMax from histos 1 and 8
105 // 113 //
106 G4double stepMax = csdaRange; << 114 G4double stepMax = fCsdaRange;
107 G4int ih = 1; 115 G4int ih = 1;
108 if (analysisManager->GetH1Activation(ih)) 116 if (analysisManager->GetH1Activation(ih))
109 stepMax = analysisManager->GetH1Width(ih) << 117 stepMax = analysisManager->GetH1Width(ih);
110 ih = 8; 118 ih = 8;
111 if (analysisManager->GetH1Activation(ih)) { 119 if (analysisManager->GetH1Activation(ih)) {
112 G4double width = analysisManager->GetH1Wid 120 G4double width = analysisManager->GetH1Width(ih);
113 stepMax = std::min(stepMax, width * csdaRa << 121 stepMax = std::min(stepMax, width*fCsdaRange);
114 } << 122 }
115 fPhysics->GetStepMaxProcess()->SetMaxStep2(s 123 fPhysics->GetStepMaxProcess()->SetMaxStep2(stepMax);
116 << 124 G4cout << "\n---> stepMax from histos 1 and 8 = "
117 /// G4cout << "\n---> stepMax from histos 1 << 125 << G4BestUnit(stepMax,"Length") << G4endl;
118 /// << G4BestUnit(stepMax,"Length") < <<
119 } 126 }
120 127
121 //....oooOO0OOooo........oooOO0OOooo........oo 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 129
123 void RunAction::EndOfRunAction(const G4Run*) << 130 void RunAction::EndOfRunAction(const G4Run* aRun)
124 { 131 {
125 if (isMaster) fRun->EndOfRun(); << 132 std::ios::fmtflags mode = G4cout.flags();
>> 133 G4cout.setf(std::ios::fixed,std::ios::floatfield);
>> 134
>> 135 G4int NbofEvents = aRun->GetNumberOfEvent();
>> 136 if (NbofEvents == 0) return;
>> 137 G4double fNbofEvents = double(NbofEvents);
126 138
>> 139 //run conditions
>> 140 //
>> 141 G4Material* material = fDetector->GetAbsorMaterial();
>> 142 G4double density = material->GetDensity();
>> 143
>> 144 G4ParticleDefinition* particle = fKinematic->GetParticleGun()
>> 145 ->GetParticleDefinition();
>> 146 G4String partName = particle->GetParticleName();
>> 147 G4double energy = fKinematic->GetParticleGun()->GetParticleEnergy();
>> 148
>> 149 G4cout << "\n ======================== run summary ======================\n";
>> 150
>> 151 G4int prec = G4cout.precision(2);
>> 152
>> 153 G4cout << "\n The run consists of " << NbofEvents << " "<< partName << " of "
>> 154 << G4BestUnit(energy,"Energy") << " through "
>> 155 << G4BestUnit(fDetector->GetAbsorRadius(),"Length") << " of "
>> 156 << material->GetName() << " (density: "
>> 157 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
>> 158
>> 159 G4cout << "\n ============================================================\n";
>> 160
>> 161 //compute total energy deposit
>> 162 //
>> 163 fEdeposit /= NbofEvents; fEdeposit2 /= NbofEvents;
>> 164 G4double rms = fEdeposit2 - fEdeposit*fEdeposit;
>> 165 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
>> 166
>> 167 G4cout.precision(3);
>> 168 G4cout
>> 169 << "\n Total Energy deposited = " << G4BestUnit(fEdeposit,"Energy")
>> 170 << " +- " << G4BestUnit( rms,"Energy")
>> 171 << G4endl;
>> 172
>> 173 //compute track length of primary track
>> 174 //
>> 175 fTrackLen /= NbofEvents; fTrackLen2 /= NbofEvents;
>> 176 rms = fTrackLen2 - fTrackLen*fTrackLen;
>> 177 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
>> 178
>> 179 G4cout.precision(3);
>> 180 G4cout
>> 181 << "\n Track length of primary track = " << G4BestUnit(fTrackLen,"Length")
>> 182 << " +- " << G4BestUnit( rms,"Length");
>> 183
>> 184 //compare with csda range
>> 185 //
>> 186 //G4EmCalculator emCalculator;
>> 187 //G4double fCsdaRange = 0.;
>> 188 //if (particle->GetPDGCharge() != 0.)
>> 189 // fCsdaRange = emCalculator.GetCSDARange(energy,particle,material);
>> 190 G4cout
>> 191 << "\n Range from EmCalculator = " << G4BestUnit(fCsdaRange,"Length")
>> 192 << " (from full dE/dx)" << G4endl;
>> 193
>> 194 //compute projected range of primary track
>> 195 //
>> 196 fProjRange /= NbofEvents; fProjRange2 /= NbofEvents;
>> 197 rms = fProjRange2 - fProjRange*fProjRange;
>> 198 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
>> 199
>> 200 G4cout
>> 201 << "\n Projected range = " << G4BestUnit(fProjRange,"Length")
>> 202 << " +- " << G4BestUnit( rms,"Length")
>> 203 << G4endl;
>> 204
>> 205 //nb of steps and step size of primary track
>> 206 //
>> 207 G4double fNbSteps = fNbOfSteps/fNbofEvents, fNbSteps2 = fNbOfSteps2/fNbofEvents;
>> 208 rms = fNbSteps2 - fNbSteps*fNbSteps;
>> 209 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
>> 210
>> 211 G4cout.precision(2);
>> 212 G4cout << "\n Nb of steps of primary track = " << fNbSteps << " +- " << rms;
>> 213
>> 214 fStepSize /= NbofEvents; fStepSize2 /= NbofEvents;
>> 215 rms = fStepSize2 - fStepSize*fStepSize;
>> 216 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
>> 217
>> 218 G4cout.precision(3);
>> 219 G4cout
>> 220 << "\t Step size= " << G4BestUnit(fStepSize,"Length")
>> 221 << " +- " << G4BestUnit( rms,"Length")
>> 222 << G4endl;
>> 223
>> 224 // normalize histogram of longitudinal energy profile
>> 225 //
>> 226 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
>> 227 G4int ih = 1;
>> 228 G4double binWidth = analysisManager->GetH1Width(ih);
>> 229 G4double fac = (1./(NbofEvents*binWidth))*(mm/MeV);
>> 230 analysisManager->ScaleH1(ih,fac);
>> 231
>> 232 // normalize histogram d(E/E0)/d(r/r0)
>> 233 //
>> 234 ih = 8;
>> 235 binWidth = analysisManager->GetH1Width(ih);
>> 236 fac = 1./(NbofEvents*binWidth*energy);
>> 237 analysisManager->ScaleH1(ih,fac);
>> 238
>> 239 // reset default formats
>> 240 G4cout.setf(mode,std::ios::floatfield);
>> 241 G4cout.precision(prec);
>> 242
127 // save histograms 243 // save histograms
128 G4AnalysisManager* analysisManager = G4Analy << 244 if ( analysisManager->IsActive() ) {
129 if (analysisManager->IsActive()) { <<
130 analysisManager->Write(); 245 analysisManager->Write();
131 analysisManager->CloseFile(); 246 analysisManager->CloseFile();
132 } << 247 }
133 << 248
134 // show Rndm status 249 // show Rndm status
135 if (isMaster) G4Random::showEngineStatus(); << 250 CLHEP::HepRandom::showEngineStatus();
136 } 251 }
137 252
138 //....oooOO0OOooo........oooOO0OOooo........oo 253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
139 254