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