| 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 medical/electronScattering/src/RunAc 26 /// \file medical/electronScattering/src/RunAction.cc
27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction class
28 // 28 //
>> 29 // $Id$
29 // 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 "PrimaryGeneratorAction.hh" 36 #include "PrimaryGeneratorAction.hh"
38 #include "Run.hh" << 37 #include "HistoManager.hh"
39 38
40 #include "G4PhysicalConstants.hh" <<
41 #include "G4Run.hh" 39 #include "G4Run.hh"
42 #include "G4RunManager.hh" 40 #include "G4RunManager.hh"
43 #include "G4SystemOfUnits.hh" <<
44 #include "G4UnitsTable.hh" 41 #include "G4UnitsTable.hh"
45 #include "Randomize.hh" <<
46 42
>> 43 #include "G4PhysicalConstants.hh"
>> 44 #include "G4SystemOfUnits.hh"
>> 45 #include "Randomize.hh"
47 #include <iomanip> 46 #include <iomanip>
48 47
49 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 49
51 RunAction::RunAction(DetectorConstruction* det << 50 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin,
52 : fDetector(det), fPrimary(kin), fHistoManag << 51 HistoManager* histo)
53 { << 52 :fDetector(det), fPrimary(kin), fHistoManager(histo)
54 fHistoManager = new HistoManager(); << 53 { }
55 } <<
56 54
57 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 56
59 RunAction::~RunAction() 57 RunAction::~RunAction()
60 { << 58 { }
61 delete fHistoManager; <<
62 } <<
63 59
64 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65 61
66 G4Run* RunAction::GenerateRun() << 62 void RunAction::BeginOfRunAction(const G4Run* aRun)
67 { 63 {
68 fRun = new Run(fDetector); << 64 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
69 return fRun; << 65
>> 66 fHistoManager->book();
>> 67
>> 68 InitFluence();
>> 69
>> 70 // save Rndm status
>> 71 G4RunManager::GetRunManager()->SetRandomNumberStore(false);
>> 72 CLHEP::HepRandom::showEngineStatus();
70 } 73 }
71 74
72 //....oooOO0OOooo........oooOO0OOooo........oo 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73 76
74 void RunAction::BeginOfRunAction(const G4Run*) << 77 void RunAction::EndOfRunAction(const G4Run* aRun)
75 { 78 {
76 // save Rndm status << 79 G4int TotNbofEvents = aRun->GetNumberOfEvent();
77 G4RunManager::GetRunManager()->SetRandomNumb << 80 if (TotNbofEvents == 0) return;
>> 81
>> 82 //Scatter foil
>> 83 //
>> 84 G4Material* material = fDetector->GetMaterialScatter();
>> 85 G4double length = fDetector->GetThicknessScatter();
>> 86 G4double density = material->GetDensity();
>> 87
>> 88 G4ParticleDefinition* particle = fPrimary->GetParticleGun()
>> 89 ->GetParticleDefinition();
>> 90 G4String partName = particle->GetParticleName();
>> 91 G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
>> 92
>> 93 G4cout << "\n ======================== run summary ======================\n";
>> 94
>> 95 G4int prec = G4cout.precision(3);
>> 96
>> 97 G4cout << "\n The run was " << TotNbofEvents << " " << partName << " of "
>> 98 << G4BestUnit(energy,"Energy") << " through "
>> 99 << G4BestUnit(length,"Length") << " of "
>> 100 << material->GetName() << " (density: "
>> 101 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
>> 102
>> 103 G4cout.precision(prec);
>> 104
>> 105 // normalize histograms
>> 106 //
>> 107 G4double fac = 1./(double(TotNbofEvents));
>> 108 fHistoManager->Normalize(1,fac);
>> 109 fHistoManager->Normalize(2,fac);
>> 110 fHistoManager->Normalize(3,fac);
>> 111
>> 112 ComputeFluenceError();
>> 113 PrintFluence(TotNbofEvents);
>> 114
>> 115 // save histograms
>> 116 fHistoManager->save();
>> 117
>> 118 // show Rndm status
78 CLHEP::HepRandom::showEngineStatus(); 119 CLHEP::HepRandom::showEngineStatus();
>> 120 }
79 121
80 // keep run conductions << 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81 if (fPrimary) { <<
82 G4ParticleDefinition* particle = fPrimary- <<
83 G4double energy = fPrimary->GetParticleGun <<
84 fRun->SetPrimary(particle, energy); <<
85 } <<
86 123
87 // histograms << 124 void RunAction::InitFluence()
>> 125 {
>> 126 // create Fluence histo in any case
88 // 127 //
89 G4AnalysisManager* analysisManager = G4Analy << 128 G4int ih = 4;
90 if (analysisManager->IsActive()) { << 129 if (!fHistoManager->HistoExist(ih))
91 analysisManager->OpenFile(); << 130 fHistoManager->SetHisto(ih, 120, 0*mm, 240*mm, "mm");
92 } << 131
>> 132 //construct vectors for fluence distribution
>> 133 //
>> 134 fNbBins = fHistoManager->GetNbins(ih);
>> 135 fDr = fHistoManager->GetBinWidth(ih);
>> 136 fluence.resize(fNbBins, 0.);
>> 137 fluence1.resize(fNbBins, 0.);
>> 138 fluence2.resize(fNbBins, 0.);
>> 139 fNbEntries.resize(fNbBins, 0);
93 } 140 }
94 141
95 //....oooOO0OOooo........oooOO0OOooo........oo 142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
96 143
97 void RunAction::EndOfRunAction(const G4Run*) << 144 void RunAction::SumFluence(G4double r, G4double fl)
98 { 145 {
99 // compute and print statistic << 146 G4int ibin = (int)(r/fDr);
100 if (isMaster) { << 147 if (ibin >= fNbBins) return;
101 fRun->EndOfRun(); << 148 fNbEntries[ibin]++;
102 } << 149 fluence[ibin] += fl;
>> 150 fluence2[ibin] += fl*fl;
>> 151 }
103 152
104 // save histograms << 153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 G4AnalysisManager* analysisManager = G4Analy <<
106 154
107 if (analysisManager->IsActive()) { << 155 void RunAction::ComputeFluenceError()
108 analysisManager->Write(); << 156 {
109 analysisManager->CloseFile(); << 157 //compute rms
>> 158 //
>> 159 G4double ds,variance,rms;
>> 160 G4double rmean = -0.5*fDr;
>> 161
>> 162 for (G4int bin=0; bin<fNbBins; bin++) {
>> 163 rmean += fDr;
>> 164 ds = twopi*rmean*fDr;
>> 165 fluence[bin] /= ds;
>> 166 fluence2[bin] /= (ds*ds);
>> 167 variance = 0.;
>> 168 if (fNbEntries[bin] > 0)
>> 169 variance = fluence2[bin] - (fluence[bin]*fluence[bin])/fNbEntries[bin];
>> 170 rms = 0.;
>> 171 if(variance > 0.) rms = std::sqrt(variance);
>> 172 fluence2[bin] = rms;
>> 173 }
>> 174
>> 175 //normalize to first bins, compute error and fill histo
>> 176 //
>> 177 G4double rnorm(4*mm), radius(0.), fnorm(0.), fnorm2(0.);
>> 178 G4int inorm = -1;
>> 179 do {
>> 180 inorm++; radius += fDr; fnorm += fluence[inorm]; fnorm2 += fluence2[inorm];
>> 181 } while (radius < rnorm);
>> 182 fnorm /= (inorm+1);
>> 183 fnorm2 /= (inorm+1);
>> 184 //
>> 185 G4double ratio, error;
>> 186 G4double scale = 1./fnorm;
>> 187 G4double err0 = fnorm2/fnorm, err1 = 0.;
>> 188 //
>> 189 rmean = -0.5*fDr;
>> 190
>> 191 for (G4int bin=0; bin<fNbBins; bin++) {
>> 192 ratio = fluence[bin]*scale;
>> 193 error = 0.;
>> 194 if (ratio > 0.) {
>> 195 err1 = fluence2[bin]/fluence[bin];
>> 196 error = ratio*std::sqrt(err1*err1 + err0*err0);
>> 197 }
>> 198 fluence1[bin] = ratio;
>> 199 fluence2[bin] = error;
>> 200 rmean += fDr;
>> 201 fHistoManager->FillHisto(4,rmean,ratio);
110 } 202 }
>> 203 }
111 204
112 // show Rndm status << 205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113 if (isMaster) { << 206
114 CLHEP::HepRandom::showEngineStatus(); << 207 #include <fstream>
>> 208
>> 209 void RunAction::PrintFluence(G4int TotEvents)
>> 210 {
>> 211 G4String name = fHistoManager->GetFileName();
>> 212 G4String fileName = name + ".ascii";
>> 213 std::ofstream File(fileName, std::ios::out);
>> 214
>> 215 std::ios::fmtflags mode = File.flags();
>> 216 File.setf( std::ios::scientific, std::ios::floatfield );
>> 217 G4int prec = File.precision(3);
>> 218
>> 219 File << " Fluence density distribution \n "
>> 220 << "\n ibin \t radius (mm) \t Nb \t fluence\t norma fl\t rms/nfl (%) \n"
>> 221 << G4endl;
>> 222
>> 223 G4double rmean = -0.5*fDr;
>> 224 for (G4int bin=0; bin<fNbBins; bin++) {
>> 225 rmean +=fDr;
>> 226 G4double error = 0.;
>> 227 if (fluence1[bin] > 0.) error = 100*fluence2[bin]/fluence1[bin];
>> 228 File << " " << bin << "\t " << rmean/mm << "\t " << fNbEntries[bin]
>> 229 << "\t " << fluence[bin]/double(TotEvents) << "\t " << fluence1[bin]
>> 230 << "\t " << error
>> 231 << G4endl;
115 } 232 }
>> 233
>> 234 // restaure default formats
>> 235 File.setf(mode,std::ios::floatfield);
>> 236 File.precision(prec);
116 } 237 }
117 238
118 //....oooOO0OOooo........oooOO0OOooo........oo 239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 240
119 241