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