| 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/TestEm5/src/Run.cc 26 /// \file electromagnetic/TestEm5/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 "DetectorConstruction.hh" 34 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh" <<
37 #include "PrimaryGeneratorAction.hh" 35 #include "PrimaryGeneratorAction.hh"
>> 36 #include "HistoManager.hh"
38 37
39 #include "G4EmCalculator.hh" 38 #include "G4EmCalculator.hh"
40 #include "G4SystemOfUnits.hh" 39 #include "G4SystemOfUnits.hh"
41 #include "G4UnitsTable.hh" 40 #include "G4UnitsTable.hh"
42 41
43 #include <iomanip> 42 #include <iomanip>
44 43
45 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46 45
47 Run::Run(DetectorConstruction* det) : fDetecto << 46 Run::Run(DetectorConstruction* det)
>> 47 : G4Run(),
>> 48 fDetector(det),
>> 49 fParticle(0), fEkin(0.)
>> 50 {
>> 51 fEnergyDeposit = fEnergyDeposit2 = 0.;
>> 52 fTrakLenCharged = fTrakLenCharged2 = 0.;
>> 53 fTrakLenNeutral = fTrakLenNeutral2 = 0.;
>> 54 fNbStepsCharged = fNbStepsCharged2 = 0.;
>> 55 fNbStepsNeutral = fNbStepsNeutral2 = 0.;
>> 56 fMscProjecTheta = fMscProjecTheta2 = 0.;
>> 57 fMscThetaCentral = 0.;
>> 58
>> 59 fNbGamma = fNbElect = fNbPosit = 0;
>> 60
>> 61 fTransmit[0] = fTransmit[1] = fReflect[0] = fReflect[1] = 0;
>> 62
>> 63 fMscEntryCentral = 0;
>> 64 fTypes[0] = fTypes[1] = fTypes[2] = fTypes[3] = 0;
>> 65
>> 66 fEnergyLeak[0] = fEnergyLeak[1] = fEnergyLeak2[0] = fEnergyLeak2[1] = 0.;
>> 67 }
48 68
49 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 70
51 void Run::SetPrimary(G4ParticleDefinition* par 71 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
52 { << 72 {
53 fParticle = particle; 73 fParticle = particle;
54 fEkin = energy; 74 fEkin = energy;
55 << 75
56 // compute theta0 << 76 //compute theta0
57 fMscThetaCentral = 3 * ComputeMscHighland(); << 77 fMscThetaCentral = 3*ComputeMscHighland();
58 } 78 }
59 << 79
60 //....oooOO0OOooo........oooOO0OOooo........oo 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 81
62 void Run::Merge(const G4Run* run) 82 void Run::Merge(const G4Run* run)
63 { 83 {
64 const Run* localRun = static_cast<const Run* 84 const Run* localRun = static_cast<const Run*>(run);
65 85
66 // pass information about primary particle 86 // pass information about primary particle
67 fParticle = localRun->fParticle; 87 fParticle = localRun->fParticle;
68 fEkin = localRun->fEkin; << 88 fEkin = localRun->fEkin;
69 << 89
70 fMscThetaCentral = localRun->fMscThetaCentra 90 fMscThetaCentral = localRun->fMscThetaCentral;
71 91
72 // accumulate sums 92 // accumulate sums
73 // 93 //
74 fEnergyDeposit += localRun->fEnergyDeposit; << 94 fEnergyDeposit += localRun->fEnergyDeposit;
75 fEnergyDeposit2 += localRun->fEnergyDeposit2 << 95 fEnergyDeposit2 += localRun->fEnergyDeposit2;
76 fTrakLenCharged += localRun->fTrakLenCharged << 96 fTrakLenCharged += localRun->fTrakLenCharged;
77 fTrakLenCharged2 += localRun->fTrakLenCharge << 97 fTrakLenCharged2 += localRun->fTrakLenCharged2;
78 fTrakLenNeutral += localRun->fTrakLenNeutral << 98 fTrakLenNeutral += localRun->fTrakLenNeutral;
79 fTrakLenNeutral2 += localRun->fTrakLenNeutra 99 fTrakLenNeutral2 += localRun->fTrakLenNeutral2;
80 fNbStepsCharged += localRun->fNbStepsCharged << 100 fNbStepsCharged += localRun->fNbStepsCharged;
81 fNbStepsCharged2 += localRun->fNbStepsCharge 101 fNbStepsCharged2 += localRun->fNbStepsCharged2;
82 fNbStepsNeutral += localRun->fNbStepsNeutral << 102 fNbStepsNeutral += localRun->fNbStepsNeutral;
83 fNbStepsNeutral2 += localRun->fNbStepsNeutra 103 fNbStepsNeutral2 += localRun->fNbStepsNeutral2;
84 fMscProjecTheta += localRun->fMscProjecTheta << 104 fMscProjecTheta += localRun->fMscProjecTheta;
85 fMscProjecTheta2 += localRun->fMscProjecThet 105 fMscProjecTheta2 += localRun->fMscProjecTheta2;
86 106
87 fTypes[0] += localRun->fTypes[0]; 107 fTypes[0] += localRun->fTypes[0];
88 fTypes[1] += localRun->fTypes[1]; 108 fTypes[1] += localRun->fTypes[1];
89 fTypes[2] += localRun->fTypes[2]; 109 fTypes[2] += localRun->fTypes[2];
90 fTypes[3] += localRun->fTypes[3]; 110 fTypes[3] += localRun->fTypes[3];
91 << 111
92 fNbGamma += localRun->fNbGamma; 112 fNbGamma += localRun->fNbGamma;
93 fNbElect += localRun->fNbElect; << 113 fNbElect += localRun->fNbElect;
94 fNbPosit += localRun->fNbPosit; 114 fNbPosit += localRun->fNbPosit;
95 << 115
96 fTransmit[0] += localRun->fTransmit[0]; << 116 fTransmit[0] += localRun->fTransmit[0];
97 fTransmit[1] += localRun->fTransmit[1]; 117 fTransmit[1] += localRun->fTransmit[1];
98 fReflect[0] += localRun->fReflect[0]; << 118 fReflect[0] += localRun->fReflect[0];
99 fReflect[1] += localRun->fReflect[1]; << 119 fReflect[1] += localRun->fReflect[1];
100 << 120
101 fMscEntryCentral += localRun->fMscEntryCentr 121 fMscEntryCentral += localRun->fMscEntryCentral;
102 << 122
103 fEnergyLeak[0] += localRun->fEnergyLeak[0]; << 123 fEnergyLeak[0] += localRun->fEnergyLeak[0];
104 fEnergyLeak[1] += localRun->fEnergyLeak[1]; << 124 fEnergyLeak[1] += localRun->fEnergyLeak[1];
105 fEnergyLeak2[0] += localRun->fEnergyLeak2[0] 125 fEnergyLeak2[0] += localRun->fEnergyLeak2[0];
106 fEnergyLeak2[1] += localRun->fEnergyLeak2[1] 126 fEnergyLeak2[1] += localRun->fEnergyLeak2[1];
107 << 127
108 G4Run::Merge(run); << 128 G4Run::Merge(run);
109 } << 129 }
110 130
111 //....oooOO0OOooo........oooOO0OOooo........oo 131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
112 132
113 void Run::EndOfRun() 133 void Run::EndOfRun()
114 { 134 {
115 // compute mean and rms 135 // compute mean and rms
116 // 136 //
117 G4int TotNbofEvents = numberOfEvent; 137 G4int TotNbofEvents = numberOfEvent;
118 if (TotNbofEvents == 0) return; 138 if (TotNbofEvents == 0) return;
119 << 139
120 G4double EnergyBalance = fEnergyDeposit + fE 140 G4double EnergyBalance = fEnergyDeposit + fEnergyLeak[0] + fEnergyLeak[1];
121 EnergyBalance /= TotNbofEvents; 141 EnergyBalance /= TotNbofEvents;
122 142
123 fEnergyDeposit /= TotNbofEvents; << 143 fEnergyDeposit /= TotNbofEvents; fEnergyDeposit2 /= TotNbofEvents;
124 fEnergyDeposit2 /= TotNbofEvents; << 144 G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit*fEnergyDeposit;
125 G4double rmsEdep = fEnergyDeposit2 - fEnergy << 145 if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep/TotNbofEvents);
126 if (rmsEdep > 0.) << 146 else rmsEdep = 0.;
127 rmsEdep = std::sqrt(rmsEdep / TotNbofEvent << 147
128 else << 148 fTrakLenCharged /= TotNbofEvents; fTrakLenCharged2 /= TotNbofEvents;
129 rmsEdep = 0.; << 149 G4double rmsTLCh = fTrakLenCharged2 - fTrakLenCharged*fTrakLenCharged;
130 << 150 if (rmsTLCh>0.) rmsTLCh = std::sqrt(rmsTLCh/TotNbofEvents);
131 fTrakLenCharged /= TotNbofEvents; << 151 else rmsTLCh = 0.;
132 fTrakLenCharged2 /= TotNbofEvents; << 152
133 G4double rmsTLCh = fTrakLenCharged2 - fTrakL << 153 fTrakLenNeutral /= TotNbofEvents; fTrakLenNeutral2 /= TotNbofEvents;
134 if (rmsTLCh > 0.) << 154 G4double rmsTLNe = fTrakLenNeutral2 - fTrakLenNeutral*fTrakLenNeutral;
135 rmsTLCh = std::sqrt(rmsTLCh / TotNbofEvent << 155 if (rmsTLNe>0.) rmsTLNe = std::sqrt(rmsTLNe/TotNbofEvents);
136 else << 156 else rmsTLNe = 0.;
137 rmsTLCh = 0.; << 157
138 << 158 fNbStepsCharged /= TotNbofEvents; fNbStepsCharged2 /= TotNbofEvents;
139 fTrakLenNeutral /= TotNbofEvents; << 159 G4double rmsStCh = fNbStepsCharged2 - fNbStepsCharged*fNbStepsCharged;
140 fTrakLenNeutral2 /= TotNbofEvents; << 160 if (rmsStCh>0.) rmsStCh = std::sqrt(rmsStCh/TotNbofEvents);
141 G4double rmsTLNe = fTrakLenNeutral2 - fTrakL << 161 else rmsStCh = 0.;
142 if (rmsTLNe > 0.) << 162
143 rmsTLNe = std::sqrt(rmsTLNe / TotNbofEvent << 163 fNbStepsNeutral /= TotNbofEvents; fNbStepsNeutral2 /= TotNbofEvents;
144 else << 164 G4double rmsStNe = fNbStepsNeutral2 - fNbStepsNeutral*fNbStepsNeutral;
145 rmsTLNe = 0.; << 165 if (rmsStNe>0.) rmsStNe = std::sqrt(rmsStNe/TotNbofEvents);
146 << 166 else rmsStNe = 0.;
147 fNbStepsCharged /= TotNbofEvents; << 167
148 fNbStepsCharged2 /= TotNbofEvents; << 168 G4double Gamma = (G4double)fNbGamma/TotNbofEvents;
149 G4double rmsStCh = fNbStepsCharged2 - fNbSte << 169 G4double Elect = (G4double)fNbElect/TotNbofEvents;
150 if (rmsStCh > 0.) << 170 G4double Posit = (G4double)fNbPosit/TotNbofEvents;
151 rmsStCh = std::sqrt(rmsStCh / TotNbofEvent <<
152 else <<
153 rmsStCh = 0.; <<
154 <<
155 fNbStepsNeutral /= TotNbofEvents; <<
156 fNbStepsNeutral2 /= TotNbofEvents; <<
157 G4double rmsStNe = fNbStepsNeutral2 - fNbSte <<
158 if (rmsStNe > 0.) <<
159 rmsStNe = std::sqrt(rmsStNe / TotNbofEvent <<
160 else <<
161 rmsStNe = 0.; <<
162 <<
163 G4double Gamma = (G4double)fNbGamma / TotNbo <<
164 G4double Elect = (G4double)fNbElect / TotNbo <<
165 G4double Posit = (G4double)fNbPosit / TotNbo <<
166 171
167 G4double transmit[2]; 172 G4double transmit[2];
168 transmit[0] = 100. * fTransmit[0] / TotNbofE << 173 transmit[0] = 100.*fTransmit[0]/TotNbofEvents;
169 transmit[1] = 100. * fTransmit[1] / TotNbofE << 174 transmit[1] = 100.*fTransmit[1]/TotNbofEvents;
170 175
171 G4double reflect[2]; 176 G4double reflect[2];
172 reflect[0] = 100. * fReflect[0] / TotNbofEve << 177 reflect[0] = 100.*fReflect[0]/TotNbofEvents;
173 reflect[1] = 100. * fReflect[1] / TotNbofEve << 178 reflect[1] = 100.*fReflect[1]/TotNbofEvents;
174 179
175 G4double rmsMsc = 0., tailMsc = 0.; 180 G4double rmsMsc = 0., tailMsc = 0.;
176 if (fMscEntryCentral > 0) { 181 if (fMscEntryCentral > 0) {
177 fMscProjecTheta /= fMscEntryCentral; << 182 fMscProjecTheta /= fMscEntryCentral; fMscProjecTheta2 /= fMscEntryCentral;
178 fMscProjecTheta2 /= fMscEntryCentral; << 183 rmsMsc = fMscProjecTheta2 - fMscProjecTheta*fMscProjecTheta;
179 rmsMsc = fMscProjecTheta2 - fMscProjecThet << 184 if (rmsMsc > 0.) { rmsMsc = std::sqrt(rmsMsc); }
180 if (rmsMsc > 0.) { << 185 if(fTransmit[1] > 0.0) {
181 rmsMsc = std::sqrt(rmsMsc); << 186 tailMsc = 100.- (100.*fMscEntryCentral)/(2*fTransmit[1]);
182 } << 187 }
183 if (fTransmit[1] > 0.0) { <<
184 tailMsc = 100. - (100. * fMscEntryCentra <<
185 } <<
186 } 188 }
187 << 189
188 fEnergyLeak[0] /= TotNbofEvents; << 190 fEnergyLeak[0] /= TotNbofEvents; fEnergyLeak2[0] /= TotNbofEvents;
189 fEnergyLeak2[0] /= TotNbofEvents; << 191 G4double rmsEl0 = fEnergyLeak2[0] - fEnergyLeak[0]*fEnergyLeak[0];
190 G4double rmsEl0 = fEnergyLeak2[0] - fEnergyL << 192 if (rmsEl0>0.) rmsEl0 = std::sqrt(rmsEl0/TotNbofEvents);
191 if (rmsEl0 > 0.) << 193 else rmsEl0 = 0.;
192 rmsEl0 = std::sqrt(rmsEl0 / TotNbofEvents) << 194
193 else << 195 fEnergyLeak[1] /= TotNbofEvents; fEnergyLeak2[1] /= TotNbofEvents;
194 rmsEl0 = 0.; << 196 G4double rmsEl1 = fEnergyLeak2[1] - fEnergyLeak[1]*fEnergyLeak[1];
195 << 197 if (rmsEl1>0.) rmsEl1 = std::sqrt(rmsEl1/TotNbofEvents);
196 fEnergyLeak[1] /= TotNbofEvents; << 198 else rmsEl1 = 0.;
197 fEnergyLeak2[1] /= TotNbofEvents; << 199
198 G4double rmsEl1 = fEnergyLeak2[1] - fEnergyL << 200
199 if (rmsEl1 > 0.) << 201 //Stopping Power from input Table.
200 rmsEl1 = std::sqrt(rmsEl1 / TotNbofEvents) <<
201 else <<
202 rmsEl1 = 0.; <<
203 <<
204 // Stopping Power from input Table. <<
205 // 202 //
206 const G4Material* material = fDetector->GetA 203 const G4Material* material = fDetector->GetAbsorberMaterial();
207 G4double length = fDetector->GetAbsorberThic << 204 G4double length = fDetector->GetAbsorberThickness();
208 G4double density = material->GetDensity(); << 205 G4double density = material->GetDensity();
209 G4String partName = fParticle->GetParticleNa << 206 G4String partName = fParticle->GetParticleName();
210 207
211 G4EmCalculator emCalculator; 208 G4EmCalculator emCalculator;
212 G4double dEdxTable = 0., dEdxFull = 0.; 209 G4double dEdxTable = 0., dEdxFull = 0.;
213 if (fParticle->GetPDGCharge() != 0.) { << 210 if (fParticle->GetPDGCharge()!= 0.) {
214 dEdxTable = emCalculator.GetDEDX(fEkin, fP << 211 dEdxTable = emCalculator.GetDEDX(fEkin,fParticle,material);
215 dEdxFull = emCalculator.ComputeTotalDEDX(f << 212 dEdxFull = emCalculator.ComputeTotalDEDX(fEkin,fParticle,material);
216 } 213 }
217 G4double stopTable = dEdxTable / density; << 214 G4double stopTable = dEdxTable/density;
218 G4double stopFull = dEdxFull / density; << 215 G4double stopFull = dEdxFull /density;
219 << 216
220 // Stopping Power from simulation. << 217 //Stopping Power from simulation.
221 // << 218 //
222 G4double meandEdx = fEnergyDeposit / length; << 219 G4double meandEdx = fEnergyDeposit/length;
223 G4double stopPower = meandEdx / density; << 220 G4double stopPower = meandEdx/density;
224 221
225 G4cout << "\n ======================== run s 222 G4cout << "\n ======================== run summary ======================\n";
226 223
227 G4int prec = G4cout.precision(3); 224 G4int prec = G4cout.precision(3);
228 << 225
229 G4cout << "\n The run was " << TotNbofEvents 226 G4cout << "\n The run was " << TotNbofEvents << " " << partName << " of "
230 << G4BestUnit(fEkin, "Energy") << " t << 227 << G4BestUnit(fEkin,"Energy") << " through "
231 << material->GetName() << " (density: << 228 << G4BestUnit(length,"Length") << " of "
232 << G4endl; << 229 << material->GetName() << " (density: "
233 << 230 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
>> 231
234 G4cout.precision(4); 232 G4cout.precision(4);
235 << 233
236 G4cout << "\n Total energy deposit in absorb 234 G4cout << "\n Total energy deposit in absorber per event = "
237 << G4BestUnit(fEnergyDeposit, "Energy << 235 << G4BestUnit(fEnergyDeposit,"Energy") << " +- "
>> 236 << G4BestUnit(rmsEdep, "Energy")
238 << G4endl; 237 << G4endl;
239 << 238
240 G4cout << "\n -----> Mean dE/dx = " << meand << 239 G4cout << "\n -----> Mean dE/dx = " << meandEdx/(MeV/cm) << " MeV/cm"
241 << "\t(" << stopPower / (MeV * cm2 / << 240 << "\t(" << stopPower/(MeV*cm2/g) << " MeV*cm2/g)"
242 << 241 << G4endl;
243 G4cout << "\n From formulas :" << G4endl; << 242
244 G4cout << " restricted dEdx = " << dEdxTab << 243 G4cout << "\n From formulas :" << G4endl;
245 << "\t(" << stopTable / (MeV * cm2 / << 244 G4cout << " restricted dEdx = " << dEdxTable/(MeV/cm) << " MeV/cm"
246 << 245 << "\t(" << stopTable/(MeV*cm2/g) << " MeV*cm2/g)"
247 G4cout << " full dEdx = " << dEdxFul << 246 << G4endl;
248 << "\t(" << stopFull / (MeV * cm2 / g << 247
249 << 248 G4cout << " full dEdx = " << dEdxFull/(MeV/cm) << " MeV/cm"
250 G4cout << "\n Leakage : primary = " << G4Be << 249 << "\t(" << stopFull/(MeV*cm2/g) << " MeV*cm2/g)"
251 << G4BestUnit(rmsEl0, "Energy") << 250 << G4endl;
252 << " secondaries = " << G4BestUnit( << 251
253 << G4BestUnit(rmsEl1, "Energy") << G4 << 252 G4cout << "\n Leakage : primary = "
254 << 253 << G4BestUnit(fEnergyLeak[0],"Energy") << " +- "
255 G4cout << " Energy balance : edep + eleak = << 254 << G4BestUnit(rmsEl0, "Energy")
256 << 255 << " secondaries = "
257 G4cout << "\n Total track length (charged) i << 256 << G4BestUnit(fEnergyLeak[1],"Energy") << " +- "
258 << G4BestUnit(fTrakLenCharged, "Lengt << 257 << G4BestUnit(rmsEl1, "Energy")
259 << G4endl; 258 << G4endl;
>> 259
>> 260 G4cout << " Energy balance : edep + eleak = "
>> 261 << G4BestUnit(EnergyBalance,"Energy")
>> 262 << G4endl;
>> 263
>> 264 G4cout << "\n Total track length (charged) in absorber per event = "
>> 265 << G4BestUnit(fTrakLenCharged,"Length") << " +- "
>> 266 << G4BestUnit(rmsTLCh, "Length") << G4endl;
260 267
261 G4cout << " Total track length (neutral) in 268 G4cout << " Total track length (neutral) in absorber per event = "
262 << G4BestUnit(fTrakLenNeutral, "Lengt << 269 << G4BestUnit(fTrakLenNeutral,"Length") << " +- "
263 << G4endl; << 270 << G4BestUnit(rmsTLNe, "Length") << G4endl;
264 271
265 G4cout << "\n Number of steps (charged) in a << 272 G4cout << "\n Number of steps (charged) in absorber per event = "
266 << rmsStCh << G4endl; << 273 << fNbStepsCharged << " +- " << rmsStCh << G4endl;
267 274
268 G4cout << " Number of steps (neutral) in abs << 275 G4cout << " Number of steps (neutral) in absorber per event = "
269 << rmsStNe << G4endl; << 276 << fNbStepsNeutral << " +- " << rmsStNe << G4endl;
270 277
271 G4cout << "\n Number of secondaries per even << 278 G4cout << "\n Number of secondaries per event : Gammas = " << Gamma
272 << "; positrons = " << Posit << G4e << 279 << "; electrons = " << Elect
>> 280 << "; positrons = " << Posit << G4endl;
273 281
274 G4cout << "\n Number of events with the prim << 282 G4cout << "\n Number of events with the primary particle transmitted = "
275 << G4endl; << 283 << transmit[1] << " %" << G4endl;
276 284
277 G4cout << " Number of events with at least 285 G4cout << " Number of events with at least 1 particle transmitted "
278 << "(same charge as primary) = " << t 286 << "(same charge as primary) = " << transmit[0] << " %" << G4endl;
279 287
280 G4cout << "\n Number of events with the prim << 288 G4cout << "\n Number of events with the primary particle reflected = "
281 << G4endl; << 289 << reflect[1] << " %" << G4endl;
282 290
283 G4cout << " Number of events with at least 291 G4cout << " Number of events with at least 1 particle reflected "
284 << "(same charge as primary) = " << r 292 << "(same charge as primary) = " << reflect[0] << " %" << G4endl;
285 293
286 // compute width of the Gaussian central par 294 // compute width of the Gaussian central part of the MultipleScattering
287 // 295 //
288 G4cout << "\n MultipleScattering:" << 296 G4cout << "\n MultipleScattering:"
289 << "\n rms proj angle of transmit pr << 297 << "\n rms proj angle of transmit primary particle = "
290 << " mrad (central part only)" << G4e << 298 << rmsMsc/mrad << " mrad (central part only)" << G4endl;
291 << 299
292 G4cout << " computed theta0 (Highland formu << 300 G4cout << " computed theta0 (Highland formula) = "
293 << " mrad" << G4endl; << 301 << ComputeMscHighland()/mrad << " mrad" << G4endl;
294 << 302
295 G4cout << " central part defined as +- " << << 303 G4cout << " central part defined as +- "
>> 304 << fMscThetaCentral/mrad << " mrad; "
296 << " Tail ratio = " << tailMsc << " 305 << " Tail ratio = " << tailMsc << " %" << G4endl;
297 306
298 // gamma process counts << 307 G4cout << "## Gamma process counts:" << G4endl;
299 // <<
300 G4cout << "\n Gamma process counts:" << G4en <<
301 G4cout << " Photoeffect " << fTypes[0] << 308 G4cout << " Photoeffect " << fTypes[0] << G4endl;
302 G4cout << " Compton " << fTypes[1] << 309 G4cout << " Compton " << fTypes[1] << G4endl;
303 G4cout << " Conversion " << fTypes[2] << 310 G4cout << " Conversion " << fTypes[2] << G4endl;
304 G4cout << " Rayleigh " << fTypes[3] << 311 G4cout << " Rayleigh " << fTypes[3] << G4endl;
305 << 312
306 // normalize histograms 313 // normalize histograms
307 // 314 //
308 G4AnalysisManager* analysisManager = G4Analy 315 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
309 << 316
310 G4int ih = 1; 317 G4int ih = 1;
311 G4double binWidth = analysisManager->GetH1Wi 318 G4double binWidth = analysisManager->GetH1Width(ih);
312 G4double fac = 1. / (TotNbofEvents * binWidt << 319 G4double fac = 1./(TotNbofEvents*binWidth);
313 analysisManager->ScaleH1(ih, fac); << 320 analysisManager->ScaleH1(ih,fac);
314 321
315 ih = 10; 322 ih = 10;
316 binWidth = analysisManager->GetH1Width(ih); 323 binWidth = analysisManager->GetH1Width(ih);
317 fac = 1. / (TotNbofEvents * binWidth); << 324 fac = 1./(TotNbofEvents*binWidth);
318 analysisManager->ScaleH1(ih, fac); << 325 analysisManager->ScaleH1(ih,fac);
319 326
320 ih = 12; 327 ih = 12;
321 analysisManager->ScaleH1(ih, 1. / TotNbofEve << 328 analysisManager->ScaleH1(ih,1./TotNbofEvents);
322 329
323 // reset default precision 330 // reset default precision
324 G4cout.precision(prec); 331 G4cout.precision(prec);
325 } << 332 }
326 333
327 //....oooOO0OOooo........oooOO0OOooo........oo 334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
328 335
329 G4double Run::ComputeMscHighland() 336 G4double Run::ComputeMscHighland()
330 { 337 {
331 // compute the width of the Gaussian central << 338 //compute the width of the Gaussian central part of the MultipleScattering
332 // projected angular distribution. << 339 //projected angular distribution.
333 // Eur. Phys. Jour. C15 (2000) page 166, for << 340 //Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
334 341
335 G4double t = << 342 G4double t = (fDetector->GetAbsorberThickness())
336 (fDetector->GetAbsorberThickness()) / (fDe << 343 /(fDetector->GetAbsorberMaterial()->GetRadlen());
337 if (t < DBL_MIN) return 0.; 344 if (t < DBL_MIN) return 0.;
338 345
339 G4double T = fEkin; 346 G4double T = fEkin;
340 G4double M = fParticle->GetPDGMass(); 347 G4double M = fParticle->GetPDGMass();
341 G4double z = std::abs(fParticle->GetPDGCharg << 348 G4double z = std::abs(fParticle->GetPDGCharge()/eplus);
342 349
343 G4double bpc = T * (T + 2 * M) / (T + M); << 350 G4double bpc = T*(T+2*M)/(T+M);
344 G4double teta0 = 13.6 * MeV * z * std::sqrt( << 351 G4double teta0 = 13.6*MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc;
345 return teta0; 352 return teta0;
346 } 353 }
347 354
348 //....oooOO0OOooo........oooOO0OOooo........oo 355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
349 356