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