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/TestEm7/src/RunActio 26 /// \file electromagnetic/TestEm7/src/RunAction.cc 27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction class 28 // 28 // 29 // << 29 // $Id$ >> 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 #include "DetectorConstruction.hh" 35 #include "DetectorConstruction.hh" 36 #include "PhysicsList.hh" 36 #include "PhysicsList.hh" 37 #include "PrimaryGeneratorAction.hh" << 38 #include "StepMax.hh" 37 #include "StepMax.hh" >> 38 #include "PrimaryGeneratorAction.hh" 39 39 40 #include "G4Run.hh" 40 #include "G4Run.hh" 41 #include "G4RunManager.hh" 41 #include "G4RunManager.hh" 42 #include "G4SystemOfUnits.hh" << 43 #include "G4UnitsTable.hh" 42 #include "G4UnitsTable.hh" >> 43 #include "G4SystemOfUnits.hh" 44 #include "G4ios.hh" 44 #include "G4ios.hh" >> 45 45 #include "Randomize.hh" 46 #include "Randomize.hh" 46 47 47 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 48 49 49 RunAction::RunAction(DetectorConstruction* det << 50 RunAction::RunAction(DetectorConstruction* det, PhysicsList* phys, 50 : G4UserRunAction(), << 51 PrimaryGeneratorAction* kin) 51 fAnalysisManager(0), << 52 : fAnalysisManager(0), fDetector(det), fPhysics(phys), fKinematic(kin), 52 fDetector(det), << 53 fTallyEdep(new G4double[MaxTally]), fProjRange(0.), fProjRange2(0.), 53 fPhysics(phys), << 54 fEdeptot(0.), fEniel(0.), fNbPrimarySteps(0), fRange(0) 54 fKinematic(kin), << 55 { 55 fTallyEdep(new G4double[kMaxTally]), << 56 fProjRange(0.), << 57 fProjRange2(0.), << 58 fEdeptot(0.), << 59 fEniel(0.), << 60 fNbPrimarySteps(0), << 61 fRange(0) << 62 { << 63 // Book predefined histograms 56 // Book predefined histograms 64 BookHisto(); << 57 BookHisto(); 65 } 58 } 66 59 67 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 68 61 69 RunAction::~RunAction() 62 RunAction::~RunAction() 70 { 63 { 71 delete[] fTallyEdep; << 64 delete [] fTallyEdep; 72 } 65 } 73 66 74 //....oooOO0OOooo........oooOO0OOooo........oo 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 75 68 76 void RunAction::BeginOfRunAction(const G4Run* 69 void RunAction::BeginOfRunAction(const G4Run* aRun) 77 { << 70 { 78 G4cout << "### Run " << aRun->GetRunID() << 71 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl; 79 << 72 80 if (!fAnalysisManager) { << 73 // save Rndm status 81 BookHisto(); << 74 ////G4RunManager::GetRunManager()->SetRandomNumberStore(true); 82 } << 83 << 84 CLHEP::HepRandom::showEngineStatus(); 75 CLHEP::HepRandom::showEngineStatus(); 85 << 76 86 // initialize projected range, tallies, Ebea << 77 //initialize projected range, tallies, Ebeam, and book histograms 87 // 78 // 88 fNbPrimarySteps = 0; 79 fNbPrimarySteps = 0; 89 fRange = 0; 80 fRange = 0; 90 fProjRange = fProjRange2 = 0.; 81 fProjRange = fProjRange2 = 0.; 91 fEdeptot = fEniel = 0.; 82 fEdeptot = fEniel = 0.; 92 for (G4int j = 0; j < kMaxTally; ++j) { << 83 for (G4int j=0; j<MaxTally; j++) fTallyEdep[j] = 0.; 93 fTallyEdep[j] = 0.; << 94 } << 95 fKinematic->ResetEbeamCumul(); 84 fKinematic->ResetEbeamCumul(); 96 << 85 97 if (fAnalysisManager->IsActive()) { 86 if (fAnalysisManager->IsActive()) { 98 fAnalysisManager->OpenFile(); << 87 fAnalysisManager->OpenFile(); 99 88 100 // histogram "1" is defined by the length 89 // histogram "1" is defined by the length of the target 101 // zoomed histograms are defined by UI com << 90 // zoomed histograms are defined by UI command 102 G4double length = fDetector->GetAbsorSizeX << 91 G4double length = fDetector->GetAbsorSizeX(); 103 G4double stepMax = fPhysics->GetStepMaxPro 92 G4double stepMax = fPhysics->GetStepMaxProcess()->GetMaxStep(); 104 G4int nbBins = 100; << 93 G4int nbmin = 100; 105 if (stepMax < DBL_MAX) { << 94 G4int nbBins = (G4int)(0.5 + length/stepMax); 106 G4int nb = (G4int)(0.5 + length / stepMa << 95 if (nbBins < nbmin) nbBins = nbmin; 107 nbBins = std::min(std::max(nbBins, nb), << 108 } << 109 fAnalysisManager->SetH1(1, nbBins, 0., len 96 fAnalysisManager->SetH1(1, nbBins, 0., length, "mm"); 110 } 97 } 111 } 98 } 112 99 113 //....oooOO0OOooo........oooOO0OOooo........oo 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 114 101 115 void RunAction::EndOfRunAction(const G4Run* aR 102 void RunAction::EndOfRunAction(const G4Run* aRun) 116 { 103 { 117 G4int nbofEvents = aRun->GetNumberOfEvent(); 104 G4int nbofEvents = aRun->GetNumberOfEvent(); 118 if (nbofEvents == 0) return; 105 if (nbofEvents == 0) return; 119 106 120 // run conditions << 107 //run conditions 121 // << 108 // 122 const G4Material* material = fDetector->GetA << 109 G4Material* material = fDetector->GetAbsorMaterial(); 123 G4double density = material->GetDensity(); 110 G4double density = material->GetDensity(); 124 << 111 125 G4String particle = fKinematic->GetParticleG << 112 G4String particle = fKinematic->GetParticleGun()->GetParticleDefinition() >> 113 ->GetParticleName(); 126 G4double energy = fKinematic->GetParticleGun 114 G4double energy = fKinematic->GetParticleGun()->GetParticleEnergy(); 127 G4cout << "\n The run consists of " << nbofE << 115 G4cout << "\n The run consists of " << nbofEvents << " "<< particle << " of " 128 << G4BestUnit(energy, "Energy") << " << 116 << G4BestUnit(energy,"Energy") << " through " 129 << G4BestUnit(fDetector->GetAbsorSize << 117 << G4BestUnit(fDetector->GetAbsorSizeX(),"Length") << " of " 130 << " (density: " << G4BestUnit(densit << 118 << material->GetName() << " (density: " 131 << 119 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl; 132 // compute projected range and straggling << 120 >> 121 //compute projected range and straggling 133 // 122 // 134 if (fRange > 0) { << 123 if(fRange > 0) { 135 fProjRange /= fRange; << 124 fProjRange /= fRange; 136 fProjRange2 /= fRange; 125 fProjRange2 /= fRange; 137 } 126 } 138 G4double rms = fProjRange2 - fProjRange * fP << 127 G4double rms = fProjRange2 - fProjRange*fProjRange; 139 if (rms > 0.) << 128 if (rms>0.) rms = std::sqrt(rms); else rms = 0.; 140 rms = std::sqrt(rms); << 129 141 else << 130 G4double nstep = G4double(fNbPrimarySteps)/G4double(nbofEvents); 142 rms = 0.; << 143 << 144 G4double nstep = G4double(fNbPrimarySteps) / << 145 << 146 G4cout.precision(6); << 147 G4cout << "\n Projected Range= " << G4BestUn << 148 << " rms= " << G4BestUnit(rms, "Len << 149 G4cout << " Mean number of primary steps = " << 150 << 151 // compute energy deposition and niel << 152 // << 153 fEdeptot /= nbofEvents; << 154 G4cout << " Total energy deposit= " << G4Bes << 155 fEniel /= nbofEvents; << 156 G4cout << " niel energy deposit = " << G4Bes << 157 131 158 // print dose in tallies << 132 G4cout.precision(6); >> 133 G4cout << "\n Projected Range= "<< G4BestUnit(fProjRange,"Length") >> 134 << " rms= " << G4BestUnit( rms,"Length") >> 135 << G4endl; >> 136 G4cout << " Mean number of primary steps = "<< nstep << G4endl; >> 137 >> 138 //compute energy deposition and niel >> 139 // >> 140 fEdeptot /= nbofEvents; >> 141 G4cout << " Total energy deposit= "<< G4BestUnit(fEdeptot,"Energy") >> 142 << G4endl; >> 143 fEniel /= nbofEvents; >> 144 G4cout << " niel energy deposit = "<< G4BestUnit(fEniel,"Energy") >> 145 << G4endl; >> 146 >> 147 //print dose in tallies 159 // 148 // 160 G4int tallyNumber = fDetector->GetTallyNumbe 149 G4int tallyNumber = fDetector->GetTallyNumber(); 161 if (tallyNumber > 0) { 150 if (tallyNumber > 0) { 162 G4double Ebeam = fKinematic->GetEbeamCumul 151 G4double Ebeam = fKinematic->GetEbeamCumul(); 163 G4cout << "\n----------------------------- 152 G4cout << "\n---------------------------------------------------------\n"; 164 G4cout << " Cumulated Doses : \tEdep 153 G4cout << " Cumulated Doses : \tEdep \tEdep/Ebeam \tDose" << G4endl; 165 for (G4int j = 0; j < tallyNumber; ++j) { << 154 for (G4int j=1; j <= tallyNumber; j++) { 166 G4double Edep = fTallyEdep[j], ratio = 1 << 155 G4double Edep = fTallyEdep[j], ratio = 100*Edep/Ebeam; 167 G4double tallyMass = fDetector->GetTally << 156 G4double tallyMass = fDetector->GetTallyMass(j); 168 G4double Dose = Edep / tallyMass; << 157 G4double Dose = Edep/tallyMass; 169 G4cout << " tally " << j << ": \t \t" << << 158 G4cout << " tally " << j << ": \t \t" 170 << " % \t" << G4BestUnit(Dose, "D << 159 << G4BestUnit(Edep,"Energy") << "\t" >> 160 << ratio << " % \t" >> 161 << G4BestUnit(Dose,"Dose") << G4endl; 171 } 162 } 172 G4cout << "\n----------------------------- 163 G4cout << "\n---------------------------------------------------------\n"; 173 G4cout << G4endl; << 164 G4cout << G4endl; 174 } 165 } 175 166 176 if (fAnalysisManager->IsActive()) { << 167 // normalize histograms 177 // normalize histograms << 168 // 178 // << 169 for (G4int j=1; j<3; j++) { 179 for (G4int j = 1; j < 3; ++j) { << 170 G4double binWidth = fAnalysisManager->GetH1Width(j); 180 G4double binWidth = fAnalysisManager->Ge << 171 G4double fac = (mm/MeV)/(nbofEvents * binWidth); 181 G4double fac = (mm / MeV) / (nbofEvents << 172 fAnalysisManager->ScaleH1(j, fac); 182 fAnalysisManager->ScaleH1(j, fac); << 173 } 183 } << 174 fAnalysisManager->ScaleH1(3, 1./nbofEvents); 184 fAnalysisManager->ScaleH1(3, 1. / nbofEven << 175 185 << 176 // save histograms 186 // save histograms << 177 if (fAnalysisManager->IsActive() ) { 187 fAnalysisManager->Write(); 178 fAnalysisManager->Write(); 188 fAnalysisManager->CloseFile(); 179 fAnalysisManager->CloseFile(); 189 } 180 } 190 << 181 191 // show Rndm status 182 // show Rndm status 192 // 183 // 193 CLHEP::HepRandom::showEngineStatus(); 184 CLHEP::HepRandom::showEngineStatus(); 194 } 185 } 195 186 196 //....oooOO0OOooo........oooOO0OOooo........oo 187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 197 188 198 void RunAction::BookHisto() 189 void RunAction::BookHisto() 199 { 190 { 200 // Create or get analysis manager 191 // Create or get analysis manager 201 // The choice of analysis technology is done 192 // The choice of analysis technology is done via selection of a namespace 202 // in HistoManager.hh 193 // in HistoManager.hh 203 fAnalysisManager = G4AnalysisManager::Instan 194 fAnalysisManager = G4AnalysisManager::Instance(); 204 fAnalysisManager->SetDefaultFileType("root") << 205 fAnalysisManager->SetFileName("testem7"); 195 fAnalysisManager->SetFileName("testem7"); 206 fAnalysisManager->SetVerboseLevel(1); 196 fAnalysisManager->SetVerboseLevel(1); 207 fAnalysisManager->SetActivation(true); // e 197 fAnalysisManager->SetActivation(true); // enable inactivation of histograms 208 198 209 // Define histograms start values 199 // Define histograms start values 210 const G4int kMaxHisto = 4; 200 const G4int kMaxHisto = 4; 211 const G4String id[] = {"h0", "h1", "h2", "h3 << 201 const G4String id[] = { "0", "1", "2", "3" }; 212 const G4String title[] = { << 202 const G4String title[] = 213 "dummy", // 0 << 203 { "dummy", //0 214 "Edep (MeV/mm) along absorber ", // 1 << 204 "Edep (MeV/mm) along absorber ", //1 215 "Edep (MeV/mm) along absorber zoomed", // << 205 "Edep (MeV/mm) along absorber zoomed", //2 216 "projectile range" // 3 << 206 "projectile range" //3 217 }; << 207 }; 218 208 219 // Default values (to be reset via /analysis << 209 // Default values (to be reset via /analysis/h1/set command) 220 G4int nbins = 100; 210 G4int nbins = 100; 221 G4double vmin = 0.; 211 G4double vmin = 0.; 222 G4double vmax = 100.; 212 G4double vmax = 100.; 223 213 224 // Create all histograms as inactivated << 214 // Create all histograms as inactivated 225 // as we have not yet set nbins, vmin, vmax 215 // as we have not yet set nbins, vmin, vmax 226 for (G4int k = 0; k < kMaxHisto; ++k) { << 216 for (G4int k=0; k<kMaxHisto; k++) { 227 G4int ih = fAnalysisManager->CreateH1(id[k 217 G4int ih = fAnalysisManager->CreateH1(id[k], title[k], nbins, vmin, vmax); 228 G4bool activ = false; 218 G4bool activ = false; 229 if (k == 1) activ = true; 219 if (k == 1) activ = true; 230 fAnalysisManager->SetH1Activation(ih, acti << 220 fAnalysisManager->SetActivation(G4VAnalysisManager::kH1, ih, activ); 231 } 221 } 232 } 222 } 233 223 234 //....oooOO0OOooo........oooOO0OOooo........oo 224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 235 225