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