Geant4_10
Run.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
28 //
29 // $Id: Run.cc 75577 2013-11-04 12:03:26Z vnivanch $
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "Run.hh"
35 
36 #include "DetectorConstruction.hh"
37 #include "PrimaryGeneratorAction.hh"
38 #include "RunActionMessenger.hh"
39 #include "EmAcceptance.hh"
40 
41 #include "G4Run.hh"
42 #include "G4RunManager.hh"
43 #include "G4UnitsTable.hh"
44 #include "G4Threading.hh"
45 
46 #include "G4SystemOfUnits.hh"
47 #include <iomanip>
48 
49 #include "Randomize.hh"
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52 
54  :G4Run(),fDet(det),fKin(kin),
55  f_nLbin(MaxBin),f_nRbin(MaxBin)
56 {
57  Reset();
58 }
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
62 void Run::Reset()
63 {
64  f_nLbin = fDet->GetnLtot();
65  f_dEdL.resize(f_nLbin);
66  fSumELongit.resize(f_nLbin);
67  fSumELongitCumul.resize(f_nLbin);
68  fSumE2Longit.resize(f_nLbin);
69  fSumE2LongitCumul.resize(f_nLbin);
70 
71  f_nRbin = fDet->GetnRtot();
72  f_dEdR.resize(f_nRbin);
73  fSumERadial.resize(f_nRbin);
74  fSumERadialCumul.resize(f_nRbin);
75  fSumE2Radial.resize(f_nRbin);
76  fSumE2RadialCumul.resize(f_nRbin);
77 
78  fChargedStep = 0.0;
79  fNeutralStep = 0.0;
80 
81  fVerbose = 0;
82 
83  //initialize arrays of cumulative energy deposition
84  //
85  for (G4int i=0; i<f_nLbin; i++)
86  fSumELongit[i]=fSumE2Longit[i]=fSumELongitCumul[i]=fSumE2LongitCumul[i]=0.;
87 
88  for (G4int j=0; j<f_nRbin; j++)
89  fSumERadial[j]=fSumE2Radial[j]=fSumERadialCumul[j]=fSumE2RadialCumul[j]=0.;
90 
91  //initialize track length
92  fSumChargTrLength=fSum2ChargTrLength=fSumNeutrTrLength=fSum2NeutrTrLength=0.;
93 
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97 
98 Run::~Run()
99 {}
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102 
104 {
105  //initialize arrays of energy deposit per bin
106  for (G4int i=0; i<f_nLbin; i++)
107  { f_dEdL[i] = 0.; }
108 
109  for (G4int j=0; j<f_nRbin; j++)
110  { f_dEdR[j] = 0.; }
111 
112  //initialize tracklength
113  fChargTrLength = fNeutrTrLength = 0.;
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 
119 {
120  //accumulate statistic
121  //
122  G4double dLCumul = 0.;
123  for (G4int i=0; i<f_nLbin; i++)
124  {
125  fSumELongit[i] += f_dEdL[i];
126  fSumE2Longit[i] += f_dEdL[i]*f_dEdL[i];
127  dLCumul += f_dEdL[i];
128  fSumELongitCumul[i] += dLCumul;
129  fSumE2LongitCumul[i] += dLCumul*dLCumul;
130  }
131 
132  G4double dRCumul = 0.;
133  for (G4int j=0; j<f_nRbin; j++)
134  {
135  fSumERadial[j] += f_dEdR[j];
136  fSumE2Radial[j] += f_dEdR[j]*f_dEdR[j];
137  dRCumul += f_dEdR[j];
138  fSumERadialCumul[j] += dRCumul;
139  fSumE2RadialCumul[j] += dRCumul*dRCumul;
140  }
141 
142  fSumChargTrLength += fChargTrLength;
143  fSum2ChargTrLength += fChargTrLength*fChargTrLength;
144  fSumNeutrTrLength += fNeutrTrLength;
145  fSum2NeutrTrLength += fNeutrTrLength*fNeutrTrLength;
146 
147  //fill histograms
148  //
149 
150  G4double Ekin=fKin->GetParticleGun()->GetParticleEnergy();
152  G4double radl=fDet->GetMaterial()->GetRadlen();
153 
154  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
155  analysisManager->FillH1(1, 100.*dLCumul/(Ekin+mass));
156  analysisManager->FillH1(2, fChargTrLength/radl);
157  analysisManager->FillH1(3, fNeutrTrLength/radl);
158 }
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
161 
162 void Run::Merge(const G4Run* run)
163 {
164  const Run* localRun = static_cast<const Run*>(run);
165 
166  fChargedStep += localRun->fChargedStep;
167  fNeutralStep += localRun->fNeutralStep;
168 
169  for (G4int i=0; i<f_nLbin; i++) {
170  fSumELongit[i] += localRun->fSumELongit[i];
171  fSumE2Longit[i] += localRun->fSumE2Longit[i];
172  fSumELongitCumul[i] += localRun->fSumELongitCumul[i];
173  fSumE2LongitCumul[i] += localRun->fSumE2LongitCumul[i];
174  }
175 
176  for (G4int j=0; j<f_nRbin; j++) {
177  fSumERadial[j] += localRun->fSumERadial[j];
178  fSumE2Radial[j] += localRun->fSumE2Radial[j];
179  fSumERadialCumul[j] += localRun->fSumERadialCumul[j];
180  fSumE2RadialCumul[j] += localRun->fSumE2RadialCumul[j];
181  }
182 
183  fSumChargTrLength += localRun->fSumChargTrLength;
184  fSum2ChargTrLength += localRun->fSum2ChargTrLength;
185  fSumNeutrTrLength += localRun->fSumNeutrTrLength;
186  fSum2NeutrTrLength += localRun->fSum2NeutrTrLength;
187 
188  G4Run::Merge(run);
189 }
190 
191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
192 
194 {
195  G4int NbOfEvents = GetNumberOfEvent();
196 
197  G4double kinEnergy = fKin->GetParticleGun()->GetParticleEnergy();
198  assert(NbOfEvents*kinEnergy > 0);
199 
200  fChargedStep /= G4double(NbOfEvents);
201  fNeutralStep /= G4double(NbOfEvents);
202 
204  G4double norme = 100./(NbOfEvents*(kinEnergy+mass));
205 
206  //longitudinal
207  //
208  G4double dLradl = fDet->GetdLradl();
209 
210  MyVector MeanELongit(f_nLbin), rmsELongit(f_nLbin);
211  MyVector MeanELongitCumul(f_nLbin), rmsELongitCumul(f_nLbin);
212 
213  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
214 
215  G4int i;
216  for (i=0; i<f_nLbin; i++) {
217  MeanELongit[i] = norme*fSumELongit[i];
218  rmsELongit[i] =
219  norme*std::sqrt(std::abs(NbOfEvents*fSumE2Longit[i]
220  - fSumELongit[i]*fSumELongit[i]));
221 
222  MeanELongitCumul[i] = norme*fSumELongitCumul[i];
223  rmsELongitCumul[i] = norme*std::sqrt(std::abs(NbOfEvents*
224  fSumE2LongitCumul[i] - fSumELongitCumul[i]*fSumELongitCumul[i]));
225  G4double bin = (i+0.5)*dLradl;
226  analysisManager->FillH1(4, bin,MeanELongit[i]/dLradl);
227  analysisManager->FillH1(5, bin, rmsELongit[i]/dLradl);
228  bin = (i+1)*dLradl;
229  analysisManager->FillH1(6, bin,MeanELongitCumul[i]);
230  analysisManager->FillH1(7, bin, rmsELongitCumul[i]);
231  }
232 
233  //radial
234  //
235  G4double dRradl = fDet->GetdRradl();
236 
237  MyVector MeanERadial(f_nRbin), rmsERadial(f_nRbin);
238  MyVector MeanERadialCumul(f_nRbin), rmsERadialCumul(f_nRbin);
239 
240  for (i=0; i<f_nRbin; i++) {
241  MeanERadial[i] = norme*fSumERadial[i];
242  rmsERadial[i] = norme*std::sqrt(std::abs(NbOfEvents*fSumE2Radial[i]
243  - fSumERadial[i]*fSumERadial[i]));
244 
245  MeanERadialCumul[i] = norme*fSumERadialCumul[i];
246  rmsERadialCumul[i] =
247  norme*std::sqrt(std::abs(NbOfEvents*fSumE2RadialCumul[i]
248  - fSumERadialCumul[i]*fSumERadialCumul[i]));
249 
250  G4double bin = (i+0.5)*dRradl;
251  analysisManager->FillH1(8, bin,MeanERadial[i]/dRradl);
252  analysisManager->FillH1(9, bin, rmsERadial[i]/dRradl);
253  bin = (i+1)*dRradl;
254  analysisManager->FillH1(10, bin,MeanERadialCumul[i]);
255  analysisManager->FillH1(11, bin, rmsERadialCumul[i]);
256  }
257 
258  //find Moliere confinement
259  //
260  const G4double EMoliere = 90.;
261  G4double iMoliere = 0.;
262  if ((MeanERadialCumul[0] <= EMoliere) &&
263  (MeanERadialCumul[f_nRbin-1] >= EMoliere)) {
264  G4int imin = 0;
265  while( (imin < f_nRbin-1) && (MeanERadialCumul[imin] < EMoliere) )
266  { ++imin; }
267  G4double ratio = (EMoliere - MeanERadialCumul[imin]) /
268  (MeanERadialCumul[imin+1] - MeanERadialCumul[imin]);
269  iMoliere = 1. + imin + ratio;
270  }
271 
272  //track length
273  //
274  norme = 1./(NbOfEvents*(fDet->GetMaterial()->GetRadlen()));
275  G4double MeanChargTrLength = norme*fSumChargTrLength;
276  G4double rmsChargTrLength =
277  norme*std::sqrt(std::fabs(NbOfEvents*fSum2ChargTrLength
278  - fSumChargTrLength*fSumChargTrLength));
279 
280  G4double MeanNeutrTrLength = norme*fSumNeutrTrLength;
281  G4double rmsNeutrTrLength =
282  norme*std::sqrt(std::fabs(NbOfEvents*fSum2NeutrTrLength
283  - fSumNeutrTrLength*fSumNeutrTrLength));
284 
285  //print
286  std::ios::fmtflags mode = G4cout.flags();
287  G4cout.setf(std::ios::fixed,std::ios::floatfield);
288  G4int prec = G4cout.precision(2);
289 
290  if (fVerbose) {
291 
292  G4cout << " LOGITUDINAL PROFILE "
293  << " CUMULATIVE LOGITUDINAL PROFILE" << G4endl << G4endl;
294 
295  G4cout << " bin " << " Mean rms "
296  << " bin " << " Mean rms \n" << G4endl;
297 
298  for (i=0; i<f_nLbin; i++) {
299  G4double inf=i*dLradl, sup=inf+dLradl;
300 
301  G4cout << std::setw(8) << inf << "->"
302  << std::setw(5) << sup << " radl: "
303  << std::setw(7) << MeanELongit[i] << "% "
304  << std::setw(9) << rmsELongit[i] << "% "
305  << " 0->" << std::setw(5) << sup << " radl: "
306  << std::setw(7) << MeanELongitCumul[i] << "% "
307  << std::setw(7) << rmsELongitCumul[i] << "% "
308  <<G4endl;
309  }
310 
311  G4cout << G4endl << G4endl << G4endl;
312 
313  G4cout << " RADIAL PROFILE "
314  << " CUMULATIVE RADIAL PROFILE" << G4endl << G4endl;
315 
316  G4cout << " bin " << " Mean rms "
317  << " bin " << " Mean rms \n" << G4endl;
318 
319  for (i=0; i<f_nRbin; i++) {
320  G4double inf=i*dRradl, sup=inf+dRradl;
321 
322  G4cout << std::setw(8) << inf << "->"
323  << std::setw(5) << sup << " radl: "
324  << std::setw(7) << MeanERadial[i] << "% "
325  << std::setw(9) << rmsERadial[i] << "% "
326  << " 0->" << std::setw(5) << sup << " radl: "
327  << std::setw(7) << MeanERadialCumul[i] << "% "
328  << std::setw(7) << rmsERadialCumul[i] << "% "
329  <<G4endl;
330  }
331  }
332 
333  G4cout << "\n ===== SUMMARY ===== \n" << G4endl;
334 
335  G4cout << " Total number pf events: " << NbOfEvents << "\n"
336  << " Mean number of charged steps: " << fChargedStep << G4endl;
337  G4cout << " Mean number of neutral steps: " << fNeutralStep
338  << "\n" << G4endl;
339 
340  G4cout << " energy deposit : "
341  << std::setw(7) << MeanELongitCumul[f_nLbin-1] << " % E0 +- "
342  << std::setw(7) << rmsELongitCumul[f_nLbin-1] << " % E0" << G4endl;
343  G4cout << " charged traklen: "
344  << std::setw(7) << MeanChargTrLength << " radl +- "
345  << std::setw(7) << rmsChargTrLength << " radl" << G4endl;
346  G4cout << " neutral traklen: "
347  << std::setw(7) << MeanNeutrTrLength << " radl +- "
348  << std::setw(7) << rmsNeutrTrLength << " radl" << G4endl;
349 
350  if (iMoliere > 0. ) {
351  G4double RMoliere1 = iMoliere*fDet->GetdRradl();
352  G4double RMoliere2 = iMoliere*fDet->GetdRlength();
353  G4cout << "\n " << EMoliere << " % confinement: radius = "
354  << RMoliere1 << " radl ("
355  << G4BestUnit( RMoliere2, "Length") << ")" << "\n" << G4endl;
356  }
357 
358  G4cout.setf(mode,std::ios::floatfield);
359  G4cout.precision(prec);
360 
361  // Acceptance
362 
363  G4int nLbin = fDet->GetnLtot();
364  if (limit < DBL_MAX) {
365  EmAcceptance acc;
366  acc.BeginOfAcceptance("Total Energy in Absorber",NbOfEvents);
367  G4double e = MeanELongitCumul[nLbin-1]/100.;
368  G4double r = rmsELongitCumul[nLbin-1]/100.;
369  acc.EmAcceptanceGauss("Edep",NbOfEvents,e,edep,rms,limit);
370  acc.EmAcceptanceGauss("Erms",NbOfEvents,r,rms,rms,2.0*limit);
371  acc.EndOfAcceptance();
372  }
373  limit = DBL_MAX;
374 }
375 
376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
virtual void Merge(const G4Run *)
Definition: G4Run.cc:54
void EmAcceptanceGauss(const G4String &title, G4int stat, G4double avr, G4double avr0, G4double rms, G4double limit)
Definition: EmAcceptance.cc:68
#define assert(x)
Definition: mymalloc.cc:1309
float bin[41]
Definition: plottest35.C:14
void InitializePerEvent()
Definition: Run.cc:103
void BeginOfAcceptance(const G4String &title, G4int stat)
Definition: EmAcceptance.cc:49
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
int G4int
Definition: G4Types.hh:78
std::vector< G4double > MyVector
Definition: Run.hh:46
Double_t edep
Definition: macro.C:13
G4GLOB_DLL std::ostream G4cout
void EndOfAcceptance()
Definition: EmAcceptance.cc:58
G4int GetNumberOfEvent() const
Definition: G4Run.hh:79
const G4int MaxBin
Definition: G4Run.hh:46
ExG4HbookAnalysisManager G4AnalysisManager
Definition: g4hbook_defs.hh:46
G4double GetRadlen() const
Definition: G4Material.hh:218
jump r
Definition: plot.C:36
G4double GetPDGMass() const
void ComputeStatistics()
Definition: Run.cc:97
G4ParticleGun * GetParticleGun()
G4ParticleDefinition * GetParticleDefinition() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
Definition: Run.hh:46
void FillPerEvent()
Definition: Run.cc:118
#define DBL_MAX
Definition: templates.hh:83
virtual void Merge(const G4Run *)
Definition: Run.cc:180
G4double GetParticleEnergy() const
~Run()
Definition: Run.cc:60
Run(DetectorConstruction *)
Definition: Run.cc:46