| 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 polarisation/Pol01/src/RunAction.cc 26 /// \file polarisation/Pol01/src/RunAction.cc
27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction class
28 // 28 //
29 // << 29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 32
33 #include "RunAction.hh" 33 #include "RunAction.hh"
34 34
35 #include "DetectorConstruction.hh" 35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh" 36 #include "PrimaryGeneratorAction.hh"
37 #include "ProcessesCount.hh" <<
38 <<
39 #include "G4AccumulableManager.hh" <<
40 #include "G4Electron.hh" <<
41 #include "G4EmCalculator.hh" <<
42 #include "G4Gamma.hh" <<
43 #include "G4ParticleDefinition.hh" 37 #include "G4ParticleDefinition.hh"
44 #include "G4PhysicalConstants.hh" << 38
45 #include "G4Positron.hh" <<
46 #include "G4Run.hh" 39 #include "G4Run.hh"
47 #include "G4RunManager.hh" 40 #include "G4RunManager.hh"
48 #include "G4SystemOfUnits.hh" <<
49 #include "G4UnitsTable.hh" 41 #include "G4UnitsTable.hh"
>> 42 #include "G4EmCalculator.hh"
>> 43
>> 44 #include "G4Gamma.hh"
>> 45 #include "G4Electron.hh"
>> 46 #include "G4Positron.hh"
50 47
>> 48 #include "G4SystemOfUnits.hh"
>> 49 #include "G4PhysicalConstants.hh"
51 #include <iomanip> 50 #include <iomanip>
52 51
53 //....oooOO0OOooo........oooOO0OOooo........oo 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 53
55 RunAction::RunAction(DetectorConstruction* det 54 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
56 : fDetector(det), fPrimary(prim), fAnalysisM << 55 : G4UserRunAction(),
>> 56 fGamma(0), fElectron(0), fPositron(0),
>> 57 fDetector(det), fPrimary(prim), fProcCounter(0), fAnalysisManager(0),
>> 58 fTotalEventCount(0),
>> 59 fPhotonStats(), fElectronStats(), fPositronStats()
57 { 60 {
58 fGamma = G4Gamma::Gamma(); 61 fGamma = G4Gamma::Gamma();
59 fElectron = G4Electron::Electron(); 62 fElectron = G4Electron::Electron();
60 fPositron = G4Positron::Positron(); 63 fPositron = G4Positron::Positron();
61 64
62 auto accumulableManager = G4AccumulableManag <<
63 auto fPhotonStats = new ParticleStatistics(" <<
64 auto fElectronStats = new ParticleStatistics <<
65 auto fPositronStats = new ParticleStatistics <<
66 auto fProcCounter = new ProcessesCount("Proc <<
67 <<
68 accumulableManager->Register(fPhotonStats); <<
69 accumulableManager->Register(fElectronStats) <<
70 accumulableManager->Register(fPositronStats) <<
71 <<
72 accumulableManager->Register(fProcCounter); <<
73 <<
74 BookHisto(); 65 BookHisto();
75 } 66 }
76 67
77 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 69
79 RunAction::~RunAction() {} << 70 RunAction::~RunAction()
>> 71 {}
80 72
81 //....oooOO0OOooo........oooOO0OOooo........oo 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 74
83 void RunAction::BeginOfRunAction(const G4Run* 75 void RunAction::BeginOfRunAction(const G4Run* aRun)
84 { << 76 {
85 G4cout << "### Run " << aRun->GetRunID() << 77 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
86 << 78
87 auto accumulableManager = G4AccumulableManag <<
88 accumulableManager->Reset(); <<
89 <<
90 // save Rndm status 79 // save Rndm status
91 // G4RunManager::GetRunManager()->SetRandom 80 // G4RunManager::GetRunManager()->SetRandomNumberStore(false);
92 // CLHEP::HepRandom::showEngineStatus(); 81 // CLHEP::HepRandom::showEngineStatus();
93 << 82
>> 83 if (fProcCounter) delete fProcCounter;
>> 84 fProcCounter = new ProcessesCount;
94 fTotalEventCount = 0; 85 fTotalEventCount = 0;
95 << 86 fPhotonStats.Clear();
96 // Open file histogram file << 87 fElectronStats.Clear();
97 fAnalysisManager->OpenFile(); << 88 fPositronStats.Clear();
98 } 89 }
99 90
100 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 92
102 void RunAction::FillData(const G4ParticleDefin << 93 void RunAction::FillData(const G4ParticleDefinition* particle,
103 G4double costheta, G4 << 94 G4double kinEnergy, G4double costheta,
>> 95 G4double phi,
>> 96 G4double longitudinalPolarization)
104 { 97 {
105 auto accManager = G4AccumulableManager::Inst <<
106 G4int id = -1; 98 G4int id = -1;
107 if (particle == fGamma) { << 99 if (particle == fGamma) {
108 dynamic_cast<ParticleStatistics*>(accManag << 100 fPhotonStats.FillData(kinEnergy, costheta, longitudinalPolarization);
109 ->FillData(kinEnergy, costheta, longitud << 101 if(fAnalysisManager) { id = 1; }
110 if (fAnalysisManager) { << 102 } else if (particle == fElectron) {
111 id = 1; << 103 fElectronStats.FillData(kinEnergy, costheta, longitudinalPolarization);
112 } << 104 if(fAnalysisManager) { id = 5; }
>> 105 } else if (particle == fPositron) {
>> 106 fPositronStats.FillData(kinEnergy, costheta, longitudinalPolarization);
>> 107 if(fAnalysisManager) { id = 9; }
113 } 108 }
114 else if (particle == fElectron) { << 109 if(id > 0) {
115 dynamic_cast<ParticleStatistics*>(accManag << 110 fAnalysisManager->FillH1(id,kinEnergy,1.0);
116 ->FillData(kinEnergy, costheta, longitud << 111 fAnalysisManager->FillH1(id+1,costheta,1.0);
117 if (fAnalysisManager) { << 112 fAnalysisManager->FillH1(id+2,phi,1.0);
118 id = 5; << 113 fAnalysisManager->FillH1(id+3,longitudinalPolarization,1.0);
119 } <<
120 } <<
121 else if (particle == fPositron) { <<
122 dynamic_cast<ParticleStatistics*>(accManag <<
123 ->FillData(kinEnergy, costheta, longitud <<
124 if (fAnalysisManager) { <<
125 id = 9; <<
126 } <<
127 } <<
128 if (id > 0) { <<
129 fAnalysisManager->FillH1(id, kinEnergy, 1. <<
130 fAnalysisManager->FillH1(id + 1, costheta, <<
131 fAnalysisManager->FillH1(id + 2, phi, 1.0) <<
132 fAnalysisManager->FillH1(id + 3, longitudi <<
133 } 114 }
134 } 115 }
135 116
136 //....oooOO0OOooo........oooOO0OOooo........oo 117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137 118
138 void RunAction::BookHisto() 119 void RunAction::BookHisto()
139 { 120 {
140 // Always creating analysis manager 121 // Always creating analysis manager
141 fAnalysisManager = G4AnalysisManager::Instan 122 fAnalysisManager = G4AnalysisManager::Instance();
142 fAnalysisManager->SetDefaultFileType("root") 123 fAnalysisManager->SetDefaultFileType("root");
143 fAnalysisManager->SetActivation(true); 124 fAnalysisManager->SetActivation(true);
144 fAnalysisManager->SetVerboseLevel(1); 125 fAnalysisManager->SetVerboseLevel(1);
145 126
146 // Open file histogram file 127 // Open file histogram file
147 fAnalysisManager->SetFileName("pol01"); << 128 fAnalysisManager->OpenFile("pol01");
148 <<
149 fAnalysisManager->SetFirstHistoId(1); 129 fAnalysisManager->SetFirstHistoId(1);
150 130
151 // Creating an 1-dimensional histograms in t 131 // Creating an 1-dimensional histograms in the root directory of the tree
152 const G4String id[] = {"h1", "h2", "h3", "h4 << 132 const G4String id[] = { "h1", "h2", "h3", "h4", "h5",
153 const G4String title[] = { << 133 "h6", "h7", "h8", "h9", "h10", "h11", "h12"};
154 "Gamma Energy distribution", // 1 << 134 const G4String title[] =
155 "Gamma Cos(Theta) distribution", // 2 << 135 { "Gamma Energy distribution", //1
156 "Gamma Phi angular distribution", // 3 << 136 "Gamma Cos(Theta) distribution", //2
157 "Gamma longitudinal Polarization", // 4 << 137 "Gamma Phi angular distribution", //3
158 "Electron Energy distribution", // 5 << 138 "Gamma longitudinal Polarization", //4
159 "Electron Cos(Theta) distribution", // 6 << 139 "Electron Energy distribution", //5
160 "Electron Phi angular distribution", // 7 << 140 "Electron Cos(Theta) distribution", //6
161 "Electron longitudinal Polarization", // << 141 "Electron Phi angular distribution", //7
162 "Positron Energy distribution", // 9 << 142 "Electron longitudinal Polarization", //8
163 "Positron Cos(Theta) distribution", // 10 << 143 "Positron Energy distribution", //9
164 "Positron Phi angular distribution", // 1 << 144 "Positron Cos(Theta) distribution", //10
165 "Positron longitudinal Polarization" // 1 << 145 "Positron Phi angular distribution", //11
166 }; << 146 "Positron longitudinal Polarization" //12
>> 147 };
167 G4double vmin, vmax; 148 G4double vmin, vmax;
168 G4int nbins = 120; 149 G4int nbins = 120;
169 for (int i = 0; i < 12; ++i) { << 150 for(int i=0; i<12; ++i) {
170 G4int j = i - i / 4 * 4; << 151 G4int j = i - i/4*4;
171 if (0 == j) { << 152 if(0==j) { vmin = 0.; vmax = 12.*MeV; }
172 vmin = 0.; << 153 else if(1==j) { vmin = -1.; vmax = 1.; }
173 vmax = 12. * MeV; << 154 else if(2==j) { vmin = 0.; vmax = pi; }
174 } << 155 else { vmin = -1.5; vmax = 1.5; }
175 else if (1 == j) { << 156 G4int ih = fAnalysisManager->CreateH1(id[i],title[i],nbins,vmin,vmax);
176 vmin = -1.; <<
177 vmax = 1.; <<
178 } <<
179 else if (2 == j) { <<
180 vmin = 0.; <<
181 vmax = pi; <<
182 } <<
183 else { <<
184 vmin = -1.5; <<
185 vmax = 1.5; <<
186 } <<
187 G4int ih = fAnalysisManager->CreateH1(id[i <<
188 fAnalysisManager->SetH1Activation(ih, fals 157 fAnalysisManager->SetH1Activation(ih, false);
189 } 158 }
190 } 159 }
191 160
192 //....oooOO0OOooo........oooOO0OOooo........oo 161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
193 162
194 void RunAction::SaveHisto(G4int nevents) 163 void RunAction::SaveHisto(G4int nevents)
195 { 164 {
196 if (fAnalysisManager) { << 165 if(fAnalysisManager) {
197 if (IsMaster()) { << 166 G4double norm = 1.0/G4double(nevents);
198 G4double norm = 1.0 / G4double(nevents); << 167 for(int i=0; i<12; ++i) {
199 for (int i = 0; i < 12; ++i) { << 168 fAnalysisManager->ScaleH1(i, norm);
200 fAnalysisManager->ScaleH1(i, norm); <<
201 } <<
202 } 169 }
203 fAnalysisManager->Write(); 170 fAnalysisManager->Write();
204 fAnalysisManager->CloseFile(); 171 fAnalysisManager->CloseFile();
205 } << 172 }
206 } 173 }
207 174
208 //....oooOO0OOooo........oooOO0OOooo........oo 175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
209 176
210 void RunAction::CountProcesses(G4String& procN << 177 void RunAction::CountProcesses(G4String procName)
211 { 178 {
212 auto accManager = G4AccumulableManager::Inst << 179 // is the process already counted ?
213 dynamic_cast<ProcessesCount*>(accManager->Ge << 180 // *AS* change to std::map?!
>> 181 size_t nbProc = fProcCounter->size();
>> 182 size_t i = 0;
>> 183 while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++;
>> 184 if (i == nbProc) fProcCounter->push_back( new ProcessCount(procName));
>> 185
>> 186 (*fProcCounter)[i]->Count();
214 } 187 }
215 188
216 //....oooOO0OOooo........oooOO0OOooo........oo 189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
217 190
218 void RunAction::EndOfRunAction(const G4Run* aR 191 void RunAction::EndOfRunAction(const G4Run* aRun)
219 { 192 {
220 // Total number of events in run (all threa << 193 G4int NbOfEvents = aRun->GetNumberOfEvent();
221 G4int NbOfEvents = aRun->GetNumberOfEventToB <<
222 // G4int NbOfEvents = aRun->GetNumberOfEven <<
223 <<
224 if (NbOfEvents == 0) return; 194 if (NbOfEvents == 0) return;
225 << 195
226 G4int prec = G4cout.precision(5); << 196 G4int prec = G4cout.precision(5);
227 << 197
228 G4Material* material = fDetector->GetMateria 198 G4Material* material = fDetector->GetMaterial();
229 G4double density = material->GetDensity(); 199 G4double density = material->GetDensity();
230 << 200
231 if (fPrimary != nullptr) { << 201 G4ParticleDefinition* particle =
232 G4ParticleDefinition* particle = fPrimary- << 202 fPrimary->GetParticleGun()->GetParticleDefinition();
233 G4String Particle = particle->GetParticleN << 203 G4String Particle = particle->GetParticleName();
234 G4double energy = fPrimary->GetParticleGun << 204 G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
235 G4cout << "\n The run consists of " << fTo << 205 G4cout << "\n The run consists of " << NbOfEvents << " "<< Particle << " of "
236 << G4BestUnit(energy, "Energy") << << 206 << G4BestUnit(energy,"Energy") << " through "
237 << G4BestUnit(fDetector->GetBoxSize << 207 << G4BestUnit(fDetector->GetBoxSizeZ(),"Length") << " of "
238 << " (density: " << G4BestUnit(dens << 208 << material->GetName() << " (density: "
239 } << 209 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
240 // cross check from G4EmCalculator << 210
241 // G4cout << "\n Verification from G4EmCal << 211 //frequency of processes
242 // G4EmCalculator emCal; << 212 G4cout << "\n Process calls frequency --->\n";
243 << 213 for (size_t i=0; i< fProcCounter->size();i++) {
244 auto accManager = G4AccumulableManager::Inst << 214 G4String procName = (*fProcCounter)[i]->GetName();
245 accManager->Merge(); << 215 G4int count = (*fProcCounter)[i]->GetCounter();
246 << 216 G4cout << "\t" << procName << " = " << count<<"\n";
247 if (IsMaster()) { <<
248 // frequency of processes <<
249 G4cout << "\n Process calls frequency ---> <<
250 dynamic_cast<ProcessesCount*>(accManager-> <<
251 G4cout << " Gamma: \n"; <<
252 dynamic_cast<ParticleStatistics*>(accManag <<
253 ->PrintResults(NbOfEvents); <<
254 G4cout << " Electron: \n"; <<
255 dynamic_cast<ParticleStatistics*>(accManag <<
256 ->PrintResults(NbOfEvents); <<
257 G4cout << " Positron: \n"; <<
258 dynamic_cast<ParticleStatistics*>(accManag <<
259 ->PrintResults(NbOfEvents); <<
260 G4cout << G4endl; <<
261 } 217 }
>> 218
>> 219 if (fTotalEventCount == 0) return;
>> 220
>> 221 G4cout<<" Gamma: \n";
>> 222 fPhotonStats.PrintResults(fTotalEventCount);
>> 223 G4cout<<" Electron: \n";
>> 224 fElectronStats.PrintResults(fTotalEventCount);
>> 225 G4cout<<" Positron: \n";
>> 226 fPositronStats.PrintResults(fTotalEventCount);
>> 227
>> 228 //cross check from G4EmCalculator
>> 229 // G4cout << "\n Verification from G4EmCalculator. \n";
>> 230 // G4EmCalculator emCal;
262 231
263 // restore default format << 232 //restore default format
264 G4cout.precision(prec); << 233 G4cout.precision(prec);
265 234
266 // write out histograms << 235 // write out histograms
267 SaveHisto(NbOfEvents); 236 SaveHisto(NbOfEvents);
268 237
269 if (IsMaster()) { << 238 // show Rndm status
270 // show Rndm status << 239 CLHEP::HepRandom::showEngineStatus();
271 CLHEP::HepRandom::showEngineStatus(); <<
272 } <<
273 } 240 }
274 241
275 //....oooOO0OOooo........oooOO0OOooo........oo 242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
276 243
277 void RunAction::EventFinished() 244 void RunAction::EventFinished()
278 { 245 {
279 auto accManager = G4AccumulableManager::Inst <<
280 <<
281 ++fTotalEventCount; 246 ++fTotalEventCount;
282 << 247 fPhotonStats.EventFinished();
283 dynamic_cast<ParticleStatistics*>(accManager << 248 fElectronStats.EventFinished();
284 dynamic_cast<ParticleStatistics*>(accManager << 249 fPositronStats.EventFinished();
285 dynamic_cast<ParticleStatistics*>(accManager <<
286 } 250 }
287 251
288 //....oooOO0OOooo........oooOO0OOooo........oo 252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289 253
290 RunAction::ParticleStatistics::ParticleStatist << 254 RunAction::ParticleStatistics::ParticleStatistics()
291 : G4VAccumulable(name), << 255 : fCurrentNumber(0),
292 fCurrentNumber(0), << 256 fTotalNumber(0), fTotalNumber2(0),
293 fTotalNumber(0), << 257 fSumEnergy(0), fSumEnergy2(0),
294 fTotalNumber2(0), << 258 fSumPolarization(0), fSumPolarization2(0),
295 fSumEnergy(0), << 259 fSumCosTheta(0), fSumCosTheta2(0)
296 fSumEnergy2(0), <<
297 fSumPolarization(0), <<
298 fSumPolarization2(0), <<
299 fSumCosTheta(0), <<
300 fSumCosTheta2(0) <<
301 {} 260 {}
302 261
303 //....oooOO0OOooo........oooOO0OOooo........oo 262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
304 263
305 RunAction::ParticleStatistics::~ParticleStatis << 264 RunAction::ParticleStatistics::~ParticleStatistics()
>> 265 {}
306 266
307 //....oooOO0OOooo........oooOO0OOooo........oo 267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308 268
309 void RunAction::ParticleStatistics::EventFinis 269 void RunAction::ParticleStatistics::EventFinished()
310 { 270 {
311 fTotalNumber += fCurrentNumber; << 271 fTotalNumber+=fCurrentNumber;
312 fTotalNumber2 += fCurrentNumber * fCurrentNu << 272 fTotalNumber2+=fCurrentNumber*fCurrentNumber;
313 fCurrentNumber = 0; << 273 fCurrentNumber=0;
314 } 274 }
315 //....oooOO0OOooo........oooOO0OOooo........oo 275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
316 276
317 void RunAction::ParticleStatistics::FillData(G << 277 void RunAction::ParticleStatistics:: FillData(G4double kinEnergy,
318 G << 278 G4double costheta,
>> 279 G4double longitudinalPolarization)
319 { 280 {
320 ++fCurrentNumber; 281 ++fCurrentNumber;
321 fSumEnergy += kinEnergy; << 282 fSumEnergy+=kinEnergy;
322 fSumEnergy2 += kinEnergy * kinEnergy; << 283 fSumEnergy2+=kinEnergy*kinEnergy;
323 fSumPolarization += longitudinalPolarization << 284 fSumPolarization+=longitudinalPolarization;
324 fSumPolarization2 += longitudinalPolarizatio << 285 fSumPolarization2+=longitudinalPolarization*longitudinalPolarization;
325 fSumCosTheta += costheta; << 286 fSumCosTheta+=costheta;
326 fSumCosTheta2 += costheta * costheta; << 287 fSumCosTheta2+=costheta*costheta;
327 } 288 }
328 289
329 //....oooOO0OOooo........oooOO0OOooo........oo 290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
330 291
331 void RunAction::ParticleStatistics::PrintResul 292 void RunAction::ParticleStatistics::PrintResults(G4int totalNumberOfEvents)
332 { 293 {
333 G4cout << "Mean Number per Event :" << G4dou << 294 G4cout<<"Mean Number per Event :"
334 << "\n"; << 295 <<G4double(fTotalNumber)/G4double(totalNumberOfEvents)<<"\n";
335 if (fTotalNumber == 0) fTotalNumber = 1; << 296 if (fTotalNumber==0) fTotalNumber=1;
336 G4double energyMean = fSumEnergy / fTotalNum << 297 G4double energyMean=fSumEnergy/fTotalNumber;
337 G4double energyRms = std::sqrt(fSumEnergy2 / << 298 G4double energyRms=std::sqrt(fSumEnergy2/fTotalNumber-energyMean*energyMean);
338 G4cout << "Mean Energy :" << G4BestUnit(ener << 299 G4cout<<"Mean Energy :"<< G4BestUnit(energyMean,"Energy")
339 << G4BestUnit(energyRms, "Energy") << << 300 <<" +- "<<G4BestUnit(energyRms,"Energy")<<"\n";
340 G4double polarizationMean = fSumPolarization << 301 G4double polarizationMean=fSumPolarization/fTotalNumber;
341 G4double polarizationRms = << 302 G4double polarizationRms=
342 std::sqrt(fSumPolarization2 / fTotalNumber << 303 std::sqrt(fSumPolarization2/fTotalNumber-polarizationMean*polarizationMean);
343 G4cout << "Mean Polarization :" << polarizat << 304 G4cout<<"Mean Polarization :"<< polarizationMean
344 } << 305 <<" +- "<<polarizationRms<<"\n";
345 <<
346 //....oooOO0OOooo........oooOO0OOooo........oo <<
347 <<
348 void RunAction::ParticleStatistics::Reset() <<
349 { <<
350 fCurrentNumber = 0; <<
351 fTotalNumber = fTotalNumber2 = 0; <<
352 fSumEnergy = fSumEnergy2 = 0; <<
353 fSumPolarization = fSumPolarization2 = 0; <<
354 fSumCosTheta = fSumCosTheta2 = 0; <<
355 } 306 }
356 307
357 //....oooOO0OOooo........oooOO0OOooo........oo 308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
358 309
359 void RunAction::ParticleStatistics::Merge(cons << 310 void RunAction::ParticleStatistics::Clear()
360 { << 311 {
361 auto rstat = dynamic_cast<const RunAction::P << 312 fCurrentNumber=0;
362 << 313 fTotalNumber=fTotalNumber2=0;
363 fCurrentNumber += rstat.fCurrentNumber; << 314 fSumEnergy=fSumEnergy2=0;
364 fTotalNumber += rstat.fTotalNumber; << 315 fSumPolarization=fSumPolarization2=0;
365 fTotalNumber2 += rstat.fTotalNumber2; << 316 fSumCosTheta=fSumCosTheta2=0;
366 fSumEnergy += rstat.fSumEnergy; <<
367 fSumEnergy2 += rstat.fSumEnergy2; <<
368 fSumPolarization += rstat.fSumPolarization; <<
369 fSumPolarization2 += rstat.fSumPolarization2 <<
370 fSumCosTheta += rstat.fSumCosTheta; <<
371 fSumCosTheta2 += rstat.fSumCosTheta2; <<
372 } 317 }
373 318
374 //....oooOO0OOooo........oooOO0OOooo........oo 319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
375 320