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