Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 /// \file electromagnetic/TestEm11/src/RunActi << 23 // $Id: RunAction.cc,v 1.1 2005/06/03 15:20:32 maire Exp $ 27 /// \brief Implementation of the RunAction cla << 24 // GEANT4 tag $Name: geant4-07-01-patch-01 $ 28 // << 25 // 29 // << 30 //....oooOO0OOooo........oooOO0OOooo........oo 26 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oo 27 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 28 33 #include "RunAction.hh" 29 #include "RunAction.hh" 34 << 30 #include "DetectorConstruction.hh" 35 #include "HistoManager.hh" << 36 #include "PhysicsList.hh" 31 #include "PhysicsList.hh" 37 #include "PrimaryGeneratorAction.hh" << 38 #include "Run.hh" << 39 #include "StepMax.hh" 32 #include "StepMax.hh" >> 33 #include "PrimaryGeneratorAction.hh" >> 34 #include "HistoManager.hh" 40 35 41 #include "G4EmCalculator.hh" << 42 #include "G4EmParameters.hh" << 43 #include "G4Run.hh" 36 #include "G4Run.hh" 44 #include "G4SystemOfUnits.hh" << 37 #include "G4RunManager.hh" 45 #include "G4UnitsTable.hh" 38 #include "G4UnitsTable.hh" >> 39 #include "G4EmCalculator.hh" >> 40 46 #include "Randomize.hh" 41 #include "Randomize.hh" >> 42 #include <iomanip> 47 43 48 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 49 45 50 RunAction::RunAction(DetectorConstruction* det << 46 RunAction::RunAction(DetectorConstruction* det, PhysicsList* phys, 51 : fDetector(det), fPhysics(phys), fPrimary(k << 47 PrimaryGeneratorAction* kin, HistoManager* histo) 52 { << 48 :detector(det),physics(phys),kinematic(kin),histoManager(histo) 53 // Book predefined histograms << 49 { } 54 fHistoManager = new HistoManager(); << 55 } << 56 50 57 //....oooOO0OOooo........oooOO0OOooo........oo 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 58 52 59 RunAction::~RunAction() 53 RunAction::~RunAction() 60 { << 54 { } 61 delete fHistoManager; << 62 } << 63 55 64 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 65 57 66 G4Run* RunAction::GenerateRun() << 58 void RunAction::BeginOfRunAction(const G4Run* aRun) 67 { << 59 { 68 fRun = new Run(fDetector); << 60 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl; 69 return fRun; << 61 >> 62 // save Rndm status >> 63 G4RunManager::GetRunManager()->SetRandomNumberStore(true); >> 64 HepRandom::showEngineStatus(); >> 65 >> 66 //initialize total energy deposit >> 67 // >> 68 Edeposit = Edeposit2 = 0.; >> 69 >> 70 //initialize track legth of primary >> 71 // >> 72 trackLen = trackLen2 = 0.; >> 73 >> 74 //initialize projected range >> 75 // >> 76 projRange = projRange2 = 0.; >> 77 >> 78 //initialize mean step size >> 79 // >> 80 nbOfSteps = nbOfSteps2 = 0; stepSize = stepSize2 = 0.; >> 81 >> 82 //initialize track status >> 83 // >> 84 status[0] = status[1] = status[2] = 0; >> 85 >> 86 //histograms >> 87 // >> 88 histoManager->book(); 70 } 89 } 71 90 72 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 73 92 74 void RunAction::BeginOfRunAction(const G4Run*) << 93 void RunAction::EndOfRunAction(const G4Run* aRun) 75 { 94 { 76 // show Rndm status << 95 std::ios::fmtflags mode = G4cout.flags(); 77 if (isMaster) { << 96 G4cout.setf(std::ios::fixed,std::ios::floatfield); 78 G4Random::showEngineStatus(); << 97 79 G4EmParameters::Instance()->Dump(); << 98 G4int NbofEvents = aRun->GetNumberOfEvent(); 80 } << 99 if (NbofEvents == 0) return; 81 << 100 G4double fNbofEvents = double(NbofEvents); 82 // keep run condition << 101 83 if (fPrimary) { << 102 //run conditions 84 G4ParticleDefinition* particle = fPrimary- << 103 // 85 G4double energy = fPrimary->GetParticleGun << 104 G4Material* material = detector->GetAbsorMaterial(); 86 fRun->SetPrimary(particle, energy); << 105 G4double density = material->GetDensity(); 87 } << 106 88 << 107 G4ParticleDefinition* particle = kinematic->GetParticleGun() 89 // histograms << 108 ->GetParticleDefinition(); 90 // << 109 G4String partName = particle->GetParticleName(); 91 G4AnalysisManager* analysisManager = G4Analy << 110 G4double energy = kinematic->GetParticleGun()->GetParticleEnergy(); 92 if (analysisManager->IsActive()) { << 111 93 analysisManager->OpenFile(); << 112 G4cout << "\n ======================== run summary ======================\n"; 94 } << 113 95 << 114 G4int prec = G4cout.precision(2); 96 if (!fPrimary) return; << 115 97 << 116 G4cout << "\n The run consists of " << NbofEvents << " "<< partName << " of " 98 // set StepMax from histos 1 and 8 << 117 << G4BestUnit(energy,"Energy") << " through " >> 118 << G4BestUnit(detector->GetAbsorSizeX(),"Length") << " of " >> 119 << material->GetName() << " (density: " >> 120 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl; >> 121 >> 122 G4cout << "\n ============================================================\n"; >> 123 >> 124 //compute total energy deposit >> 125 // >> 126 Edeposit /= NbofEvents; Edeposit2 /= NbofEvents; >> 127 G4double rms = Edeposit2 - Edeposit*Edeposit; >> 128 if (rms>0.) rms = std::sqrt(rms); else rms = 0.; >> 129 >> 130 G4cout.precision(3); >> 131 G4cout >> 132 << "\n Total Energy deposited = " << G4BestUnit(Edeposit,"Energy") >> 133 << " +- " << G4BestUnit( rms,"Energy") >> 134 << G4endl; >> 135 >> 136 //compute track length of primary track >> 137 // >> 138 trackLen /= NbofEvents; trackLen2 /= NbofEvents; >> 139 rms = trackLen2 - trackLen*trackLen; >> 140 if (rms>0.) rms = std::sqrt(rms); else rms = 0.; >> 141 >> 142 G4cout.precision(3); >> 143 G4cout >> 144 << "\n Track length of primary track = " << G4BestUnit(trackLen,"Length") >> 145 << " +- " << G4BestUnit( rms,"Length"); >> 146 >> 147 //compare with csda range >> 148 // >> 149 G4EmCalculator emCalculator; >> 150 G4double csdaRange = 0.; >> 151 if (particle->GetPDGCharge() != 0.) >> 152 csdaRange = emCalculator.GetRange(energy,particle,material); >> 153 G4cout >> 154 << "\n csda Range from EmCalculator = " << G4BestUnit(csdaRange,"Length") >> 155 << G4endl; >> 156 >> 157 //compute projected range of primary track >> 158 // >> 159 projRange /= NbofEvents; projRange2 /= NbofEvents; >> 160 rms = projRange2 - projRange*projRange; >> 161 if (rms>0.) rms = std::sqrt(rms); else rms = 0.; >> 162 >> 163 G4cout >> 164 << "\n Projected range = " << G4BestUnit(projRange,"Length") >> 165 << " +- " << G4BestUnit( rms,"Length") >> 166 << G4endl; >> 167 >> 168 //nb of steps and step size of primary track >> 169 // >> 170 G4double fNbSteps = nbOfSteps/fNbofEvents, fNbSteps2 = nbOfSteps2/fNbofEvents; >> 171 rms = fNbSteps2 - fNbSteps*fNbSteps; >> 172 if (rms>0.) rms = std::sqrt(rms); else rms = 0.; >> 173 >> 174 G4cout.precision(2); >> 175 G4cout << "\n Nb of steps of primary track = " << fNbSteps << " +- " << rms; >> 176 >> 177 stepSize /= NbofEvents; stepSize2 /= NbofEvents; >> 178 rms = stepSize2 - stepSize*stepSize; >> 179 if (rms>0.) rms = std::sqrt(rms); else rms = 0.; >> 180 >> 181 G4cout.precision(3); >> 182 G4cout >> 183 << "\t Step size= " << G4BestUnit(stepSize,"Length") >> 184 << " +- " << G4BestUnit( rms,"Length") >> 185 << G4endl; >> 186 >> 187 //transmission coefficients >> 188 // >> 189 G4double absorbed = 100.*status[0]/fNbofEvents; >> 190 G4double transmit = 100.*status[1]/fNbofEvents; >> 191 G4double reflected = 100.*status[2]/fNbofEvents; >> 192 >> 193 G4cout.precision(2); >> 194 G4cout >> 195 << "\n absorbed = " << absorbed << " %" >> 196 << " transmit = " << transmit << " %" >> 197 << " reflected = " << reflected << " %" << G4endl; >> 198 >> 199 // normalize histogram of longitudinal energy profile 99 // 200 // 100 G4double stepMax = DBL_MAX; << 101 G4int ih = 1; 201 G4int ih = 1; 102 if (analysisManager->GetH1Activation(ih)) { << 202 G4double binWidth = histoManager->GetBinWidth(ih); 103 stepMax = analysisManager->GetH1Width(ih) << 203 G4double fac = (1./(NbofEvents*binWidth))*(mm/MeV); 104 } << 204 histoManager->Scale(ih,fac); 105 << 205 106 ih = 8; << 206 // reset default formats 107 G4ParticleDefinition* particle = fPrimary->G << 207 G4cout.setf(mode,std::ios::floatfield); 108 if (particle->GetPDGCharge() != 0.) { << 208 G4cout.precision(prec); 109 G4double width = analysisManager->GetH1Wid << 209 110 if (width == 0.) width = 1.; << 111 G4EmCalculator emCalculator; << 112 G4double energy = fPrimary->GetParticleGun << 113 G4int nbOfAbsor = fDetector->GetNbOfAbsor( << 114 for (G4int i = 1; i <= nbOfAbsor; i++) { << 115 G4Material* material = fDetector->GetAbs << 116 G4double newCsdaRange = emCalculator.Get << 117 fRun->SetCsdaRange(i, newCsdaRange); << 118 if (analysisManager->GetH1Activation(ih) << 119 if (i > 1) { << 120 G4double thickness = fDetector->GetAbs << 121 G4double xfrontNorm = fRun->GetXfrontN << 122 G4double csdaRange = fRun->GetCsdaRang << 123 G4double newXfrontNorm = xfrontNorm + << 124 fRun->SetXfrontNorm(i, newXfrontNorm); << 125 } << 126 } << 127 } << 128 fPhysics->GetStepMaxProcess()->SetMaxStep2(s << 129 } << 130 << 131 //....oooOO0OOooo........oooOO0OOooo........oo << 132 << 133 void RunAction::EndOfRunAction(const G4Run*) << 134 { << 135 if (isMaster) fRun->EndOfRun(); << 136 << 137 // save histograms 210 // save histograms 138 G4AnalysisManager* analysisManager = G4Analy << 211 histoManager->save(); 139 if (analysisManager->IsActive()) { << 212 140 analysisManager->Write(); << 141 analysisManager->CloseFile(); << 142 } << 143 << 144 // show Rndm status 213 // show Rndm status 145 if (isMaster) G4Random::showEngineStatus(); << 214 HepRandom::showEngineStatus(); 146 } 215 } 147 216 148 //....oooOO0OOooo........oooOO0OOooo........oo 217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 149 218