Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 "EmAcceptance.hh"
39 
40 #include "G4Run.hh"
41 #include "G4UnitsTable.hh"
42 #include "G4SystemOfUnits.hh"
43 #include <iomanip>
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46 
48  :G4Run(),fDet(det),fKin(kin),
49  f_nLbin(kMaxBin),f_nRbin(kMaxBin)
50 {
51  Reset();
52 }
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 
56 void Run::Reset()
57 {
58  f_nLbin = fDet->GetnLtot();
59  f_dEdL.resize(f_nLbin);
60  fSumELongit.resize(f_nLbin);
61  fSumELongitCumul.resize(f_nLbin);
62  fSumE2Longit.resize(f_nLbin);
63  fSumE2LongitCumul.resize(f_nLbin);
64 
65  f_nRbin = fDet->GetnRtot();
66  f_dEdR.resize(f_nRbin);
67  fSumERadial.resize(f_nRbin);
68  fSumERadialCumul.resize(f_nRbin);
69  fSumE2Radial.resize(f_nRbin);
70  fSumE2RadialCumul.resize(f_nRbin);
71 
72  fChargedStep = 0.0;
73  fNeutralStep = 0.0;
74 
75  fVerbose = 0;
76 
77  //initialize arrays of cumulative energy deposition
78  //
79  for (G4int i=0; i<f_nLbin; i++)
80  fSumELongit[i]=fSumE2Longit[i]=fSumELongitCumul[i]=fSumE2LongitCumul[i]=0.;
81 
82  for (G4int j=0; j<f_nRbin; j++)
83  fSumERadial[j]=fSumE2Radial[j]=fSumERadialCumul[j]=fSumE2RadialCumul[j]=0.;
84 
85  //initialize track length
86  fSumChargTrLength=fSum2ChargTrLength=fSumNeutrTrLength=fSum2NeutrTrLength=0.;
87 
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 
92 Run::~Run()
93 {}
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
96 
98 {
99  //initialize arrays of energy deposit per bin
100  for (G4int i=0; i<f_nLbin; i++)
101  { f_dEdL[i] = 0.; }
102 
103  for (G4int j=0; j<f_nRbin; j++)
104  { f_dEdR[j] = 0.; }
105 
106  //initialize tracklength
107  fChargTrLength = fNeutrTrLength = 0.;
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111 
113 {
114  //accumulate statistic
115  //
116  G4double dLCumul = 0.;
117  for (G4int i=0; i<f_nLbin; i++)
118  {
119  fSumELongit[i] += f_dEdL[i];
120  fSumE2Longit[i] += f_dEdL[i]*f_dEdL[i];
121  dLCumul += f_dEdL[i];
122  fSumELongitCumul[i] += dLCumul;
123  fSumE2LongitCumul[i] += dLCumul*dLCumul;
124  }
125 
126  G4double dRCumul = 0.;
127  for (G4int j=0; j<f_nRbin; j++)
128  {
129  fSumERadial[j] += f_dEdR[j];
130  fSumE2Radial[j] += f_dEdR[j]*f_dEdR[j];
131  dRCumul += f_dEdR[j];
132  fSumERadialCumul[j] += dRCumul;
133  fSumE2RadialCumul[j] += dRCumul*dRCumul;
134  }
135 
136  fSumChargTrLength += fChargTrLength;
137  fSum2ChargTrLength += fChargTrLength*fChargTrLength;
138  fSumNeutrTrLength += fNeutrTrLength;
139  fSum2NeutrTrLength += fNeutrTrLength*fNeutrTrLength;
140 
141  //fill histograms
142  //
143 
144  G4double Ekin=fKin->GetParticleGun()->GetParticleEnergy();
146  G4double radl=fDet->GetMaterial()->GetRadlen();
147 
148  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
149  analysisManager->FillH1(1, 100.*dLCumul/(Ekin+mass));
150  analysisManager->FillH1(2, fChargTrLength/radl);
151  analysisManager->FillH1(3, fNeutrTrLength/radl);
152 
153  //profiles
154  G4double norm = 100./(Ekin+mass);
155  G4double dLradl = fDet->GetdLradl();
156  for (G4int i=0; i<f_nLbin; i++) {
157  G4double bin = (i+0.5)*dLradl;
158  analysisManager->FillP1(0, bin, norm*f_dEdL[i]/dLradl);
159  }
160  G4double dRradl = fDet->GetdRradl();
161  for (G4int j=0; j<f_nRbin; j++) {
162  G4double bin = (j+0.5)*dRradl;
163  analysisManager->FillP1(1, bin, norm*f_dEdR[j]/dRradl);
164  }
165 }
166 
167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
168 
169 void Run::Merge(const G4Run* run)
170 {
171  const Run* localRun = static_cast<const Run*>(run);
172 
173  fChargedStep += localRun->fChargedStep;
174  fNeutralStep += localRun->fNeutralStep;
175 
176  for (G4int i=0; i<f_nLbin; i++) {
177  fSumELongit[i] += localRun->fSumELongit[i];
178  fSumE2Longit[i] += localRun->fSumE2Longit[i];
179  fSumELongitCumul[i] += localRun->fSumELongitCumul[i];
180  fSumE2LongitCumul[i] += localRun->fSumE2LongitCumul[i];
181  }
182 
183  for (G4int j=0; j<f_nRbin; j++) {
184  fSumERadial[j] += localRun->fSumERadial[j];
185  fSumE2Radial[j] += localRun->fSumE2Radial[j];
186  fSumERadialCumul[j] += localRun->fSumERadialCumul[j];
187  fSumE2RadialCumul[j] += localRun->fSumE2RadialCumul[j];
188  }
189 
190  fSumChargTrLength += localRun->fSumChargTrLength;
191  fSum2ChargTrLength += localRun->fSum2ChargTrLength;
192  fSumNeutrTrLength += localRun->fSumNeutrTrLength;
193  fSum2NeutrTrLength += localRun->fSum2NeutrTrLength;
194 
195  G4Run::Merge(run);
196 }
197 
198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
199 
200 void Run::EndOfRun(G4double edep, G4double rms, G4double& limit)
201 {
202  G4int NbOfEvents = GetNumberOfEvent();
203 
204  G4double kinEnergy = fKin->GetParticleGun()->GetParticleEnergy();
205  assert(NbOfEvents*kinEnergy > 0);
206 
207  fChargedStep /= G4double(NbOfEvents);
208  fNeutralStep /= G4double(NbOfEvents);
209 
211  G4double norme = 100./(NbOfEvents*(kinEnergy+mass));
212 
213  //longitudinal
214  //
215  G4double dLradl = fDet->GetdLradl();
216 
217  MyVector MeanELongit(f_nLbin), rmsELongit(f_nLbin);
218  MyVector MeanELongitCumul(f_nLbin), rmsELongitCumul(f_nLbin);
219 
220  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
221 
222  G4int i;
223  for (i=0; i<f_nLbin; i++) {
224  MeanELongit[i] = norme*fSumELongit[i];
225  rmsELongit[i] =
226  norme*std::sqrt(std::abs(NbOfEvents*fSumE2Longit[i]
227  - fSumELongit[i]*fSumELongit[i]));
228 
229  MeanELongitCumul[i] = norme*fSumELongitCumul[i];
230  rmsELongitCumul[i] = norme*std::sqrt(std::abs(NbOfEvents*
231  fSumE2LongitCumul[i] - fSumELongitCumul[i]*fSumELongitCumul[i]));
232  G4double bin = (i+0.5)*dLradl;
233  analysisManager->FillH1(4, bin,MeanELongit[i]/dLradl);
234  analysisManager->FillH1(5, bin, rmsELongit[i]/dLradl);
235  bin = (i+1)*dLradl;
236  analysisManager->FillH1(6, bin,MeanELongitCumul[i]);
237  analysisManager->FillH1(7, bin, rmsELongitCumul[i]);
238  }
239 
240  //radial
241  //
242  G4double dRradl = fDet->GetdRradl();
243 
244  MyVector MeanERadial(f_nRbin), rmsERadial(f_nRbin);
245  MyVector MeanERadialCumul(f_nRbin), rmsERadialCumul(f_nRbin);
246 
247  for (i=0; i<f_nRbin; i++) {
248  MeanERadial[i] = norme*fSumERadial[i];
249  rmsERadial[i] = norme*std::sqrt(std::abs(NbOfEvents*fSumE2Radial[i]
250  - fSumERadial[i]*fSumERadial[i]));
251 
252  MeanERadialCumul[i] = norme*fSumERadialCumul[i];
253  rmsERadialCumul[i] =
254  norme*std::sqrt(std::abs(NbOfEvents*fSumE2RadialCumul[i]
255  - fSumERadialCumul[i]*fSumERadialCumul[i]));
256 
257  G4double bin = (i+0.5)*dRradl;
258  analysisManager->FillH1(8, bin,MeanERadial[i]/dRradl);
259  analysisManager->FillH1(9, bin, rmsERadial[i]/dRradl);
260  bin = (i+1)*dRradl;
261  analysisManager->FillH1(10, bin,MeanERadialCumul[i]);
262  analysisManager->FillH1(11, bin, rmsERadialCumul[i]);
263  }
264 
265  //find Moliere confinement
266  //
267  const G4double EMoliere = 90.;
268  G4double iMoliere = 0.;
269  if ((MeanERadialCumul[0] <= EMoliere) &&
270  (MeanERadialCumul[f_nRbin-1] >= EMoliere)) {
271  G4int imin = 0;
272  while( (imin < f_nRbin-1) && (MeanERadialCumul[imin] < EMoliere) )
273  { ++imin; }
274  G4double ratio = (EMoliere - MeanERadialCumul[imin]) /
275  (MeanERadialCumul[imin+1] - MeanERadialCumul[imin]);
276  iMoliere = 1. + imin + ratio;
277  }
278 
279  //track length
280  //
281  norme = 1./(NbOfEvents*(fDet->GetMaterial()->GetRadlen()));
282  G4double MeanChargTrLength = norme*fSumChargTrLength;
283  G4double rmsChargTrLength =
284  norme*std::sqrt(std::fabs(NbOfEvents*fSum2ChargTrLength
285  - fSumChargTrLength*fSumChargTrLength));
286 
287  G4double MeanNeutrTrLength = norme*fSumNeutrTrLength;
288  G4double rmsNeutrTrLength =
289  norme*std::sqrt(std::fabs(NbOfEvents*fSum2NeutrTrLength
290  - fSumNeutrTrLength*fSumNeutrTrLength));
291 
292  //print
293  std::ios::fmtflags mode = G4cout.flags();
294  G4cout.setf(std::ios::fixed,std::ios::floatfield);
295  G4int prec = G4cout.precision(2);
296 
297  if (fVerbose) {
298 
299  G4cout << " LOGITUDINAL PROFILE "
300  << " CUMULATIVE LOGITUDINAL PROFILE" << G4endl << G4endl;
301 
302  G4cout << " bin " << " Mean rms "
303  << " bin " << " Mean rms \n" << G4endl;
304 
305  for (i=0; i<f_nLbin; i++) {
306  G4double inf=i*dLradl, sup=inf+dLradl;
307 
308  G4cout << std::setw(8) << inf << "->"
309  << std::setw(5) << sup << " radl: "
310  << std::setw(7) << MeanELongit[i] << "% "
311  << std::setw(9) << rmsELongit[i] << "% "
312  << " 0->" << std::setw(5) << sup << " radl: "
313  << std::setw(7) << MeanELongitCumul[i] << "% "
314  << std::setw(7) << rmsELongitCumul[i] << "% "
315  <<G4endl;
316  }
317 
318  G4cout << G4endl << G4endl << G4endl;
319 
320  G4cout << " RADIAL PROFILE "
321  << " CUMULATIVE RADIAL PROFILE" << G4endl << G4endl;
322 
323  G4cout << " bin " << " Mean rms "
324  << " bin " << " Mean rms \n" << G4endl;
325 
326  for (i=0; i<f_nRbin; i++) {
327  G4double inf=i*dRradl, sup=inf+dRradl;
328 
329  G4cout << std::setw(8) << inf << "->"
330  << std::setw(5) << sup << " radl: "
331  << std::setw(7) << MeanERadial[i] << "% "
332  << std::setw(9) << rmsERadial[i] << "% "
333  << " 0->" << std::setw(5) << sup << " radl: "
334  << std::setw(7) << MeanERadialCumul[i] << "% "
335  << std::setw(7) << rmsERadialCumul[i] << "% "
336  <<G4endl;
337  }
338  }
339 
340  G4cout << "\n ===== SUMMARY ===== \n" << G4endl;
341 
342  G4cout << " Total number of events: " << NbOfEvents << "\n"
343  << " Mean number of charged steps: " << fChargedStep << G4endl;
344  G4cout << " Mean number of neutral steps: " << fNeutralStep
345  << "\n" << G4endl;
346 
347  G4cout << " energy deposit : "
348  << std::setw(7) << MeanELongitCumul[f_nLbin-1] << " % E0 +- "
349  << std::setw(7) << rmsELongitCumul[f_nLbin-1] << " % E0" << G4endl;
350  G4cout << " charged traklen: "
351  << std::setw(7) << MeanChargTrLength << " radl +- "
352  << std::setw(7) << rmsChargTrLength << " radl" << G4endl;
353  G4cout << " neutral traklen: "
354  << std::setw(7) << MeanNeutrTrLength << " radl +- "
355  << std::setw(7) << rmsNeutrTrLength << " radl" << G4endl;
356 
357  if (iMoliere > 0. ) {
358  G4double RMoliere1 = iMoliere*fDet->GetdRradl();
359  G4double RMoliere2 = iMoliere*fDet->GetdRlength();
360  G4cout << "\n " << EMoliere << " % confinement: radius = "
361  << RMoliere1 << " radl ("
362  << G4BestUnit( RMoliere2, "Length") << ")" << "\n" << G4endl;
363  }
364 
365  G4cout.setf(mode,std::ios::floatfield);
366  G4cout.precision(prec);
367 
368  // Acceptance
369 
370  G4int nLbin = fDet->GetnLtot();
371  if (limit < DBL_MAX) {
372  EmAcceptance acc;
373  acc.BeginOfAcceptance("Total Energy in Absorber",NbOfEvents);
374  G4double e = MeanELongitCumul[nLbin-1]/100.;
375  G4double r = rmsELongitCumul[nLbin-1]/100.;
376  acc.EmAcceptanceGauss("Edep",NbOfEvents,e,edep,rms,limit);
377  acc.EmAcceptanceGauss("Erms",NbOfEvents,r,rms,rms,2.0*limit);
378  acc.EndOfAcceptance();
379  }
380  limit = DBL_MAX;
381 }
382 
383 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
The primary generator action class with particle gun.
virtual void Merge(const G4Run *)
Definition: G4Run.cc:54
tuple bin
Definition: plottest35.py:22
Run()
Definition: Run.cc:43
void EmAcceptanceGauss(const G4String &title, G4int stat, G4double avr, G4double avr0, G4double rms, G4double limit)
Definition: EmAcceptance.cc:68
void InitializePerEvent()
Definition: Run.cc:97
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
static const double prec
Definition: RanecuEngine.cc:58
G4GLOB_DLL std::ostream G4cout
void EndOfAcceptance()
Definition: EmAcceptance.cc:58
G4int GetNumberOfEvent() const
Definition: G4Run.hh:79
Definition: G4Run.hh:46
std::vector< G4double > MyVector
Definition: Run.hh:42
G4double GetRadlen() const
Definition: G4Material.hh:220
G4double GetPDGMass() const
G4ParticleGun * GetParticleGun()
G4ParticleDefinition * GetParticleDefinition() const
Detector construction class to define materials and geometry.
void EndOfRun()
Definition: Run.cc:147
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
Definition: Run.hh:46
void FillPerEvent()
Definition: Run.cc:112
#define DBL_MAX
Definition: templates.hh:83
virtual void Merge(const G4Run *)
Definition: Run.cc:115
G4CsvAnalysisManager G4AnalysisManager
Definition: g4csv_defs.hh:77
const G4int kMaxBin
G4double GetParticleEnergy() const
~Run()
Definition: Run.cc:72