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