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 medical/fanoCavity2/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 "DetectorConstruction.hh" 36 #include "HistoManager.hh" 37 #include "PrimaryGeneratorAction.hh" 38 39 #include "G4Electron.hh" 40 #include "G4EmCalculator.hh" 41 #include "G4Run.hh" 42 #include "G4RunManager.hh" 43 #include "G4SystemOfUnits.hh" 44 #include "G4UnitsTable.hh" 45 #include "Randomize.hh" 46 47 #include <iomanip> 48 49 //....oooOO0OOooo........oooOO0OOooo........oo 50 51 Run::Run(DetectorConstruction* det, PrimaryGen 52 : G4Run(), 53 fDetector(det), 54 fKinematic(kin), 55 fProcCounter(0), 56 fEdepCavity(0.), 57 fEdepCavity2(0.), 58 fTrkSegmCavity(0.), 59 fNbEventCavity(0), 60 fStepWall(0.), 61 fStepWall2(0.), 62 fStepCavity(0.), 63 fStepCavity2(0.), 64 fNbStepWall(0), 65 fNbStepCavity(0), 66 fEnergyGun(0.), 67 fMassWall(0.), 68 fMassCavity(0.), 69 fIsMaster(isMaster) 70 { 71 // run conditions 72 // 73 G4ParticleDefinition* particleGun = fKinemat 74 G4String partName = particleGun->GetParticle 75 fEnergyGun = fKinematic->GetParticleGun()->G 76 77 // geometry : effective wall volume 78 // 79 G4double cavityThickness = fDetector->GetCav 80 G4Material* mateCavity = fDetector->GetCavit 81 G4double densityCavity = mateCavity->GetDens 82 fMassCavity = cavityThickness * densityCavit 83 84 G4double wallThickness = fDetector->GetWallT 85 G4Material* mateWall = fDetector->GetWallMat 86 G4double densityWall = mateWall->GetDensity( 87 88 G4EmCalculator emCal; 89 G4double RangeWall = emCal.GetCSDARange(fEne 90 G4double factor = 1.2; 91 G4double effWallThick = factor * RangeWall; 92 if ((effWallThick > wallThickness) || (effWa 93 fMassWall = 2 * effWallThick * densityWall; 94 95 G4double massTotal = fMassWall + fMassCavity 96 G4double fMassWallRatio = fMassWall / massTo 97 fKinematic->RunInitialisation(effWallThick, 98 99 G4double massRatio = fMassCavity / fMassWall 100 101 // check radius 102 // 103 G4double worldRadius = fDetector->GetWallRad 104 G4double RangeCavity = emCal.GetCSDARange(fE 105 106 // G4String partName = particleGun->GetPa 107 108 std::ios::fmtflags mode = G4cout.flags(); 109 G4cout.setf(std::ios::fixed, std::ios::float 110 G4int prec = G4cout.precision(3); 111 112 G4cout << "\n ===================== run cond 113 114 G4cout << "\n The run will be " << numberOfE 115 << G4BestUnit(fEnergyGun, "Energy") < 116 << " of " << mateWall->GetName() 117 << " (density: " << G4BestUnit(densit 118 << "); Mass/cm2 = " << G4BestUnit(fMa 119 << "\n csdaRange: " << G4BestUnit(Ran 120 121 G4cout << "\n the cavity is " << G4BestUnit( 122 << mateCavity->GetName() << " (densit 123 << "); Mass/cm2 = " << G4BestUnit(fMa 124 << " --> massRatio = " << std::setpre 125 126 G4cout.precision(3); 127 G4cout << " Wall radius: " << G4BestUnit(wor 128 << "; range in cavity: " << G4BestUni 129 130 G4cout << "\n ============================== 131 132 // stopping power from EmCalculator 133 // 134 G4double dedxWall = emCal.GetDEDX(fEnergyGun 135 dedxWall /= densityWall; 136 G4double dedxCavity = emCal.GetDEDX(fEnergyG 137 dedxCavity /= densityCavity; 138 139 G4cout << std::setprecision(4) 140 << "\n StoppingPower in wall = " << 141 << "\n in cavity = " << 142 << G4endl; 143 144 // process counter 145 // 146 fProcCounter = new ProcessesCount; 147 148 // charged particles and energy flow in cavi 149 // 150 fPartFlowCavity[0] = fPartFlowCavity[1] = 0; 151 fEnerFlowCavity[0] = fEnerFlowCavity[1] = 0. 152 153 // total energy deposit and charged track se 154 // 155 fEdepCavity = fEdepCavity2 = fTrkSegmCavity 156 fNbEventCavity = 0; 157 158 // stepLenth of charged particles 159 // 160 fStepWall = fStepWall2 = fStepCavity = fStep 161 fNbStepWall = fNbStepCavity = 0; 162 163 // reset default formats 164 G4cout.setf(mode, std::ios::floatfield); 165 G4cout.precision(prec); 166 167 // histograms 168 // 169 G4AnalysisManager* analysisManager = G4Analy 170 if (analysisManager->IsActive()) { 171 analysisManager->OpenFile(); 172 } 173 } 174 175 //....oooOO0OOooo........oooOO0OOooo........oo 176 177 Run::~Run() {} 178 179 //....oooOO0OOooo........oooOO0OOooo........oo 180 181 void Run::CountProcesses(G4String procName) 182 { 183 // does the process already encounted ? 184 size_t nbProc = fProcCounter->size(); 185 size_t i = 0; 186 while ((i < nbProc) && ((*fProcCounter)[i]-> 187 i++; 188 if (i == nbProc) fProcCounter->push_back(new 189 190 (*fProcCounter)[i]->Count(); 191 } 192 193 //....oooOO0OOooo........oooOO0OOooo........oo 194 195 void Run::SurveyConvergence(G4int NbofEvents) 196 { 197 if (NbofEvents == 0) return; 198 199 // beam fluence 200 // 201 G4int Nwall = fKinematic->GetWallCount(); 202 G4int Ncavity = fKinematic->GetCavityCount() 203 G4double Iwall = Nwall / fMassWall; 204 G4double Icavity = Ncavity / fMassCavity; 205 G4double Iratio = Icavity / Iwall; 206 G4double Itot = NbofEvents / (fMassWall + fM 207 G4double energyFluence = fEnergyGun * Itot; 208 209 // total dose in cavity 210 // 211 G4double doseCavity = fEdepCavity / fMassCav 212 G4double ratio = doseCavity / energyFluence; 213 G4double err = 100 * (ratio - 1.); 214 215 std::ios::fmtflags mode = G4cout.flags(); 216 G4cout.setf(std::ios::fixed, std::ios::float 217 G4int prec = G4cout.precision(5); 218 219 G4cout << "--->evntNb= " << NbofEvents << " 220 << " Ic/Iw= " << Iratio << " Ne-_cav= 221 << " doseCavity/Ebeam= " << ratio << 222 << G4endl; 223 224 // reset default formats 225 G4cout.setf(mode, std::ios::floatfield); 226 G4cout.precision(prec); 227 } 228 229 //....oooOO0OOooo........oooOO0OOooo........oo 230 231 void Run::EndOfRun() 232 { // Only call by Master thread 233 234 std::ios::fmtflags mode = G4cout.flags(); 235 G4cout.setf(std::ios::fixed, std::ios::float 236 G4int prec = G4cout.precision(3); 237 238 if (numberOfEvent == 0) return; 239 240 // frequency of processes 241 // 242 G4cout << "\n Process calls frequency --->"; 243 for (size_t i = 0; i < fProcCounter->size(); 244 G4String procName = (*fProcCounter)[i]->Ge 245 G4int count = (*fProcCounter)[i]->GetCount 246 G4cout << " " << procName << "= " << coun 247 } 248 G4cout << G4endl; 249 250 // charged particle flow in cavity 251 // 252 G4cout << "\n Charged particle flow in cavit 253 << "\n Enter --> nbParticles = " 254 << "\t Energy = " << G4BestUnit(fEner 255 << "\n Exit --> nbParticles = " 256 << "\t Energy = " << G4BestUnit(fEner 257 258 if (fPartFlowCavity[0] == 0) return; 259 260 G4int Nwall = fKinematic->GetWallCount(); 261 G4int Ncavity = fKinematic->GetCavityCount() 262 263 G4double Iwall = Nwall / fMassWall; 264 G4double Icavity = Ncavity / fMassCavity; 265 G4double Iratio = Icavity / Iwall; 266 G4double Itot = numberOfEvent / (fMassWall + 267 G4double energyFluence = fEnergyGun * Itot; 268 269 G4cout.precision(5); 270 G4cout << "\n beamFluence in wall = " << Nwa 271 << "\t Icav/Iwall = " << Iratio 272 << "\t energyFluence = " << energyFlu 273 274 // error on Edep in cavity 275 // 276 if (fNbEventCavity == 0) return; 277 G4double meanEdep = fEdepCavity / fNbEventCa 278 G4double meanEdep2 = fEdepCavity2 / fNbEvent 279 G4double varianceEdep = meanEdep2 - meanEdep 280 G4double dEoverE = 0.; 281 if (varianceEdep > 0.) dEoverE = std::sqrt(v 282 283 // total dose in cavity 284 // 285 G4double doseCavity = fEdepCavity / fMassCav 286 287 G4double ratio = doseCavity / energyFluence, 288 289 G4cout << "\n Total edep in cavity = " << G4 290 << 100 * dEoverE << " %" 291 << "\n Total dose in cavity = " << do 292 << " +- " << 100 * dEoverE << " %" 293 << "\n\n DoseCavity/EnergyFluence = " 294 295 // track length in cavity 296 G4double meantrack = fTrkSegmCavity / fPartF 297 298 G4cout.precision(4); 299 G4cout << "\n Total charged trackLength in c 300 << " (mean value = " << G4BestUnit( 301 302 // compute mean step size of charged particl 303 // 304 fStepWall /= fNbStepWall; 305 fStepWall2 /= fNbStepWall; 306 G4double rms = fStepWall2 - fStepWall * fSte 307 if (rms > 0.) 308 rms = std::sqrt(rms); 309 else 310 rms = 0.; 311 G4double nbTrackWall = fKinematic->GetWallCo 312 313 G4cout << "\n StepSize of ch. tracks in wall 314 << G4BestUnit(rms, "Length") << "\t ( 315 << ")"; 316 317 fStepCavity /= fNbStepCavity; 318 fStepCavity2 /= fNbStepCavity; 319 rms = fStepCavity2 - fStepCavity * fStepCavi 320 if (rms > 0.) 321 rms = std::sqrt(rms); 322 else 323 rms = 0.; 324 325 G4cout << "\n StepSize of ch. tracks in cavi 326 << G4BestUnit(rms, "Length") 327 << "\t (nbSteps/track = " << double(f 328 329 G4cout << G4endl; 330 331 // reset default formats 332 G4cout.setf(mode, std::ios::floatfield); 333 G4cout.precision(prec); 334 335 // delete and remove all contents in fProcCo 336 while (fProcCounter->size() > 0) { 337 OneProcessCount* aProcCount = fProcCounter 338 fProcCounter->pop_back(); 339 delete aProcCount; 340 } 341 delete fProcCounter; 342 343 // show Rndm status 344 CLHEP::HepRandom::showEngineStatus(); 345 } 346 347 //....oooOO0OOooo........oooOO0OOooo........oo 348 349 void Run::Merge(const G4Run* run) 350 { 351 const Run* localRun = static_cast<const Run* 352 353 // Merge Run variables 354 fPartFlowCavity[0] += localRun->fPartFlowCav 355 fPartFlowCavity[1] += localRun->fPartFlowCav 356 fEnerFlowCavity[0] += localRun->fEnerFlowCav 357 fEnerFlowCavity[1] += localRun->fEnerFlowCav 358 fEdepCavity += localRun->fEdepCavity; 359 fEdepCavity2 += localRun->fEdepCavity2; 360 fTrkSegmCavity += localRun->fTrkSegmCavity; 361 fNbEventCavity += localRun->fNbEventCavity; 362 fStepWall += localRun->fStepWall; 363 fStepWall2 += localRun->fStepWall2; 364 fStepCavity += localRun->fStepCavity; 365 fStepCavity2 += localRun->fStepCavity2; 366 fNbStepWall += localRun->fNbStepWall; 367 fNbStepCavity += localRun->fNbStepCavity; 368 369 // Merge PrimaryGenerator variables 370 fKinematic->AddWallCount(localRun->fKinemati 371 fKinematic->AddCavityCount(localRun->fKinema 372 373 // Merge ProcessCount varaibles 374 std::vector<OneProcessCount*>::iterator it; 375 for (it = localRun->fProcCounter->begin(); i 376 OneProcessCount* process = *it; 377 for (G4int i = 0; i < process->GetCounter( 378 this->CountProcesses(process->GetName()) 379 } 380 381 G4Run::Merge(run); 382 } 383 384 //....oooOO0OOooo........oooOO0OOooo........oo 385