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/TestEm3/src/RunActio 26 /// \file electromagnetic/TestEm3/src/RunAction.cc 27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction class 28 // 28 // >> 29 // $Id$ 29 // 30 // 30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 33 33 #include "RunAction.hh" 34 #include "RunAction.hh" 34 35 35 #include "DetectorConstruction.hh" << 36 #include "HistoManager.hh" << 37 #include "PrimaryGeneratorAction.hh" 36 #include "PrimaryGeneratorAction.hh" 38 #include "Run.hh" << 39 #include "RunActionMessenger.hh" 37 #include "RunActionMessenger.hh" >> 38 #include "HistoManager.hh" >> 39 #include "EmAcceptance.hh" 40 40 >> 41 #include "G4Run.hh" 41 #include "G4RunManager.hh" 42 #include "G4RunManager.hh" 42 #include "G4Timer.hh" << 43 >> 44 #include "G4ParticleTable.hh" >> 45 #include "G4ParticleDefinition.hh" >> 46 #include "G4Track.hh" >> 47 #include "G4Gamma.hh" >> 48 #include "G4Electron.hh" >> 49 #include "G4Positron.hh" >> 50 #include "G4ProductionCutsTable.hh" >> 51 #include "G4LossTableManager.hh" >> 52 >> 53 #include "G4UnitsTable.hh" >> 54 #include "G4SystemOfUnits.hh" >> 55 43 #include "Randomize.hh" 56 #include "Randomize.hh" 44 57 45 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 46 59 47 RunAction::RunAction(DetectorConstruction* det 60 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim) 48 : fDetector(det), fPrimary(prim) << 61 :fDetector(det), fPrimary(prim), fRunMessenger(0), fHistoManager(0) 49 { 62 { 50 fRunMessenger = new RunActionMessenger(this) 63 fRunMessenger = new RunActionMessenger(this); 51 fHistoManager = new HistoManager(); 64 fHistoManager = new HistoManager(); >> 65 fApplyLimit = false; >> 66 >> 67 fChargedStep = fNeutralStep = 0.0; >> 68 >> 69 for (G4int k=0; k<MaxAbsor; k++) { fEdeptrue[k] = fRmstrue[k] = 1.; >> 70 fLimittrue[k] = DBL_MAX; >> 71 } 52 } 72 } 53 73 54 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 55 75 56 RunAction::~RunAction() 76 RunAction::~RunAction() 57 { 77 { 58 delete fRunMessenger; 78 delete fRunMessenger; 59 } 79 } 60 80 61 //....oooOO0OOooo........oooOO0OOooo........oo 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 62 82 63 G4Run* RunAction::GenerateRun() << 83 void RunAction::BeginOfRunAction(const G4Run* aRun) 64 { 84 { 65 fRun = new Run(fDetector); << 85 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl; 66 return fRun; << 67 } << 68 86 69 //....oooOO0OOooo........oooOO0OOooo........oo << 87 // save Rndm status >> 88 // >> 89 ////G4RunManager::GetRunManager()->SetRandomNumberStore(true); >> 90 CLHEP::HepRandom::showEngineStatus(); 70 91 71 void RunAction::BeginOfRunAction(const G4Run*) << 92 //initialize cumulative quantities 72 { << 93 // 73 // keep run condition << 94 for (G4int k=0; k<MaxAbsor; k++) { 74 if (fPrimary) { << 95 fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.; 75 G4ParticleDefinition* particle = fPrimary- << 96 fEnergyDeposit[k].clear(); 76 G4double energy = fPrimary->GetParticleGun << 77 fRun->SetPrimary(particle, energy); << 78 } 97 } 79 98 80 // histograms << 99 fChargedStep = fNeutralStep = 0.0; >> 100 >> 101 fN_gamma = 0; >> 102 fN_elec = 0; >> 103 fN_pos = 0; >> 104 >> 105 //initialize Eflow >> 106 // >> 107 G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2; >> 108 fEnergyFlow.resize(nbPlanes); >> 109 fLateralEleak.resize(nbPlanes); >> 110 for (G4int k=0; k<nbPlanes; k++) {fEnergyFlow[k] = fLateralEleak[k] = 0.; } >> 111 >> 112 //histograms 81 // 113 // 82 G4AnalysisManager* analysis = G4AnalysisMana 114 G4AnalysisManager* analysis = G4AnalysisManager::Instance(); 83 if (analysis->IsActive()) analysis->OpenFile 115 if (analysis->IsActive()) analysis->OpenFile(); >> 116 >> 117 //example of print dEdx tables >> 118 // >> 119 ////PrintDedxTables(); >> 120 } 84 121 85 // save Rndm status and open the timer << 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 86 123 87 if (isMaster) { << 124 void RunAction::FillPerEvent(G4int kAbs, G4double EAbs, G4double LAbs) 88 // G4Random::showEngineStatus(); << 125 { 89 fTimer = new G4Timer(); << 126 //accumulate statistic with restriction 90 fTimer->Start(); << 127 // 91 } << 128 if(fApplyLimit) fEnergyDeposit[kAbs].push_back(EAbs); >> 129 fSumEAbs[kAbs] += EAbs; fSum2EAbs[kAbs] += EAbs*EAbs; >> 130 fSumLAbs[kAbs] += LAbs; fSum2LAbs[kAbs] += LAbs*LAbs; 92 } 131 } 93 132 94 //....oooOO0OOooo........oooOO0OOooo........oo 133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 95 134 96 void RunAction::EndOfRunAction(const G4Run*) << 135 >> 136 void RunAction::EndOfRunAction(const G4Run* aRun) 97 { 137 { 98 // compute and print statistic << 138 G4int nEvt = aRun->GetNumberOfEvent(); 99 if (isMaster) { << 139 G4double norm = G4double(nEvt); 100 fTimer->Stop(); << 140 if(norm > 0) norm = 1./norm; 101 if (!((G4RunManager::GetRunManager()->GetR << 141 G4double qnorm = std::sqrt(norm); 102 G4cout << "\n" << 142 103 << "Total number of events: " << << 143 fChargedStep *= norm; 104 G4cout << "Master thread time: " << *fT << 144 fNeutralStep *= norm; >> 145 >> 146 //compute and print statistic >> 147 // >> 148 G4double beamEnergy = fPrimary->GetParticleGun()->GetParticleEnergy(); >> 149 G4double sqbeam = std::sqrt(beamEnergy/GeV); >> 150 >> 151 G4double MeanEAbs,MeanEAbs2,rmsEAbs,resolution,rmsres; >> 152 G4double MeanLAbs,MeanLAbs2,rmsLAbs; >> 153 >> 154 std::ios::fmtflags mode = G4cout.flags(); >> 155 G4int prec = G4cout.precision(2); >> 156 G4cout << "\n------------------------------------------------------------\n"; >> 157 G4cout << std::setw(14) << "material" >> 158 << std::setw(17) << "Edep RMS" >> 159 << std::setw(33) << "sqrt(E0(GeV))*rmsE/Emean" >> 160 << std::setw(23) << "total tracklen \n \n"; >> 161 >> 162 for (G4int k=1; k<=fDetector->GetNbOfAbsor(); k++) >> 163 { >> 164 MeanEAbs = fSumEAbs[k]*norm; >> 165 MeanEAbs2 = fSum2EAbs[k]*norm; >> 166 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs)); >> 167 //G4cout << "k= " << k << " RMS= " << rmsEAbs >> 168 // << " fApplyLimit: " << fApplyLimit << G4endl; >> 169 if(fApplyLimit) { >> 170 G4int nn = 0; >> 171 G4double sume = 0.0; >> 172 G4double sume2 = 0.0; >> 173 // compute trancated means >> 174 G4double lim = rmsEAbs * 2.5; >> 175 for(G4int i=0; i<nEvt; i++) { >> 176 G4double e = (fEnergyDeposit[k])[i]; >> 177 if(std::abs(e - MeanEAbs) < lim) { >> 178 sume += e; >> 179 sume2 += e*e; >> 180 nn++; >> 181 } >> 182 } >> 183 G4double norm1 = G4double(nn); >> 184 if(norm1 > 0.0) norm1 = 1.0/norm1; >> 185 MeanEAbs = sume*norm1; >> 186 MeanEAbs2 = sume2*norm1; >> 187 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs)); >> 188 } >> 189 >> 190 resolution= 100.*sqbeam*rmsEAbs/MeanEAbs; >> 191 rmsres = resolution*qnorm; >> 192 >> 193 // Save mean and RMS >> 194 fSumEAbs[k] = MeanEAbs; >> 195 fSum2EAbs[k] = rmsEAbs; >> 196 >> 197 MeanLAbs = fSumLAbs[k]*norm; >> 198 MeanLAbs2 = fSum2LAbs[k]*norm; >> 199 rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs*MeanLAbs)); >> 200 >> 201 //print >> 202 // >> 203 G4cout >> 204 << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName() << ": " >> 205 << std::setprecision(5) >> 206 << std::setw(6) << G4BestUnit(MeanEAbs,"Energy") << " : " >> 207 << std::setprecision(4) >> 208 << std::setw(5) << G4BestUnit( rmsEAbs,"Energy") >> 209 << std::setw(10) << resolution << " +- " >> 210 << std::setw(5) << rmsres << " %" >> 211 << std::setprecision(3) >> 212 << std::setw(10) << G4BestUnit(MeanLAbs,"Length") << " +- " >> 213 << std::setw(4) << G4BestUnit( rmsLAbs,"Length") >> 214 << G4endl; 105 } 215 } 106 delete fTimer; << 216 G4cout << "\n------------------------------------------------------------\n"; 107 fRun->EndOfRun(); << 217 108 } << 218 G4cout << " Beam particle " 109 // save histograms << 219 << fPrimary->GetParticleGun()-> >> 220 GetParticleDefinition()->GetParticleName() >> 221 << " E = " << G4BestUnit(beamEnergy,"Energy") << G4endl; >> 222 G4cout << " Mean number of gamma " << (G4double)fN_gamma*norm << G4endl; >> 223 G4cout << " Mean number of e- " << (G4double)fN_elec*norm << G4endl; >> 224 G4cout << " Mean number of e+ " << (G4double)fN_pos*norm << G4endl; >> 225 G4cout << std::setprecision(6) >> 226 << " Mean number of charged steps " << fChargedStep << G4endl; >> 227 G4cout << " Mean number of neutral steps " << fNeutralStep << G4endl; >> 228 G4cout << "------------------------------------------------------------\n"; >> 229 >> 230 //Energy flow >> 231 // 110 G4AnalysisManager* analysis = G4AnalysisMana 232 G4AnalysisManager* analysis = G4AnalysisManager::Instance(); 111 if (analysis->IsActive()) { << 233 G4int Idmax = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()); >> 234 for (G4int Id=1; Id<=Idmax+1; Id++) { >> 235 analysis->FillH1(2*MaxAbsor+1, (G4double)Id, fEnergyFlow[Id]); >> 236 analysis->FillH1(2*MaxAbsor+2, (G4double)Id, fLateralEleak[Id]); >> 237 } >> 238 >> 239 //Energy deposit from energy flow balance >> 240 // >> 241 G4double EdepTot[MaxAbsor]; >> 242 for (G4int k=0; k<MaxAbsor; k++) EdepTot[k] = 0.; >> 243 >> 244 G4int nbOfAbsor = fDetector->GetNbOfAbsor(); >> 245 for (G4int Id=1; Id<=Idmax; Id++) { >> 246 G4int iAbsor = Id%nbOfAbsor; if (iAbsor==0) iAbsor = nbOfAbsor; >> 247 EdepTot [iAbsor] += (fEnergyFlow[Id] - fEnergyFlow[Id+1] - fLateralEleak[Id]); >> 248 } >> 249 >> 250 G4cout << std::setprecision(3) >> 251 << "\n Energy deposition from Energy flow balance : \n" >> 252 << std::setw(10) << " material \t Total Edep \n \n"; >> 253 G4cout.precision(6); >> 254 >> 255 for (G4int k=1; k<=nbOfAbsor; k++) { >> 256 EdepTot [k] *= norm; >> 257 G4cout << std::setw(10) << fDetector->GetAbsorMaterial(k)->GetName() << ":" >> 258 << "\t " << G4BestUnit(EdepTot [k],"Energy") << "\n"; >> 259 } >> 260 >> 261 G4cout << "\n------------------------------------------------------------\n" >> 262 << G4endl; >> 263 >> 264 G4cout.setf(mode,std::ios::floatfield); >> 265 G4cout.precision(prec); >> 266 >> 267 // Acceptance >> 268 EmAcceptance acc; >> 269 G4bool isStarted = false; >> 270 for (G4int j=1; j<=fDetector->GetNbOfAbsor(); j++) { >> 271 if (fLimittrue[j] < DBL_MAX) { >> 272 if (!isStarted) { >> 273 acc.BeginOfAcceptance("Sampling Calorimeter",nEvt); >> 274 isStarted = true; >> 275 } >> 276 MeanEAbs = fSumEAbs[j]; >> 277 rmsEAbs = fSum2EAbs[j]; >> 278 G4String mat = fDetector->GetAbsorMaterial(j)->GetName(); >> 279 acc.EmAcceptanceGauss("Edep"+mat, nEvt, MeanEAbs, >> 280 fEdeptrue[j], fRmstrue[j], fLimittrue[j]); >> 281 acc.EmAcceptanceGauss("Erms"+mat, nEvt, rmsEAbs, >> 282 fRmstrue[j], fRmstrue[j], 2.0*fLimittrue[j]); >> 283 } >> 284 } >> 285 if(isStarted) acc.EndOfAcceptance(); >> 286 >> 287 //normalize histograms >> 288 // >> 289 for (G4int ih = MaxAbsor+1; ih < MaxHisto; ih++) { >> 290 analysis->ScaleH1(ih,norm/MeV); >> 291 } >> 292 >> 293 //save histograms >> 294 if (analysis->IsActive()) { 112 analysis->Write(); 295 analysis->Write(); 113 analysis->CloseFile(); 296 analysis->CloseFile(); 114 } << 297 } 115 298 116 // show Rndm status 299 // show Rndm status 117 // if (isMaster) G4Random::showEngineStatu << 300 CLHEP::HepRandom::showEngineStatus(); 118 } 301 } 119 302 120 //....oooOO0OOooo........oooOO0OOooo........oo 303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 121 304 122 void RunAction::SetEdepAndRMS(G4int i, G4doubl << 305 void RunAction::PrintDedxTables() >> 306 { >> 307 //Print dE/dx tables with binning identical to the Geant3 JMATE bank. >> 308 //The printout is readable as Geant3 ffread data cards (by the program g4mat). >> 309 // >> 310 const G4double tkmin=10*keV, tkmax=10*TeV; >> 311 const G4int nbin=90; >> 312 G4double tk[nbin]; >> 313 >> 314 const G4int ncolumn = 5; >> 315 >> 316 //compute the kinetic energies >> 317 // >> 318 const G4double dp = std::log10(tkmax/tkmin)/nbin; >> 319 const G4double dt = std::pow(10.,dp); >> 320 tk[0] = tkmin; >> 321 for (G4int i=1; i<nbin; ++i) tk[i] = tk[i-1]*dt; >> 322 >> 323 //print the kinetic energies >> 324 // >> 325 std::ios::fmtflags mode = G4cout.flags(); >> 326 G4cout.setf(std::ios::fixed,std::ios::floatfield); >> 327 G4int prec = G4cout.precision(3); >> 328 >> 329 G4cout << "\n kinetic energies \n "; >> 330 for (G4int j=0; j<nbin; ++j) { >> 331 G4cout << G4BestUnit(tk[j],"Energy") << "\t"; >> 332 if ((j+1)%ncolumn == 0) G4cout << "\n "; >> 333 } >> 334 G4cout << G4endl; >> 335 >> 336 //print the dE/dx tables >> 337 // >> 338 G4cout.setf(std::ios::scientific,std::ios::floatfield); >> 339 >> 340 G4ParticleDefinition* >> 341 part = fPrimary->GetParticleGun()->GetParticleDefinition(); >> 342 >> 343 G4ProductionCutsTable* theCoupleTable = >> 344 G4ProductionCutsTable::GetProductionCutsTable(); >> 345 size_t numOfCouples = theCoupleTable->GetTableSize(); >> 346 const G4MaterialCutsCouple* couple = 0; >> 347 >> 348 for (G4int iab=1;iab <= fDetector->GetNbOfAbsor(); iab++) >> 349 { >> 350 G4Material* mat = fDetector->GetAbsorMaterial(iab); >> 351 G4int index = 0; >> 352 for (size_t i=0; i<numOfCouples; i++) { >> 353 couple = theCoupleTable->GetMaterialCutsCouple(i); >> 354 if (couple->GetMaterial() == mat) {index = i; break;} >> 355 } >> 356 G4cout << "\nLIST"; >> 357 G4cout << "\nC \nC dE/dx (MeV/cm) for " << part->GetParticleName() >> 358 << " in " << mat ->GetName() << "\nC"; >> 359 G4cout << "\nKINE (" << part->GetParticleName() << ")"; >> 360 G4cout << "\nMATE (" << mat ->GetName() << ")"; >> 361 G4cout.precision(2); >> 362 G4cout << "\nERAN " << tkmin/GeV << " (ekmin)\t" >> 363 << tkmax/GeV << " (ekmax)\t" >> 364 << nbin << " (nekbin)"; >> 365 G4double cutgam = >> 366 (*(theCoupleTable->GetEnergyCutsVector(idxG4GammaCut)))[index]; >> 367 if (cutgam < tkmin) cutgam = tkmin; >> 368 if (cutgam > tkmax) cutgam = tkmax; >> 369 G4double cutele = >> 370 (*(theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut)))[index]; >> 371 if (cutele < tkmin) cutele = tkmin; >> 372 if (cutele > tkmax) cutele = tkmax; >> 373 G4cout << "\nCUTS " << cutgam/GeV << " (cutgam)\t" >> 374 << cutele/GeV << " (cutele)"; >> 375 >> 376 G4cout.precision(6); >> 377 G4cout << "\nG4VAL \n "; >> 378 for (G4int l=0;l<nbin; ++l) >> 379 { >> 380 G4double dedx = G4LossTableManager::Instance() >> 381 ->GetDEDX(part,tk[l],couple); >> 382 G4cout << dedx/(MeV/cm) << "\t"; >> 383 if ((l+1)%ncolumn == 0) G4cout << "\n "; >> 384 } >> 385 G4cout << G4endl; >> 386 } >> 387 >> 388 G4cout.precision(prec); >> 389 G4cout.setf(mode,std::ios::floatfield); >> 390 } >> 391 >> 392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 393 >> 394 void RunAction::AddSecondaryTrack(const G4Track* track) 123 { 395 { 124 if (fRun) fRun->SetEdepAndRMS(i, edep, rms, << 396 const G4ParticleDefinition* d = track->GetDefinition(); >> 397 if(d == G4Gamma::Gamma()) { ++fN_gamma; } >> 398 else if (d == G4Electron::Electron()) { ++fN_elec; } >> 399 else if (d == G4Positron::Positron()) { ++fN_pos; } 125 } 400 } 126 401 127 //....oooOO0OOooo........oooOO0OOooo........oo 402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 128 403 129 void RunAction::SetApplyLimit(G4bool val) << 404 void RunAction::SetEdepAndRMS(G4int i, G4double edep, G4double rms, G4double lim) 130 { 405 { 131 if (fRun) fRun->SetApplyLimit(val); << 406 if (i>=0 && i<MaxAbsor) { >> 407 fEdeptrue [i] = edep; >> 408 fRmstrue [i] = rms; >> 409 fLimittrue[i] = lim; >> 410 } 132 } 411 } 133 412 134 //....oooOO0OOooo........oooOO0OOooo........oo 413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 135 414