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