| 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 // $Id: RunAction.cc 98772 2016-08-09 14:25:31Z 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
35 #include "DetectorConstruction.hh" 36 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh" 37 #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" 38 #include "G4ParticleDefinition.hh"
44 #include "G4PhysicalConstants.hh" << 39
45 #include "G4Positron.hh" <<
46 #include "G4Run.hh" 40 #include "G4Run.hh"
47 #include "G4RunManager.hh" 41 #include "G4RunManager.hh"
48 #include "G4SystemOfUnits.hh" <<
49 #include "G4UnitsTable.hh" 42 #include "G4UnitsTable.hh"
>> 43 #include "G4EmCalculator.hh"
>> 44
>> 45 #include "G4Gamma.hh"
>> 46 #include "G4Electron.hh"
>> 47 #include "G4Positron.hh"
50 48
>> 49 #include "G4SystemOfUnits.hh"
>> 50 #include "G4PhysicalConstants.hh"
51 #include <iomanip> 51 #include <iomanip>
52 52
53 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 54
55 RunAction::RunAction(DetectorConstruction* det 55 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
56 : fDetector(det), fPrimary(prim), fAnalysisM << 56 : G4UserRunAction(),
>> 57 fGamma(0), fElectron(0), fPositron(0),
>> 58 fDetector(det), fPrimary(prim), fProcCounter(0), fAnalysisManager(0),
>> 59 fTotalEventCount(0),
>> 60 fPhotonStats(), fElectronStats(), fPositronStats()
57 { 61 {
58 fGamma = G4Gamma::Gamma(); 62 fGamma = G4Gamma::Gamma();
59 fElectron = G4Electron::Electron(); 63 fElectron = G4Electron::Electron();
60 fPositron = G4Positron::Positron(); 64 fPositron = G4Positron::Positron();
61 65
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(); 66 BookHisto();
75 } 67 }
76 68
77 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 70
79 RunAction::~RunAction() {} << 71 RunAction::~RunAction()
>> 72 {}
80 73
81 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 75
83 void RunAction::BeginOfRunAction(const G4Run* 76 void RunAction::BeginOfRunAction(const G4Run* aRun)
84 { << 77 {
85 G4cout << "### Run " << aRun->GetRunID() << 78 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
86 << 79
87 auto accumulableManager = G4AccumulableManag <<
88 accumulableManager->Reset(); <<
89 <<
90 // save Rndm status 80 // save Rndm status
91 // G4RunManager::GetRunManager()->SetRandom 81 // G4RunManager::GetRunManager()->SetRandomNumberStore(false);
92 // CLHEP::HepRandom::showEngineStatus(); 82 // CLHEP::HepRandom::showEngineStatus();
93 << 83
>> 84 if (fProcCounter) delete fProcCounter;
>> 85 fProcCounter = new ProcessesCount;
94 fTotalEventCount = 0; 86 fTotalEventCount = 0;
95 << 87 fPhotonStats.Clear();
96 // Open file histogram file << 88 fElectronStats.Clear();
97 fAnalysisManager->OpenFile(); << 89 fPositronStats.Clear();
98 } 90 }
99 91
100 //....oooOO0OOooo........oooOO0OOooo........oo 92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 93
102 void RunAction::FillData(const G4ParticleDefin << 94 void RunAction::FillData(const G4ParticleDefinition* particle,
103 G4double costheta, G4 << 95 G4double kinEnergy, G4double costheta,
>> 96 G4double phi,
>> 97 G4double longitudinalPolarization)
104 { 98 {
105 auto accManager = G4AccumulableManager::Inst <<
106 G4int id = -1; 99 G4int id = -1;
107 if (particle == fGamma) { << 100 if (particle == fGamma) {
108 dynamic_cast<ParticleStatistics*>(accManag << 101 fPhotonStats.FillData(kinEnergy, costheta, longitudinalPolarization);
109 ->FillData(kinEnergy, costheta, longitud << 102 if(fAnalysisManager) { id = 1; }
110 if (fAnalysisManager) { << 103 } else if (particle == fElectron) {
111 id = 1; << 104 fElectronStats.FillData(kinEnergy, costheta, longitudinalPolarization);
112 } << 105 if(fAnalysisManager) { id = 5; }
>> 106 } else if (particle == fPositron) {
>> 107 fPositronStats.FillData(kinEnergy, costheta, longitudinalPolarization);
>> 108 if(fAnalysisManager) { id = 9; }
113 } 109 }
114 else if (particle == fElectron) { << 110 if(id > 0) {
115 dynamic_cast<ParticleStatistics*>(accManag << 111 fAnalysisManager->FillH1(id,kinEnergy,1.0);
116 ->FillData(kinEnergy, costheta, longitud << 112 fAnalysisManager->FillH1(id+1,costheta,1.0);
117 if (fAnalysisManager) { << 113 fAnalysisManager->FillH1(id+2,phi,1.0);
118 id = 5; << 114 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 } 115 }
134 } 116 }
135 117
136 //....oooOO0OOooo........oooOO0OOooo........oo 118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137 119
138 void RunAction::BookHisto() 120 void RunAction::BookHisto()
139 { 121 {
140 // Always creating analysis manager 122 // Always creating analysis manager
141 fAnalysisManager = G4AnalysisManager::Instan 123 fAnalysisManager = G4AnalysisManager::Instance();
142 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 delete fAnalysisManager;
>> 173 fAnalysisManager = 0;
>> 174 }
206 } 175 }
207 176
208 //....oooOO0OOooo........oooOO0OOooo........oo 177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
209 178
210 void RunAction::CountProcesses(G4String& procN << 179 void RunAction::CountProcesses(G4String procName)
211 { 180 {
212 auto accManager = G4AccumulableManager::Inst << 181 // is the process already counted ?
213 dynamic_cast<ProcessesCount*>(accManager->Ge << 182 // *AS* change to std::map?!
>> 183 size_t nbProc = fProcCounter->size();
>> 184 size_t i = 0;
>> 185 while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++;
>> 186 if (i == nbProc) fProcCounter->push_back( new ProcessCount(procName));
>> 187
>> 188 (*fProcCounter)[i]->Count();
214 } 189 }
215 190
216 //....oooOO0OOooo........oooOO0OOooo........oo 191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
217 192
218 void RunAction::EndOfRunAction(const G4Run* aR 193 void RunAction::EndOfRunAction(const G4Run* aRun)
219 { 194 {
220 // Total number of events in run (all threa << 195 G4int NbOfEvents = aRun->GetNumberOfEvent();
221 G4int NbOfEvents = aRun->GetNumberOfEventToB <<
222 // G4int NbOfEvents = aRun->GetNumberOfEven <<
223 <<
224 if (NbOfEvents == 0) return; 196 if (NbOfEvents == 0) return;
225 << 197
226 G4int prec = G4cout.precision(5); << 198 G4int prec = G4cout.precision(5);
227 << 199
228 G4Material* material = fDetector->GetMateria 200 G4Material* material = fDetector->GetMaterial();
229 G4double density = material->GetDensity(); 201 G4double density = material->GetDensity();
230 << 202
231 if (fPrimary != nullptr) { << 203 G4ParticleDefinition* particle =
232 G4ParticleDefinition* particle = fPrimary- << 204 fPrimary->GetParticleGun()->GetParticleDefinition();
233 G4String Particle = particle->GetParticleN << 205 G4String Particle = particle->GetParticleName();
234 G4double energy = fPrimary->GetParticleGun << 206 G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
235 G4cout << "\n The run consists of " << fTo << 207 G4cout << "\n The run consists of " << NbOfEvents << " "<< Particle << " of "
236 << G4BestUnit(energy, "Energy") << << 208 << G4BestUnit(energy,"Energy") << " through "
237 << G4BestUnit(fDetector->GetBoxSize << 209 << G4BestUnit(fDetector->GetBoxSizeZ(),"Length") << " of "
238 << " (density: " << G4BestUnit(dens << 210 << material->GetName() << " (density: "
239 } << 211 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
240 // cross check from G4EmCalculator << 212
241 // G4cout << "\n Verification from G4EmCal << 213 //frequency of processes
242 // G4EmCalculator emCal; << 214 G4cout << "\n Process calls frequency --->\n";
243 << 215 for (size_t i=0; i< fProcCounter->size();i++) {
244 auto accManager = G4AccumulableManager::Inst << 216 G4String procName = (*fProcCounter)[i]->GetName();
245 accManager->Merge(); << 217 G4int count = (*fProcCounter)[i]->GetCounter();
246 << 218 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 } 219 }
>> 220
>> 221 if (fTotalEventCount == 0) return;
>> 222
>> 223 G4cout<<" Gamma: \n";
>> 224 fPhotonStats.PrintResults(fTotalEventCount);
>> 225 G4cout<<" Electron: \n";
>> 226 fElectronStats.PrintResults(fTotalEventCount);
>> 227 G4cout<<" Positron: \n";
>> 228 fPositronStats.PrintResults(fTotalEventCount);
>> 229
>> 230 //cross check from G4EmCalculator
>> 231 // G4cout << "\n Verification from G4EmCalculator. \n";
>> 232 // G4EmCalculator emCal;
262 233
263 // restore default format << 234 //restore default format
264 G4cout.precision(prec); << 235 G4cout.precision(prec);
265 236
266 // write out histograms << 237 // write out histograms
267 SaveHisto(NbOfEvents); 238 SaveHisto(NbOfEvents);
268 239
269 if (IsMaster()) { << 240 // show Rndm status
270 // show Rndm status << 241 CLHEP::HepRandom::showEngineStatus();
271 CLHEP::HepRandom::showEngineStatus(); <<
272 } <<
273 } 242 }
274 243
275 //....oooOO0OOooo........oooOO0OOooo........oo 244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
276 245
277 void RunAction::EventFinished() 246 void RunAction::EventFinished()
278 { 247 {
279 auto accManager = G4AccumulableManager::Inst <<
280 <<
281 ++fTotalEventCount; 248 ++fTotalEventCount;
282 << 249 fPhotonStats.EventFinished();
283 dynamic_cast<ParticleStatistics*>(accManager << 250 fElectronStats.EventFinished();
284 dynamic_cast<ParticleStatistics*>(accManager << 251 fPositronStats.EventFinished();
285 dynamic_cast<ParticleStatistics*>(accManager <<
286 } 252 }
287 253
288 //....oooOO0OOooo........oooOO0OOooo........oo 254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289 255
290 RunAction::ParticleStatistics::ParticleStatist << 256 RunAction::ParticleStatistics::ParticleStatistics()
291 : G4VAccumulable(name), << 257 : fCurrentNumber(0),
292 fCurrentNumber(0), << 258 fTotalNumber(0), fTotalNumber2(0),
293 fTotalNumber(0), << 259 fSumEnergy(0), fSumEnergy2(0),
294 fTotalNumber2(0), << 260 fSumPolarization(0), fSumPolarization2(0),
295 fSumEnergy(0), << 261 fSumCosTheta(0), fSumCosTheta2(0)
296 fSumEnergy2(0), <<
297 fSumPolarization(0), <<
298 fSumPolarization2(0), <<
299 fSumCosTheta(0), <<
300 fSumCosTheta2(0) <<
301 {} 262 {}
302 263
303 //....oooOO0OOooo........oooOO0OOooo........oo 264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
304 265
305 RunAction::ParticleStatistics::~ParticleStatis << 266 RunAction::ParticleStatistics::~ParticleStatistics()
>> 267 {}
306 268
307 //....oooOO0OOooo........oooOO0OOooo........oo 269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308 270
309 void RunAction::ParticleStatistics::EventFinis 271 void RunAction::ParticleStatistics::EventFinished()
310 { 272 {
311 fTotalNumber += fCurrentNumber; << 273 fTotalNumber+=fCurrentNumber;
312 fTotalNumber2 += fCurrentNumber * fCurrentNu << 274 fTotalNumber2+=fCurrentNumber*fCurrentNumber;
313 fCurrentNumber = 0; << 275 fCurrentNumber=0;
314 } 276 }
315 //....oooOO0OOooo........oooOO0OOooo........oo 277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
316 278
317 void RunAction::ParticleStatistics::FillData(G << 279 void RunAction::ParticleStatistics:: FillData(G4double kinEnergy,
318 G << 280 G4double costheta,
>> 281 G4double longitudinalPolarization)
319 { 282 {
320 ++fCurrentNumber; 283 ++fCurrentNumber;
321 fSumEnergy += kinEnergy; << 284 fSumEnergy+=kinEnergy;
322 fSumEnergy2 += kinEnergy * kinEnergy; << 285 fSumEnergy2+=kinEnergy*kinEnergy;
323 fSumPolarization += longitudinalPolarization << 286 fSumPolarization+=longitudinalPolarization;
324 fSumPolarization2 += longitudinalPolarizatio << 287 fSumPolarization2+=longitudinalPolarization*longitudinalPolarization;
325 fSumCosTheta += costheta; << 288 fSumCosTheta+=costheta;
326 fSumCosTheta2 += costheta * costheta; << 289 fSumCosTheta2+=costheta*costheta;
327 } 290 }
328 291
329 //....oooOO0OOooo........oooOO0OOooo........oo 292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
330 293
331 void RunAction::ParticleStatistics::PrintResul 294 void RunAction::ParticleStatistics::PrintResults(G4int totalNumberOfEvents)
332 { 295 {
333 G4cout << "Mean Number per Event :" << G4dou << 296 G4cout<<"Mean Number per Event :"
334 << "\n"; << 297 <<G4double(fTotalNumber)/G4double(totalNumberOfEvents)<<"\n";
335 if (fTotalNumber == 0) fTotalNumber = 1; << 298 if (fTotalNumber==0) fTotalNumber=1;
336 G4double energyMean = fSumEnergy / fTotalNum << 299 G4double energyMean=fSumEnergy/fTotalNumber;
337 G4double energyRms = std::sqrt(fSumEnergy2 / << 300 G4double energyRms=std::sqrt(fSumEnergy2/fTotalNumber-energyMean*energyMean);
338 G4cout << "Mean Energy :" << G4BestUnit(ener << 301 G4cout<<"Mean Energy :"<< G4BestUnit(energyMean,"Energy")
339 << G4BestUnit(energyRms, "Energy") << << 302 <<" +- "<<G4BestUnit(energyRms,"Energy")<<"\n";
340 G4double polarizationMean = fSumPolarization << 303 G4double polarizationMean=fSumPolarization/fTotalNumber;
341 G4double polarizationRms = << 304 G4double polarizationRms=
342 std::sqrt(fSumPolarization2 / fTotalNumber << 305 std::sqrt(fSumPolarization2/fTotalNumber-polarizationMean*polarizationMean);
343 G4cout << "Mean Polarization :" << polarizat << 306 G4cout<<"Mean Polarization :"<< polarizationMean
344 } << 307 <<" +- "<<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 } 308 }
356 309
357 //....oooOO0OOooo........oooOO0OOooo........oo 310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
358 311
359 void RunAction::ParticleStatistics::Merge(cons << 312 void RunAction::ParticleStatistics::Clear()
360 { << 313 {
361 auto rstat = dynamic_cast<const RunAction::P << 314 fCurrentNumber=0;
362 << 315 fTotalNumber=fTotalNumber2=0;
363 fCurrentNumber += rstat.fCurrentNumber; << 316 fSumEnergy=fSumEnergy2=0;
364 fTotalNumber += rstat.fTotalNumber; << 317 fSumPolarization=fSumPolarization2=0;
365 fTotalNumber2 += rstat.fTotalNumber2; << 318 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 } 319 }
373 320
374 //....oooOO0OOooo........oooOO0OOooo........oo 321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
375 322