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