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/TestEm5/src/RunActio << 26 // $Id: RunAction.cc,v 1.8 2009/09/18 17:34:54 maire Exp $ 27 /// \brief Implementation of the RunAction cla << 27 // GEANT4 tag $Name: geant4-09-03-patch-02 $ 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 #include "DetectorConstruction.hh" << 33 35 #include "PrimaryGeneratorAction.hh" 34 #include "PrimaryGeneratorAction.hh" 36 #include "HistoManager.hh" 35 #include "HistoManager.hh" 37 #include "Run.hh" << 38 36 39 #include "G4Run.hh" 37 #include "G4Run.hh" 40 #include "G4RunManager.hh" 38 #include "G4RunManager.hh" 41 #include "G4UnitsTable.hh" 39 #include "G4UnitsTable.hh" >> 40 #include "G4Geantino.hh" 42 41 43 #include "Randomize.hh" 42 #include "Randomize.hh" 44 #include "G4SystemOfUnits.hh" << 45 #include <iomanip> << 46 43 47 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 48 45 49 RunAction::RunAction(DetectorConstruction* det << 46 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim, 50 :G4UserRunAction(),fDetector(det), fPrimary(ki << 47 HistoManager* hist) 51 { << 48 :detector(det), primary(prim), histoManager(hist) 52 // Book predefined histograms << 49 { 53 fHistoManager = new HistoManager(); << 50 writeFile = false; 54 } 51 } 55 52 56 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 57 54 58 RunAction::~RunAction() 55 RunAction::~RunAction() 59 { << 56 { } 60 delete fHistoManager; << 57 >> 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 59 >> 60 void RunAction::BeginOfRunAction(const G4Run* aRun) >> 61 { >> 62 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl; >> 63 >> 64 // save Rndm status >> 65 // >> 66 G4RunManager::GetRunManager()->SetRandomNumberStore(true); >> 67 CLHEP::HepRandom::showEngineStatus(); >> 68 >> 69 //initialize cumulative quantities >> 70 // >> 71 G4int nbPixels = detector->GetSizeVectorPixels(); >> 72 G4int size = totalEnergy.size(); >> 73 if (size < nbPixels) { >> 74 visibleEnergy.resize(nbPixels); >> 75 totalEnergy.resize(nbPixels); >> 76 >> 77 visibleEnergy2.resize(nbPixels); >> 78 totalEnergy2.resize(nbPixels); >> 79 } >> 80 >> 81 for (G4int k=0; k<nbPixels; k++) { >> 82 visibleEnergy[k] = visibleEnergy2[k] = totalEnergy[k]= totalEnergy2[k] = 0.0; >> 83 } >> 84 >> 85 G4int n1pxl = detector->GetN1Pixels(); >> 86 size = layerEtot.size(); >> 87 if (size < n1pxl) { >> 88 layerEvis.resize(n1pxl); >> 89 layerEtot.resize(n1pxl); >> 90 layerEvis2.resize(n1pxl); >> 91 layerEtot2.resize(n1pxl); >> 92 } >> 93 >> 94 for (G4int k=0; k<n1pxl; k++) { >> 95 layerEvis[k] = layerEvis2[k] = layerEtot[k]= layerEtot2[k] = 0.0; >> 96 } >> 97 >> 98 nbEvents = 0; >> 99 calorEvis = calorEvis2 = calorEtot = calorEtot2 = Eleak = Eleak2 = 0.; >> 100 EdLeak[0] = EdLeak[1] = EdLeak[2] = 0.; >> 101 nbRadLen = nbRadLen2 = 0.; >> 102 >> 103 //histograms >> 104 // >> 105 histoManager->book(); >> 106 >> 107 //create ascii file for pixels >> 108 // >> 109 if (writeFile) CreateFilePixels(); 61 } 110 } 62 111 63 //....oooOO0OOooo........oooOO0OOooo........oo 112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 64 113 65 G4Run* RunAction::GenerateRun() << 114 void RunAction::fillPerEvent_1(G4int pixel, G4double Evis, G4double Etot) 66 { << 115 { 67 fRun = new Run(fDetector); << 116 //accumulate statistic per pixel 68 return fRun; << 117 // >> 118 visibleEnergy[pixel] += Evis; visibleEnergy2[pixel] += Evis*Evis; >> 119 totalEnergy[pixel] += Etot; totalEnergy2[pixel] += Etot*Etot; 69 } 120 } 70 121 71 //....oooOO0OOooo........oooOO0OOooo........oo 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 72 123 73 void RunAction::BeginOfRunAction(const G4Run*) << 124 void RunAction::fillPerEvent_2(G4int layer, G4double Evis, G4double Etot) 74 { 125 { 75 // save Rndm status << 126 //accumulate statistic per layer 76 //// G4RunManager::GetRunManager()->SetRand << 127 // 77 if (isMaster) G4Random::showEngineStatus(); << 128 layerEvis[layer] += Evis; layerEvis2[layer] += Evis*Evis; 78 << 129 layerEtot[layer] += Etot; layerEtot2[layer] += Etot*Etot; 79 // keep run condition << 130 } 80 if ( fPrimary ) { << 131 81 G4ParticleDefinition* particle << 132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 82 = fPrimary->GetParticleGun()->GetParticl << 133 83 G4double energy = fPrimary->GetParticleGun << 134 void RunAction::fillPerEvent_3(G4double calEvis, G4double calEtot, 84 fRun->SetPrimary(particle, energy); << 135 G4double eleak) 85 } << 136 { 86 << 137 //accumulate statistic 87 //histograms << 138 // 88 // << 139 nbEvents++; 89 G4AnalysisManager* analysisManager = G4Analy << 140 calorEvis += calEvis; calorEvis2 += calEvis*calEvis; 90 if ( analysisManager->IsActive() ) { << 141 calorEtot += calEtot; calorEtot2 += calEtot*calEtot; 91 analysisManager->OpenFile(); << 142 Eleak += eleak; Eleak2 += eleak*eleak; 92 } << 143 } >> 144 >> 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 146 >> 147 void RunAction::fillDetailedLeakage(G4int icase, G4double energy) >> 148 { >> 149 //forward, backward, lateral leakage >> 150 // >> 151 EdLeak[icase] += energy; 93 } 152 } 94 153 95 //....oooOO0OOooo........oooOO0OOooo........oo 154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 96 155 >> 156 void RunAction::fillNbRadLen(G4double dn) >> 157 { >> 158 //total number of radiation length >> 159 // >> 160 nbRadLen += dn; nbRadLen2 += dn*dn; >> 161 } >> 162 >> 163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 164 >> 165 97 void RunAction::EndOfRunAction(const G4Run*) 166 void RunAction::EndOfRunAction(const G4Run*) 98 { << 167 { 99 // print Run summary << 168 //calorimeter 100 // 169 // 101 if (isMaster) fRun->EndOfRun(); << 170 detector->PrintCalorParameters(); 102 171 103 // save histograms << 172 //run conditions 104 G4AnalysisManager* analysisManager = G4Analy << 173 // 105 if ( analysisManager->IsActive() ) { << 174 G4ParticleDefinition* particle = primary->GetParticleGun() 106 analysisManager->Write(); << 175 ->GetParticleDefinition(); 107 analysisManager->CloseFile(); << 176 G4String partName = particle->GetParticleName(); 108 } << 177 G4double energy = primary->GetParticleGun()->GetParticleEnergy(); >> 178 >> 179 G4int prec = G4cout.precision(3); >> 180 >> 181 G4cout << " The run was " << nbEvents << " " << partName << " of " >> 182 << G4BestUnit(energy,"Energy") << " through the calorimeter" << G4endl; >> 183 >> 184 G4cout << "------------------------------------------------------------" >> 185 << G4endl; >> 186 >> 187 //if no events, return >> 188 // >> 189 if (nbEvents == 0) return; >> 190 >> 191 //compute and print statistic >> 192 // >> 193 std::ios::fmtflags mode = G4cout.flags(); >> 194 >> 195 //number of radiation length >> 196 // >> 197 if (particle == G4Geantino::Geantino() ) { >> 198 G4double meanNbRadL = nbRadLen/ nbEvents; >> 199 G4double meanNbRadL2 = nbRadLen2/nbEvents; >> 200 G4double varNbRadL = meanNbRadL2 - meanNbRadL*meanNbRadL; >> 201 G4double rmsNbRadL = 0.; >> 202 if (varNbRadL > 0.) rmsNbRadL = std::sqrt(varNbRadL); >> 203 G4double effRadL = (detector->GetCalorThickness())/meanNbRadL; >> 204 G4cout.precision(5); >> 205 G4cout >> 206 << "\n Calor : mean number of Rad Length = " >> 207 << meanNbRadL << " +- "<< rmsNbRadL >> 208 << " --> Effective Rad Length = " >> 209 << G4BestUnit( effRadL,"Length") << G4endl; >> 210 >> 211 G4cout << "------------------------------------------------------------" >> 212 << G4endl; >> 213 } >> 214 >> 215 G4cout.precision(prec); >> 216 G4cout << "\n " >> 217 << "visible Energy (rms/mean) " >> 218 << "total Energy (rms/mean)" << G4endl; >> 219 >> 220 G4double meanEvis,meanEvis2,varianceEvis,rmsEvis,resEvis; >> 221 G4double meanEtot,meanEtot2,varianceEtot,rmsEtot,resEtot; >> 222 >> 223 G4int n1pxl = detector->GetN1Pixels(); >> 224 >> 225 for (G4int i1=0; i1<n1pxl; i1++) { >> 226 //visible energy >> 227 meanEvis = layerEvis[i1] /nbEvents; >> 228 meanEvis2 = layerEvis2[i1]/nbEvents; >> 229 varianceEvis = meanEvis2 - meanEvis*meanEvis; >> 230 rmsEvis = 0.; >> 231 if (varianceEvis > 0.) rmsEvis = std::sqrt(varianceEvis); >> 232 resEvis = 100*rmsEvis/meanEvis; >> 233 histoManager->FillHisto(3, i1+0.5, meanEvis); >> 234 >> 235 //total energy >> 236 meanEtot = layerEtot[i1] /nbEvents; >> 237 meanEtot2 = layerEtot2[i1]/nbEvents; >> 238 varianceEtot = meanEtot2 - meanEtot*meanEtot; >> 239 rmsEtot = 0.; >> 240 if (varianceEtot > 0.) rmsEtot = std::sqrt(varianceEtot); >> 241 resEtot = 100*rmsEtot/meanEtot; >> 242 histoManager->FillHisto(4, i1+0.5, meanEtot); >> 243 >> 244 //print >> 245 // >> 246 G4cout >> 247 << "\n layer " << i1 << ": " >> 248 << std::setprecision(5) >> 249 << std::setw(6) << G4BestUnit(meanEvis,"Energy") << " +- " >> 250 << std::setprecision(4) >> 251 << std::setw(5) << G4BestUnit( rmsEvis,"Energy") << " (" >> 252 << std::setprecision(2) >> 253 << std::setw(3) << resEvis << " %)" >> 254 << " " >> 255 << std::setprecision(5) >> 256 << std::setw(6) << G4BestUnit(meanEtot,"Energy") << " +- " >> 257 << std::setprecision(4) >> 258 << std::setw(5) << G4BestUnit( rmsEtot,"Energy") << " (" >> 259 << std::setprecision(2) >> 260 << std::setw(3) << resEtot << " %)"; >> 261 } >> 262 G4cout << G4endl; >> 263 >> 264 //calorimeter: visible energy >> 265 meanEvis = calorEvis /nbEvents; >> 266 meanEvis2 = calorEvis2/nbEvents; >> 267 varianceEvis = meanEvis2 - meanEvis*meanEvis; >> 268 rmsEvis = 0.; >> 269 if (varianceEvis > 0.) rmsEvis = std::sqrt(varianceEvis); >> 270 resEvis = 100*rmsEvis/meanEvis; >> 271 >> 272 //calorimeter: total energy >> 273 meanEtot = calorEtot /nbEvents; >> 274 meanEtot2 = calorEtot2/nbEvents; >> 275 varianceEtot = meanEtot2 - meanEtot*meanEtot; >> 276 rmsEtot = 0.; >> 277 if (varianceEtot > 0.) rmsEtot = std::sqrt(varianceEtot); >> 278 resEtot = 100*rmsEtot/meanEtot; >> 279 >> 280 //print >> 281 // >> 282 G4cout >> 283 << "\n total calor : " >> 284 << std::setprecision(5) >> 285 << std::setw(6) << G4BestUnit(meanEvis,"Energy") << " +- " >> 286 << std::setprecision(4) >> 287 << std::setw(5) << G4BestUnit( rmsEvis,"Energy") << " (" >> 288 << std::setprecision(2) >> 289 << std::setw(3) << resEvis << " %)" >> 290 << " " >> 291 << std::setprecision(5) >> 292 << std::setw(6) << G4BestUnit(meanEtot,"Energy") << " +- " >> 293 << std::setprecision(4) >> 294 << std::setw(5) << G4BestUnit( rmsEtot,"Energy") << " (" >> 295 << std::setprecision(2) >> 296 << std::setw(3) << resEtot << " %)"; >> 297 >> 298 G4cout << "\n------------------------------------------------------------" >> 299 << G4endl; >> 300 >> 301 //leakage >> 302 G4double meanEleak,meanEleak2,varianceEleak,rmsEleak,ratio; >> 303 meanEleak = Eleak /nbEvents; >> 304 meanEleak2 = Eleak2/nbEvents; >> 305 varianceEleak = meanEleak2 - meanEleak*meanEleak; >> 306 rmsEleak = 0.; >> 307 if (varianceEleak > 0.) rmsEleak = std::sqrt(varianceEleak); >> 308 ratio = 100*meanEleak/energy; >> 309 >> 310 G4double forward = 100*EdLeak[0]/(nbEvents*energy); >> 311 G4double bakward = 100*EdLeak[1]/(nbEvents*energy); >> 312 G4double lateral = 100*EdLeak[2]/(nbEvents*energy); >> 313 //print >> 314 // >> 315 G4cout >> 316 << "\n Leakage : " >> 317 << std::setprecision(5) >> 318 << std::setw(6) << G4BestUnit(meanEleak,"Energy") << " +- " >> 319 << std::setprecision(4) >> 320 << std::setw(5) << G4BestUnit( rmsEleak,"Energy") >> 321 << "\n Eleak/Ebeam =" >> 322 << std::setprecision(3) >> 323 << std::setw(4) << ratio << " % ( forward =" >> 324 << std::setw(4) << forward << " %; backward =" >> 325 << std::setw(4) << bakward << " %; lateral =" >> 326 << std::setw(4) << lateral << " %)" >> 327 << G4endl; >> 328 >> 329 G4cout.setf(mode,std::ios::floatfield); >> 330 G4cout.precision(prec); >> 331 >> 332 //save histograms >> 333 histoManager->save(); 109 334 110 // show Rndm status 335 // show Rndm status 111 if (isMaster) G4Random::showEngineStatus(); << 336 CLHEP::HepRandom::showEngineStatus(); 112 } 337 } 113 338 114 //....oooOO0OOooo........oooOO0OOooo........oo 339 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 340 >> 341 #include <fstream> >> 342 >> 343 void RunAction::CreateFilePixels() >> 344 { >> 345 //create file and write run header >> 346 // >> 347 G4String name = histoManager->GetFileName(); >> 348 G4String fileName = name + ".pixels.ascii"; >> 349 >> 350 std::ofstream File(fileName, std::ios::out); >> 351 >> 352 G4int n1pxl = detector->GetN1Pixels(); >> 353 G4int n2pxl = detector->GetN2Pixels(); >> 354 G4int n1shift = detector->GetN1Shift(); >> 355 G4int nbEvents = G4RunManager::GetRunManager()->GetCurrentRun() >> 356 ->GetNumberOfEventToBeProcessed(); >> 357 File << nbEvents << " " << n1pxl << " " << n2pxl << " " << n1shift >> 358 << G4endl; >> 359 } >> 360 >> 361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 362 >> 363 115 364