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