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/TestEm9/src/HistoMan 26 /// \file electromagnetic/TestEm9/src/HistoManager.cc 27 /// \brief Implementation of the HistoManager 27 /// \brief Implementation of the HistoManager class 28 // 28 // >> 29 // $Id$ 29 // 30 // 30 //-------------------------------------------- 31 //--------------------------------------------------------------------------- 31 // 32 // 32 // ClassName: HistoManager 33 // ClassName: HistoManager 33 // 34 // 34 // 35 // 35 // Author: V.Ivanchenko 30/01/01 36 // Author: V.Ivanchenko 30/01/01 36 // 37 // 37 //-------------------------------------------- 38 //---------------------------------------------------------------------------- 38 // 39 // 39 40 40 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 41 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 42 43 43 #include "HistoManager.hh" 44 #include "HistoManager.hh" 44 << 45 #include "G4MaterialCutsCouple.hh" 45 #include "EmAcceptance.hh" << 46 #include "Histo.hh" << 47 << 48 #include "G4Electron.hh" << 49 #include "G4EmProcessSubType.hh" 46 #include "G4EmProcessSubType.hh" >> 47 #include "G4VProcess.hh" >> 48 #include "G4VEmProcess.hh" >> 49 #include "G4VEnergyLossProcess.hh" >> 50 #include "G4UnitsTable.hh" >> 51 #include "Histo.hh" >> 52 #include "EmAcceptance.hh" 50 #include "G4Gamma.hh" 53 #include "G4Gamma.hh" 51 #include "G4GammaGeneralProcess.hh" << 54 #include "G4Electron.hh" 52 #include "G4MaterialCutsCouple.hh" << 53 #include "G4Positron.hh" 55 #include "G4Positron.hh" 54 #include "G4SystemOfUnits.hh" 56 #include "G4SystemOfUnits.hh" 55 #include "G4UnitsTable.hh" << 56 #include "G4VEmProcess.hh" << 57 #include "G4VEnergyLossProcess.hh" << 58 #include "G4VProcess.hh" << 59 57 60 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 61 59 62 HistoManager* HistoManager::fManager = nullptr << 60 HistoManager* HistoManager::fManager = 0; 63 61 64 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 65 63 66 HistoManager* HistoManager::GetPointer() 64 HistoManager* HistoManager::GetPointer() 67 { 65 { 68 if (nullptr == fManager) { << 66 if(!fManager) { 69 fManager = new HistoManager(); 67 fManager = new HistoManager(); 70 } 68 } 71 return fManager; 69 return fManager; 72 } 70 } 73 71 74 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 75 73 76 HistoManager::HistoManager() 74 HistoManager::HistoManager() 77 : fGamma(G4Gamma::Gamma()), << 75 : fGamma(0), 78 fElectron(G4Electron::Electron()), << 76 fElectron(0), 79 fPositron(G4Positron::Positron()), << 77 fPositron(0), 80 fHisto(new Histo()) << 78 fHisto(0) 81 { 79 { 82 fVerbose = 1; 80 fVerbose = 1; 83 fEvt1 = -1; << 81 fEvt1 = -1; 84 fEvt2 = -1; << 82 fEvt2 = -1; 85 fNmax = 3; << 83 fNmax = 3; 86 fMaxEnergy = 50.0 * MeV; << 84 fMaxEnergy = 50.0*MeV; 87 fBeamEnergy = 1. * GeV; << 85 fBeamEnergy= 1.*GeV; 88 fMaxEnergyAbs = 10. * MeV; << 86 fMaxEnergyAbs = 10.*MeV; 89 fBinsE = 100; 87 fBinsE = 100; 90 fBinsEA = 40; << 88 fBinsEA= 40; 91 fBinsED = 100; << 89 fBinsED= 100; 92 fNHisto = 13; 90 fNHisto = 13; 93 91 >> 92 fHisto = new Histo(); 94 BookHisto(); 93 BookHisto(); >> 94 >> 95 fGamma = G4Gamma::Gamma(); >> 96 fElectron = G4Electron::Electron(); >> 97 fPositron = G4Positron::Positron(); 95 } 98 } 96 99 97 //....oooOO0OOooo........oooOO0OOooo........oo 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 98 101 99 HistoManager::~HistoManager() 102 HistoManager::~HistoManager() 100 { 103 { 101 delete fHisto; 104 delete fHisto; 102 } 105 } 103 106 104 //....oooOO0OOooo........oooOO0OOooo........oo 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 105 108 106 void HistoManager::BookHisto() 109 void HistoManager::BookHisto() 107 { 110 { 108 fHisto->Add1D("10", "Evis/E0 in central crys << 111 fHisto->Add1D("10","Evis/E0 in central crystal",fBinsED,0.0,1,1.0); 109 fHisto->Add1D("11", "Evis/E0 in 3x3", fBinsE << 112 fHisto->Add1D("11","Evis/E0 in 3x3",fBinsED,0.0,1.0,1.0); 110 fHisto->Add1D("12", "Evis/E0 in 5x5", fBinsE << 113 fHisto->Add1D("12","Evis/E0 in 5x5",fBinsED,0.0,1.0,1.0); 111 fHisto->Add1D("13", "Energy (MeV) of delta-e << 114 fHisto->Add1D("13","Energy (MeV) of delta-electrons", 112 fHisto->Add1D("14", "Energy (MeV) of gammas" << 115 fBinsE,0.0,fMaxEnergy,MeV); 113 fHisto->Add1D("15", "Energy (MeV) in abs1", << 116 fHisto->Add1D("14","Energy (MeV) of gammas",fBinsE,0.0,fMaxEnergy,MeV); 114 fHisto->Add1D("16", "Energy (MeV) in abs2", << 117 fHisto->Add1D("15","Energy (MeV) in abs1",fBinsEA,0.0,fMaxEnergyAbs,MeV); 115 fHisto->Add1D("17", "Energy (MeV) in abs3", << 118 fHisto->Add1D("16","Energy (MeV) in abs2",fBinsEA,0.0,fMaxEnergyAbs,MeV); 116 fHisto->Add1D("18", "Energy (MeV) in abs4", << 119 fHisto->Add1D("17","Energy (MeV) in abs3",fBinsEA,0.0,fMaxEnergyAbs,MeV); 117 fHisto->Add1D("19", "Number of vertex hits", << 120 fHisto->Add1D("18","Energy (MeV) in abs4",fBinsEA,0.0,fMaxEnergyAbs,MeV); 118 fHisto->Add1D("20", "E1/E9 Ratio", fBinsED, << 121 fHisto->Add1D("19","Number of vertex hits",20,-0.5,19.5,1.0); 119 fHisto->Add1D("21", "E1/E25 Ratio", fBinsED, << 122 fHisto->Add1D("20","E1/E9 Ratio",fBinsED,0.0,1,1.0); 120 fHisto->Add1D("22", "E9/E25 Ratio", fBinsED, << 123 fHisto->Add1D("21","E1/E25 Ratio",fBinsED,0.0,1.0,1.0); >> 124 fHisto->Add1D("22","E9/E25 Ratio",fBinsED,0.0,1.0,1.0); 121 } 125 } 122 126 123 //....oooOO0OOooo........oooOO0OOooo........oo 127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 124 128 125 void HistoManager::BeginOfRun() 129 void HistoManager::BeginOfRun() 126 { 130 { 127 // initilise scoring 131 // initilise scoring 128 fEvt = 0; << 132 fEvt = 0; 129 fElec = 0; 133 fElec = 0; 130 fPosit = 0; << 134 fPosit= 0; 131 fGam = 0; << 135 fGam = 0; 132 fStep = 0; 136 fStep = 0; 133 fLowe = 0; 137 fLowe = 0; 134 138 135 for (G4int i = 0; i < 6; i++) { << 139 for(G4int i=0; i<6; i++) { 136 fStat[i] = 0; 140 fStat[i] = 0; 137 fEdep[i] = 0.0; 141 fEdep[i] = 0.0; 138 fErms[i] = 0.0; 142 fErms[i] = 0.0; 139 if (i < 3) { << 143 if(i < 3) { 140 fEdeptr[i] = 0.0; 144 fEdeptr[i] = 0.0; 141 fErmstr[i] = 0.0; 145 fErmstr[i] = 0.0; 142 } 146 } 143 } 147 } 144 148 145 // initialise counters 149 // initialise counters 146 fBrem.resize(93, 0.0); << 150 fBrem.resize(93,0.0); 147 fPhot.resize(93, 0.0); << 151 fPhot.resize(93,0.0); 148 fComp.resize(93, 0.0); << 152 fComp.resize(93,0.0); 149 fConv.resize(93, 0.0); << 153 fConv.resize(93,0.0); 150 154 151 // initialise acceptance - by default is not 155 // initialise acceptance - by default is not applied 152 for (G4int i = 0; i < fNmax; i++) { << 156 for(G4int i=0; i<fNmax; i++) { 153 fEdeptrue[i] = 1.0; 157 fEdeptrue[i] = 1.0; 154 fRmstrue[i] = 1.0; << 158 fRmstrue[i] = 1.0; 155 fLimittrue[i] = 10.; << 159 fLimittrue[i]= 10.; 156 } 160 } 157 161 158 if (fHisto->IsActive()) { << 162 if(fHisto->IsActive()) { 159 for (G4int i = 0; i < fNHisto; ++i) { << 163 for(G4int i=0; i<fNHisto; ++i) {fHisto->Activate(i, true); } 160 fHisto->Activate(i, true); << 161 } << 162 fHisto->Book(); 164 fHisto->Book(); 163 165 164 if (fVerbose > 0) { << 166 if(fVerbose > 0) { 165 G4cout << "HistoManager: Histograms are << 167 G4cout << "HistoManager: Histograms are booked and run has been started" >> 168 << G4endl; 166 } 169 } 167 } 170 } 168 } 171 } 169 172 170 //....oooOO0OOooo........oooOO0OOooo........oo 173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 171 174 172 void HistoManager::EndOfRun(G4int runID) 175 void HistoManager::EndOfRun(G4int runID) 173 { 176 { 174 G4cout << "HistoManager: End of run actions << 177 >> 178 G4cout << "HistoManager: End of run actions are started RunID# " >> 179 << runID << G4endl; 175 G4String nam[6] = {"1x1", "3x3", "5x5", "E1/ 180 G4String nam[6] = {"1x1", "3x3", "5x5", "E1/E9 ", "E1/E25", "E9/E25"}; 176 181 177 // average 182 // average 178 183 179 G4cout << "================================= << 184 G4cout<<"=================================================================" >> 185 <<G4endl; 180 G4double x = (G4double)fEvt; 186 G4double x = (G4double)fEvt; 181 if (fEvt > 0) x = 1.0 / x; << 187 if(fEvt > 0) x = 1.0/x; 182 G4int j; 188 G4int j; 183 for (j = 0; j < fNmax; j++) { << 189 for(j=0; j<fNmax; j++) { >> 190 184 // total mean 191 // total mean 185 fEdep[j] *= x; 192 fEdep[j] *= x; 186 G4double y = fErms[j] * x - fEdep[j] * fEd << 193 G4double y = fErms[j]*x - fEdep[j]*fEdep[j]; 187 if (y < 0.0) y = 0.0; << 194 if(y < 0.0) y = 0.0; 188 fErms[j] = std::sqrt(y); 195 fErms[j] = std::sqrt(y); 189 196 190 // trancated mean 197 // trancated mean 191 G4double xx = G4double(fStat[j]); 198 G4double xx = G4double(fStat[j]); 192 if (xx > 0.0) xx = 1.0 / xx; << 199 if(xx > 0.0) xx = 1.0/xx; 193 fEdeptr[j] *= xx; 200 fEdeptr[j] *= xx; 194 y = fErmstr[j] * xx - fEdeptr[j] * fEdeptr << 201 y = fErmstr[j]*xx - fEdeptr[j]*fEdeptr[j]; 195 if (y < 0.0) y = 0.0; << 202 if(y < 0.0) y = 0.0; 196 fErmstr[j] = std::sqrt(y); 203 fErmstr[j] = std::sqrt(y); 197 } 204 } 198 G4double xe = x * (G4double)fElec; << 205 G4double xe = x*(G4double)fElec; 199 G4double xg = x * (G4double)fGam; << 206 G4double xg = x*(G4double)fGam; 200 G4double xp = x * (G4double)fPosit; << 207 G4double xp = x*(G4double)fPosit; 201 G4double xs = x * fStep; << 208 G4double xs = x*fStep; 202 209 203 G4double f = 100. * std::sqrt(fBeamEnergy / << 210 G4double f = 100.*std::sqrt(fBeamEnergy/GeV); 204 211 205 G4cout << "Number of events " << << 212 G4cout << "Number of events " << fEvt <<G4endl; 206 G4cout << std::setprecision(4) << "Average n 213 G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl; 207 G4cout << std::setprecision(4) << "Average n 214 G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl; 208 G4cout << std::setprecision(4) << "Average n 215 G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl; 209 G4cout << std::setprecision(4) << "Average n 216 G4cout << std::setprecision(4) << "Average number of steps " << xs << G4endl; 210 << 217 211 for (j = 0; j < 3; ++j) { << 218 for(j=0; j<3; j++) { 212 G4double ex = fEdeptr[j]; 219 G4double ex = fEdeptr[j]; 213 G4double sx = fErmstr[j]; 220 G4double sx = fErmstr[j]; 214 G4double xx = G4double(fStat[j]); << 221 G4double xx= G4double(fStat[j]); 215 if (xx > 0.0) xx = 1.0 / xx; << 222 if(xx > 0.0) xx = 1.0/xx; 216 G4double r = sx * std::sqrt(xx); << 223 G4double r = sx*std::sqrt(xx); 217 G4cout << std::setprecision(4) << "Edep " << 224 G4cout << std::setprecision(4) << "Edep " << nam[j] 218 << r; << 225 << " = " << ex 219 if (ex > 0.1) G4cout << " res= " << f * << 226 << " +- " << r; >> 227 if(ex > 0.1) G4cout << " res= " << f*sx/ex << " % " << fStat[j]; 220 G4cout << G4endl; 228 G4cout << G4endl; 221 } 229 } 222 if (fLimittrue[0] < 10. || fLimittrue[1] < 1 << 230 if(fLimittrue[0] < 10. || fLimittrue[1] < 10. || fLimittrue[2] < 10.) { 223 G4cout << "=========== Mean values withou << 231 G4cout 224 for (j = 0; j < fNmax; j++) { << 232 <<"=========== Mean values without trancating ====================="<<G4endl; >> 233 for(j=0; j<fNmax; j++) { 225 G4double ex = fEdep[j]; 234 G4double ex = fEdep[j]; 226 G4double sx = fErms[j]; 235 G4double sx = fErms[j]; 227 G4double rx = sx * std::sqrt(x); << 236 G4double rx = sx*std::sqrt(x); 228 G4cout << std::setprecision(4) << "Edep << 237 G4cout << std::setprecision(4) << "Edep " << nam[j] 229 << rx; << 238 << " = " << ex 230 if (ex > 0.0) G4cout << " res= " << f << 239 << " +- " << rx; >> 240 if(ex > 0.0) G4cout << " res= " << f*sx/ex << " %"; 231 G4cout << G4endl; 241 G4cout << G4endl; 232 } 242 } 233 } 243 } 234 G4cout << "=========== Ratios without tranc << 244 G4cout 235 for (j = 3; j < 6; ++j) { << 245 <<"=========== Ratios without trancating ==========================="<<G4endl; >> 246 for(j=3; j<6; j++) { 236 G4double e = fEdep[j]; 247 G4double e = fEdep[j]; 237 G4double xx = G4double(fStat[j]); << 248 G4double xx= G4double(fStat[j]); 238 if (xx > 0.0) xx = 1.0 / xx; << 249 if(xx > 0.0) xx = 1.0/xx; 239 e *= xx; 250 e *= xx; 240 G4double y = fErms[j] * xx - e * e; << 251 G4double y = fErms[j]*xx - e*e; 241 G4double r = 0.0; 252 G4double r = 0.0; 242 if (y > 0.0) r = std::sqrt(y * xx); << 253 if(y > 0.0) r = std::sqrt(y*xx); 243 G4cout << " " << nam[j] << " = << 254 G4cout << " " << nam[j] << " = " << e >> 255 << " +- " << r; 244 G4cout << G4endl; 256 G4cout << G4endl; 245 } 257 } 246 G4cout << std::setprecision(4) << "Beam Ener << 258 G4cout << std::setprecision(4) << "Beam Energy " << fBeamEnergy/GeV 247 << G4endl; << 259 << " GeV" << G4endl; 248 if (fLowe > 0) G4cout << "Number of events E << 260 if(fLowe > 0) G4cout << "Number of events E/E0<0.8 " << fLowe << G4endl; 249 G4cout << "================================= << 261 G4cout 250 G4cout << G4endl; << 262 <<"=================================================================="<<G4endl; >> 263 G4cout<<G4endl; 251 264 252 // normalise histograms 265 // normalise histograms 253 if (fHisto->IsActive()) { << 266 if(fHisto->IsActive()) { 254 for (G4int i = 0; i < fNHisto; ++i) { << 267 for(G4int i=0; i<fNHisto; i++) { 255 fHisto->ScaleH1(i, x); << 268 fHisto->ScaleH1(i,x); 256 } 269 } 257 fHisto->Save(); 270 fHisto->Save(); 258 } 271 } 259 if (0 < runID) { << 272 if(0 < runID) { return; } 260 return; << 261 } << 262 273 263 // Acceptance only for the first run 274 // Acceptance only for the first run 264 EmAcceptance acc; 275 EmAcceptance acc; 265 G4bool isStarted = false; 276 G4bool isStarted = false; 266 for (j = 0; j < fNmax; j++) { << 277 for (j=0; j<fNmax; j++) { >> 278 267 G4double ltrue = fLimittrue[j]; 279 G4double ltrue = fLimittrue[j]; 268 if (ltrue < DBL_MAX) { 280 if (ltrue < DBL_MAX) { 269 if (!isStarted) { 281 if (!isStarted) { 270 acc.BeginOfAcceptance("Crystal Calorim << 282 acc.BeginOfAcceptance("Crystal Calorimeter",fEvt); 271 isStarted = true; 283 isStarted = true; 272 } 284 } 273 G4double etrue = fEdeptrue[j]; 285 G4double etrue = fEdeptrue[j]; 274 G4double rtrue = fRmstrue[j]; 286 G4double rtrue = fRmstrue[j]; 275 acc.EmAcceptanceGauss("Edep" + nam[j], f << 287 acc.EmAcceptanceGauss("Edep"+nam[j],fEvt,fEdeptr[j],etrue,rtrue,ltrue); 276 acc.EmAcceptanceGauss("Erms" + nam[j], f << 288 acc.EmAcceptanceGauss("Erms"+nam[j],fEvt,fErmstr[j],rtrue,rtrue,2.0*ltrue); 277 } 289 } 278 } 290 } 279 if (isStarted) acc.EndOfAcceptance(); << 291 if(isStarted) acc.EndOfAcceptance(); 280 292 281 // atom frequency 293 // atom frequency 282 G4cout << " Z bremsstrahlung photoeffect 294 G4cout << " Z bremsstrahlung photoeffect compton conversion" << G4endl; 283 for (j = 1; j < 93; ++j) { << 295 for(j=1; j<93; ++j) { 284 G4int n1 = G4int(fBrem[j] * x); << 296 G4int n1 = G4int(fBrem[j]*x); 285 G4int n2 = G4int(fPhot[j] * x); << 297 G4int n2 = G4int(fPhot[j]*x); 286 G4int n3 = G4int(fComp[j] * x); << 298 G4int n3 = G4int(fComp[j]*x); 287 G4int n4 = G4int(fConv[j] * x); << 299 G4int n4 = G4int(fConv[j]*x); 288 if (n1 + n2 + n3 + n4 > 0) { << 300 if(n1 + n2 + n3 + n4 > 0) { 289 G4cout << std::setw(4) << j << std::setw << 301 G4cout << std::setw(4) << j << std::setw(12) << n1 << std::setw(12) << n2 290 << n3 << std::setw(12) << n4 << G << 302 << std::setw(12) << n3 << std::setw(12) << n4 << G4endl; 291 } 303 } 292 } 304 } 293 } 305 } 294 306 295 //....oooOO0OOooo........oooOO0OOooo........oo 307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 296 308 297 void HistoManager::BeginOfEvent() 309 void HistoManager::BeginOfEvent() 298 { 310 { 299 ++fEvt; 311 ++fEvt; 300 312 301 fEabs1 = 0.0; << 313 fEabs1 = 0.0; 302 fEabs2 = 0.0; << 314 fEabs2 = 0.0; 303 fEabs3 = 0.0; << 315 fEabs3 = 0.0; 304 fEabs4 = 0.0; << 316 fEabs4 = 0.0; 305 fEvertex.clear(); 317 fEvertex.clear(); 306 fNvertex.clear(); 318 fNvertex.clear(); 307 for (G4int i = 0; i < 25; i++) { << 319 for (G4int i=0; i<25; i++) { 308 fE[i] = 0.0; 320 fE[i] = 0.0; 309 } 321 } 310 } 322 } 311 323 312 //....oooOO0OOooo........oooOO0OOooo........oo 324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 313 325 314 void HistoManager::EndOfEvent() 326 void HistoManager::EndOfEvent() 315 { 327 { 316 G4double e9 = 0.0; 328 G4double e9 = 0.0; 317 G4double e25 = 0.0; << 329 G4double e25= 0.0; 318 for (G4int i = 0; i < 25; i++) { << 330 for (G4int i=0; i<25; i++) { 319 fE[i] /= fBeamEnergy; 331 fE[i] /= fBeamEnergy; 320 e25 += fE[i]; 332 e25 += fE[i]; 321 if ((6 <= i && 8 >= i) || (11 <= i && 13 > << 333 if( ( 6<=i && 8>=i) || (11<=i && 13>=i) || (16<=i && 18>=i)) e9 += fE[i]; 322 } 334 } 323 335 324 if (1 < fVerbose && e25 < 0.8) { << 336 if(1 < fVerbose && e25 < 0.8) { 325 ++fLowe; 337 ++fLowe; 326 G4cout << "### in the event# " << fEvt << 338 G4cout << "### in the event# " << fEvt << " E25= " << e25 << G4endl; 327 } 339 } 328 340 329 // compute ratios 341 // compute ratios 330 G4double e0 = fE[12]; 342 G4double e0 = fE[12]; 331 G4double e19 = 0.0; << 343 G4double e19 = 0.0; 332 G4double e125 = 0.0; 344 G4double e125 = 0.0; 333 G4double e925 = 0.0; 345 G4double e925 = 0.0; 334 if (e9 > 0.0) { << 346 if(e9 > 0.0) { 335 e19 = e0 / e9; << 347 e19 = e0/e9; 336 e125 = e0 / e25; << 348 e125 = e0/e25; 337 e925 = e9 / e25; << 349 e925 = e9/e25; 338 fEdep[3] += e19; 350 fEdep[3] += e19; 339 fErms[3] += e19 * e19; << 351 fErms[3] += e19*e19; 340 fEdep[4] += e125; 352 fEdep[4] += e125; 341 fErms[4] += e125 * e125; << 353 fErms[4] += e125*e125; 342 fEdep[5] += e925; 354 fEdep[5] += e925; 343 fErms[5] += e925 * e925; << 355 fErms[5] += e925*e925; 344 fStat[3] += 1; 356 fStat[3] += 1; 345 fStat[4] += 1; 357 fStat[4] += 1; 346 fStat[5] += 1; 358 fStat[5] += 1; 347 } 359 } 348 360 349 // Fill histo 361 // Fill histo 350 fHisto->Fill(0, e0, 1.0); << 362 fHisto->Fill(0,e0,1.0); 351 fHisto->Fill(1, e9, 1.0); << 363 fHisto->Fill(1,e9,1.0); 352 fHisto->Fill(2, e25, 1.0); << 364 fHisto->Fill(2,e25,1.0); 353 fHisto->Fill(5, fEabs1, 1.0); << 365 fHisto->Fill(5,fEabs1,1.0); 354 fHisto->Fill(6, fEabs2, 1.0); << 366 fHisto->Fill(6,fEabs2,1.0); 355 fHisto->Fill(7, fEabs3, 1.0); << 367 fHisto->Fill(7,fEabs3,1.0); 356 fHisto->Fill(8, fEabs4, 1.0); << 368 fHisto->Fill(8,fEabs4,1.0); 357 fHisto->Fill(9, G4double(fNvertex.size()), 1 << 369 fHisto->Fill(9,G4double(fNvertex.size()),1.0); 358 fHisto->Fill(10, e19, 1.0); << 370 fHisto->Fill(10,e19,1.0); 359 fHisto->Fill(11, e125, 1.0); << 371 fHisto->Fill(11,e125,1.0); 360 fHisto->Fill(12, e925, 1.0); << 372 fHisto->Fill(12,e925,1.0); 361 373 362 // compute sums 374 // compute sums 363 fEdep[0] += e0; 375 fEdep[0] += e0; 364 fErms[0] += e0 * e0; << 376 fErms[0] += e0*e0; 365 fEdep[1] += e9; 377 fEdep[1] += e9; 366 fErms[1] += e9 * e9; << 378 fErms[1] += e9*e9; 367 fEdep[2] += e25; 379 fEdep[2] += e25; 368 fErms[2] += e25 * e25; << 380 fErms[2] += e25*e25; 369 381 370 // trancated mean 382 // trancated mean 371 if (std::abs(e0 - fEdeptrue[0]) < fRmstrue[0 << 383 if(std::abs(e0-fEdeptrue[0])<fRmstrue[0]*fLimittrue[0]) { 372 fStat[0] += 1; 384 fStat[0] += 1; 373 fEdeptr[0] += e0; 385 fEdeptr[0] += e0; 374 fErmstr[0] += e0 * e0; << 386 fErmstr[0] += e0*e0; 375 } 387 } 376 if (std::abs(e9 - fEdeptrue[1]) < fRmstrue[1 << 388 if(std::abs(e9-fEdeptrue[1])<fRmstrue[1]*fLimittrue[1]) { 377 fStat[1] += 1; 389 fStat[1] += 1; 378 fEdeptr[1] += e9; 390 fEdeptr[1] += e9; 379 fErmstr[1] += e9 * e9; << 391 fErmstr[1] += e9*e9; 380 } 392 } 381 if (std::abs(e25 - fEdeptrue[2]) < fRmstrue[ << 393 if(std::abs(e25-fEdeptrue[2])<fRmstrue[2]*fLimittrue[2]) { 382 fStat[2] += 1; 394 fStat[2] += 1; 383 fEdeptr[2] += e25; 395 fEdeptr[2] += e25; 384 fErmstr[2] += e25 * e25; << 396 fErmstr[2] += e25*e25; 385 } 397 } 386 } 398 } 387 399 388 //....oooOO0OOooo........oooOO0OOooo........oo 400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 389 401 390 void HistoManager::ScoreNewTrack(const G4Track 402 void HistoManager::ScoreNewTrack(const G4Track* aTrack) 391 { 403 { 392 // Save primary parameters << 404 //Save primary parameters 393 ResetTrackLength(); 405 ResetTrackLength(); 394 const G4ParticleDefinition* particle = aTrac 406 const G4ParticleDefinition* particle = aTrack->GetDefinition(); 395 const G4DynamicParticle* dynParticle = aTrac 407 const G4DynamicParticle* dynParticle = aTrack->GetDynamicParticle(); 396 408 397 G4int pid = aTrack->GetParentID(); 409 G4int pid = aTrack->GetParentID(); 398 G4double kinE = dynParticle->GetKineticEnerg 410 G4double kinE = dynParticle->GetKineticEnergy(); 399 G4ThreeVector pos = aTrack->GetVertexPositio 411 G4ThreeVector pos = aTrack->GetVertexPosition(); 400 412 401 // primary 413 // primary 402 if (0 == pid) { << 414 if(0 == pid) { >> 415 403 G4double mass = 0.0; 416 G4double mass = 0.0; 404 if (particle) { << 417 if(particle) { 405 mass = particle->GetPDGMass(); 418 mass = particle->GetPDGMass(); 406 } 419 } 407 420 408 G4ThreeVector dir = dynParticle->GetMoment 421 G4ThreeVector dir = dynParticle->GetMomentumDirection(); 409 if (1 < fVerbose) { << 422 if(1 < fVerbose) { 410 G4cout << "TrackingAction: Primary kinE( << 423 G4cout << "TrackingAction: Primary kinE(MeV)= " << kinE/MeV 411 << "; pos= " << pos << "; dir= " << 424 << "; m(MeV)= " << mass/MeV >> 425 << "; pos= " << pos << "; dir= " << dir << G4endl; 412 } 426 } 413 427 414 // secondary 428 // secondary 415 } << 429 } else { 416 else { << 417 const G4VProcess* proc = aTrack->GetCreato 430 const G4VProcess* proc = aTrack->GetCreatorProcess(); 418 G4int type = proc->GetProcessSubType(); 431 G4int type = proc->GetProcessSubType(); 419 << 432 if(type == fBremsstrahlung) { 420 if (type == fBremsstrahlung) { << 433 const G4Element* elm = 421 auto elm = static_cast<const G4VEnergyLo << 434 static_cast<const G4VEnergyLossProcess*>(proc)->GetCurrentElement(); 422 if (nullptr != elm) { << 435 if(elm) { 423 G4int Z = elm->GetZasInt(); << 436 G4int Z = G4lrint(elm->GetZ()); 424 if (Z > 0 && Z < 93) { << 437 if(Z > 0 && Z < 93) { fBrem[Z] += 1.0; } 425 fBrem[Z] += 1.0; << 426 } << 427 } 438 } 428 } << 439 } else if(type == fPhotoElectricEffect) { 429 else if (type == fPhotoElectricEffect) { << 440 const G4Element* elm = 430 auto elm = static_cast<const G4VEmProces << 441 static_cast<const G4VEmProcess*>(proc)->GetCurrentElement(); 431 if (nullptr != elm) { << 442 if(elm) { 432 G4int Z = elm->GetZasInt(); << 443 G4int Z = G4lrint(elm->GetZ()); 433 if (Z > 0 && Z < 93) { << 444 if(Z > 0 && Z < 93) { fPhot[Z] += 1.0; } 434 fPhot[Z] += 1.0; << 435 } << 436 } 445 } 437 } << 446 } else if(type == fGammaConversion) { 438 else if (type == fGammaConversion) { << 447 const G4Element* elm = 439 auto elm = static_cast<const G4VEmProces << 448 static_cast<const G4VEmProcess*>(proc)->GetCurrentElement(); 440 if (nullptr != elm) { << 449 if(elm) { 441 G4int Z = elm->GetZasInt(); << 450 G4int Z = G4lrint(elm->GetZ()); 442 if (Z > 0 && Z < 93) { << 451 if(Z > 0 && Z < 93) { fConv[Z] += 1.0; } 443 fConv[Z] += 1.0; << 444 } << 445 } 452 } 446 } << 453 } else if(type == fComptonScattering) { 447 else if (type == fComptonScattering) { << 454 const G4Element* elm = 448 auto elm = static_cast<const G4VEmProces << 455 static_cast<const G4VEmProcess*>(proc)->GetCurrentElement(); 449 if (nullptr != elm) { << 456 if(elm) { 450 G4int Z = elm->GetZasInt(); << 457 G4int Z = G4lrint(elm->GetZ()); 451 if (Z > 0 && Z < 93) { << 458 if(Z > 0 && Z < 93) { fComp[Z] += 1.0; } 452 fComp[Z] += 1.0; << 453 } << 454 } 459 } 455 } 460 } 456 461 457 // delta-electron 462 // delta-electron 458 if (particle == fElectron) { 463 if (particle == fElectron) { 459 if (1 < fVerbose) { << 464 if(1 < fVerbose) { 460 G4cout << "TrackingAction: Secondary e 465 G4cout << "TrackingAction: Secondary electron " << G4endl; 461 } 466 } 462 AddDeltaElectron(dynParticle); 467 AddDeltaElectron(dynParticle); 463 } << 468 464 else if (particle == fPositron) { << 469 } else if (particle == fPositron) { 465 if (1 < fVerbose) { << 470 if(1 < fVerbose) { 466 G4cout << "TrackingAction: Secondary p 471 G4cout << "TrackingAction: Secondary positron " << G4endl; 467 } 472 } 468 AddPositron(dynParticle); 473 AddPositron(dynParticle); 469 } << 474 470 else if (particle == fGamma) { << 475 } else if (particle == fGamma) { 471 if (1 < fVerbose) { << 476 if(1 < fVerbose) { 472 G4cout << "TrackingAction: Secondary g << 477 G4cout << "TrackingAction: Secondary gamma; parentID= " << pid >> 478 << " E= " << kinE << G4endl; 473 } 479 } 474 AddPhoton(dynParticle); 480 AddPhoton(dynParticle); 475 } 481 } 476 } 482 } 477 } 483 } 478 484 479 //....oooOO0OOooo........oooOO0OOooo........oo 485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 480 486 481 void HistoManager::AddEnergy(G4double edep, G4 487 void HistoManager::AddEnergy(G4double edep, G4int volIndex, G4int copyNo) 482 { 488 { 483 if (1 < fVerbose) { << 489 if(1 < fVerbose) { 484 G4cout << "HistoManager::AddEnergy: e(keV) << 490 G4cout << "HistoManager::AddEnergy: e(keV)= " << edep/keV 485 << "; copyNo= " << copyNo << G4endl << 491 << "; volIdx= " << volIndex >> 492 << "; copyNo= " << copyNo >> 493 << G4endl; 486 } 494 } 487 if (0 == volIndex) { << 495 if(0 == volIndex) { 488 fE[copyNo] += edep; 496 fE[copyNo] += edep; 489 } << 497 } else if (1 == volIndex) { 490 else if (1 == volIndex) { << 491 fEabs1 += edep; 498 fEabs1 += edep; 492 } << 499 } else if (2 == volIndex) { 493 else if (2 == volIndex) { << 494 fEabs2 += edep; 500 fEabs2 += edep; 495 } << 501 } else if (3 == volIndex) { 496 else if (3 == volIndex) { << 497 fEabs3 += edep; 502 fEabs3 += edep; 498 } << 503 } else if (4 == volIndex) { 499 else if (4 == volIndex) { << 500 fEabs4 += edep; 504 fEabs4 += edep; 501 } << 505 } else if (5 == volIndex) { 502 else if (5 == volIndex) { << 503 G4int n = fNvertex.size(); 506 G4int n = fNvertex.size(); 504 G4bool newPad = true; 507 G4bool newPad = true; 505 if (n > 0) { 508 if (n > 0) { 506 for (G4int i = 0; i < n; i++) { << 509 for(G4int i=0; i<n; i++) { 507 if (copyNo == fNvertex[i]) { 510 if (copyNo == fNvertex[i]) { 508 newPad = false; 511 newPad = false; 509 fEvertex[i] += edep; 512 fEvertex[i] += edep; 510 break; 513 break; 511 } 514 } 512 } 515 } 513 } 516 } 514 if (newPad) { << 517 if(newPad) { 515 fNvertex.push_back(copyNo); 518 fNvertex.push_back(copyNo); 516 fEvertex.push_back(edep); 519 fEvertex.push_back(edep); 517 } 520 } 518 } 521 } 519 } 522 } 520 523 521 //....oooOO0OOooo........oooOO0OOooo........oo 524 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 522 525 523 void HistoManager::AddDeltaElectron(const G4Dy 526 void HistoManager::AddDeltaElectron(const G4DynamicParticle* elec) 524 { 527 { 525 G4double e = elec->GetKineticEnergy() / MeV; << 528 G4double e = elec->GetKineticEnergy()/MeV; 526 if (e > 0.0) { << 529 if(e > 0.0) fElec++; 527 ++fElec; << 530 fHisto->Fill(3,e,1.0); 528 fHisto->Fill(3, e, 1.0); << 529 } << 530 } 531 } 531 532 532 //....oooOO0OOooo........oooOO0OOooo........oo 533 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 533 534 534 void HistoManager::AddPhoton(const G4DynamicPa 535 void HistoManager::AddPhoton(const G4DynamicParticle* ph) 535 { 536 { 536 G4double e = ph->GetKineticEnergy() / MeV; << 537 G4double e = ph->GetKineticEnergy()/MeV; 537 if (e > 0.0) { << 538 if(e > 0.0) fGam++; 538 ++fGam; << 539 fHisto->Fill(4,e,1.0); 539 fHisto->Fill(4, e, 1.0); << 540 } << 541 } 540 } 542 541 543 //....oooOO0OOooo........oooOO0OOooo........oo 542 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 544 543 545 void HistoManager::SetEdepAndRMS(G4int i, cons << 544 void HistoManager::SetEdepAndRMS(G4int i, G4ThreeVector val) 546 { 545 { 547 if (i < fNmax && i >= 0) { << 546 if(i<fNmax && i>=0) { 548 if (val[0] > 0.0) fEdeptrue[i] = val[0]; << 547 if(val[0] > 0.0) fEdeptrue[i] = val[0]; 549 if (val[1] > 0.0) fRmstrue[i] = val[1]; << 548 if(val[1] > 0.0) fRmstrue[i] = val[1]; 550 if (val[2] > 0.0) fLimittrue[i] = val[2]; << 549 if(val[2] > 0.0) fLimittrue[i] = val[2]; 551 } 550 } 552 } 551 } 553 552 554 //....oooOO0OOooo........oooOO0OOooo........oo 553 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 554 555 555