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