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