| Geant4 Cross Reference |
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 //
26 /// \file electromagnetic/TestEm8/src/Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32
33 #include "Run.hh"
34
35 #include "TestParameters.hh"
36
37 #include "G4ElectronIonPair.hh"
38 #include "G4LossTableManager.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4Run.hh"
41 #include "G4Step.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "Randomize.hh"
44
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46
47 Run::Run() : G4Run(), fElIonPair(0), fParam(TestParameters::GetPointer())
48 {
49 fElIonPair = G4LossTableManager::Instance()->ElectronIonPair();
50 }
51
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53
54 void Run::BeginOfRun()
55 {
56 // initilise scoring
57 fTotStepGas = fTotCluster = fMeanCluster = fOverflow = fTotEdep = fStepGas = fCluster = 0.0;
58 fEvt = 0;
59
60 fFactorALICE = fParam->GetFactorALICE();
61 fWidthALICE = fParam->GetEnergySmear();
62
63 SetVerbose(1);
64
65 fNbins = fParam->GetNumberBins();
66 fMaxEnergy = fParam->GetMaxEnergy();
67
68 fEgas.resize(fNbins, 0.0);
69 fEdep.reset();
70
71 if (fVerbose > 0) {
72 G4int binsCluster = fParam->GetNumberBinsCluster();
73 G4cout << " BinsCluster= " << binsCluster << " BinsE= " << fNbins
74 << " Emax(keV)= " << fMaxEnergy / keV << G4endl;
75 G4cout << " WidthALICE(keV)= " << fWidthALICE / keV << " FactorALICE= " << fFactorALICE
76 << G4endl;
77 }
78 }
79
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81
82 void Run::EndOfRun()
83 {
84 G4int nEvt = GetNumberOfEvent();
85 G4double norm = (nEvt > 0) ? 1.0 / (G4double)nEvt : 0.0;
86
87 fTotStepGas *= norm;
88 fTotCluster *= norm;
89 fMeanCluster *= norm;
90 fOverflow *= norm;
91
92 G4double y1 = fEdep.mean();
93 G4double y2 = fEdep.rms();
94
95 G4double de = fMaxEnergy / G4double(fNbins);
96 G4double x1 = -de * 0.5;
97
98 fFactorALICE = fParam->GetFactorALICE();
99
100 G4cout << " ====================================================" << G4endl;
101 G4cout << " Beam Particle: " << fParam->GetBeamParticle()->GetParticleName() << G4endl
102 << " Ekin(MeV) = " << fParam->GetBeamEnergy() / MeV << G4endl
103 << " Z(mm) = " << fParam->GetPositionZ() / mm << G4endl;
104 G4cout << " ================== run summary =====================" << G4endl;
105 G4int prec = G4cout.precision(5);
106 G4cout << " End of Run TotNbofEvents = " << nEvt << G4endl;
107 G4cout << " Energy(keV) per ADC channel = " << 1.0 / (keV * fFactorALICE) << G4endl;
108
109 G4cout << G4endl;
110 G4cout << " Mean energy deposit in absorber = " << y1 / keV << " +- "
111 << y2 * std::sqrt(norm) / keV << " keV; ";
112 if (y1 > 0.0) {
113 G4cout << " RMS/Emean = " << y2 / y1;
114 }
115 G4cout << G4endl;
116 G4cout << " Mean number of steps in absorber= " << fTotStepGas
117 << "; mean number of ion-clusters = " << fTotCluster << " MeanCluster= " << fMeanCluster
118 << G4endl;
119 G4cout << G4endl;
120
121 G4cout << " ====== Energy deposit distribution Noverflows= " << fOverflow
122 << " ====== " << G4endl;
123 G4cout << " bin nb Elow entries normalized " << G4endl;
124
125 std::ofstream fileOut("distribution.out", std::ios::out);
126 fileOut.setf(std::ios::scientific, std::ios::floatfield);
127
128 x1 = 0.0;
129
130 fileOut << fNbins << G4endl;
131
132 for (G4int j = 0; j < fNbins; ++j) {
133 G4cout << std::setw(5) << j << std::setw(10) << x1 / keV << std::setw(12) << fEgas[j]
134 << std::setw(12) << fEgas[j] * norm << G4endl;
135 fileOut << x1 / keV << "\t" << fEgas[j] << G4endl;
136 x1 += de;
137 }
138 G4cout.precision(prec);
139
140 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
141 // normalize histograms
142 G4double normf = fParam->GetNormFactor();
143 analysisManager->ScaleH1(1, norm);
144 analysisManager->ScaleH1(2, norm);
145 analysisManager->ScaleH1(3, norm * normf);
146
147 G4cout << " ================== run end ==========================" << G4endl;
148 }
149
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
151
152 void Run::BeginOfEvent()
153 {
154 fTotEdep = 0.0;
155 fStepGas = 0;
156 fCluster = 0;
157 ++fEvt;
158 }
159
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
161
162 void Run::EndOfEvent()
163 {
164 fTotStepGas += fStepGas;
165 fTotCluster += fCluster;
166
167 if (fWidthALICE > 0.0) {
168 G4double x = G4RandGauss::shoot(0., fWidthALICE);
169 fTotEdep += x;
170 fTotEdep = std::max(fTotEdep, 0.0);
171 }
172
173 G4int idx = G4int(fTotEdep * fNbins / fMaxEnergy);
174
175 if (idx < 0) {
176 fEgas[0] += 1.0;
177 }
178 if (idx >= fNbins) {
179 fOverflow += 1.0;
180 }
181 else {
182 fEgas[idx] += 1.0;
183 }
184
185 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
186 // fill histo
187 analysisManager->FillH1(1, fTotEdep / keV, 1.0);
188 analysisManager->FillH1(2, fCluster, 1.0);
189 analysisManager->FillH1(3, fTotEdep * fFactorALICE, 1.0);
190 fEdep.fill(fTotEdep, 1.0);
191 }
192
193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
194
195 void Run::Merge(const G4Run* run)
196 {
197 const Run* localRun = static_cast<const Run*>(run);
198
199 fTotStepGas += localRun->fTotStepGas;
200 fTotCluster += localRun->fTotCluster;
201 fMeanCluster += localRun->fMeanCluster;
202 fOverflow += localRun->fOverflow;
203
204 G4StatDouble* stat = const_cast<G4StatDouble*>(localRun->GetStat());
205
206 fEdep.add(stat);
207
208 for (G4int j = 0; j < fNbins; ++j) {
209 fEgas[j] += localRun->fEgas[j];
210 }
211
212 G4Run::Merge(run);
213 }
214
215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
216
217 void Run::AddEnergy(G4double edep, const G4Step* step)
218 {
219 if (1 < fVerbose) {
220 G4cout << "Run::AddEnergy: e(keV)= " << edep / keV << G4endl;
221 }
222 fTotEdep += edep;
223 if (step) {
224 if (1 == step->GetTrack()->GetTrackID()) {
225 fStepGas += 1.0;
226 }
227
228 fMeanCluster += fElIonPair->MeanNumberOfIonsAlongStep(step);
229 fCluster += fElIonPair->SampleNumberOfIonsAlongStep(step);
230 }
231 }
232
233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
234