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 hadronic/Hadr01/src/HistoManager.cc << 26 // $Id: HistoManager.cc,v 1.5 2006/06/29 17:24:16 gunter Exp $ 27 /// \brief Implementation of the HistoManager << 27 // GEANT4 tag $Name: geant4-08-01-patch-01 $ 28 // 28 // 29 //-------------------------------------------- 29 //--------------------------------------------------------------------------- 30 // 30 // 31 // ClassName: HistoManager 31 // ClassName: HistoManager 32 // 32 // 33 // 33 // 34 // Author: V.Ivanchenko 30/01/01 34 // Author: V.Ivanchenko 30/01/01 35 // 35 // 36 // Modified: 36 // Modified: 37 // 04.06.2006 Adoptation of Hadr01 (V.Ivanchen << 37 // 04.06.2006 Adoptation of hadr01 (V.Ivanchenko) 38 // 16.11.2006 Add beamFlag (V.Ivanchenko) << 39 // 38 // 40 //-------------------------------------------- 39 //---------------------------------------------------------------------------- 41 // 40 // 42 41 43 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 44 //....oooOO0OOooo........oooOO0OOooo........oo 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 45 44 46 #include "HistoManager.hh" 45 #include "HistoManager.hh" 47 << 46 #include "G4UnitsTable.hh" 48 #include "Histo.hh" << 47 #include "G4Neutron.hh" 49 << 48 #include "G4Proton.hh" 50 #include "G4Alpha.hh" << 51 #include "G4AntiProton.hh" 49 #include "G4AntiProton.hh" 52 #include "G4Deuteron.hh" << 53 #include "G4Electron.hh" << 54 #include "G4Gamma.hh" 50 #include "G4Gamma.hh" 55 #include "G4He3.hh" << 51 #include "G4Electron.hh" 56 #include "G4KaonMinus.hh" << 52 #include "G4Positron.hh" 57 #include "G4KaonPlus.hh" << 58 #include "G4KaonZeroLong.hh" << 59 #include "G4KaonZeroShort.hh" << 60 #include "G4MuonMinus.hh" << 61 #include "G4MuonPlus.hh" 53 #include "G4MuonPlus.hh" 62 #include "G4Neutron.hh" << 54 #include "G4MuonMinus.hh" 63 #include "G4PhysicalConstants.hh" << 64 #include "G4PionMinus.hh" << 65 #include "G4PionPlus.hh" 55 #include "G4PionPlus.hh" >> 56 #include "G4PionMinus.hh" 66 #include "G4PionZero.hh" 57 #include "G4PionZero.hh" 67 #include "G4Positron.hh" << 58 #include "G4KaonPlus.hh" 68 #include "G4Proton.hh" << 59 #include "G4KaonMinus.hh" 69 #include "G4Step.hh" << 60 #include "G4KaonZeroShort.hh" 70 #include "G4SystemOfUnits.hh" << 61 #include "G4KaonZeroLong.hh" 71 #include "G4Track.hh" << 62 #include "G4Deuteron.hh" 72 #include "G4Triton.hh" 63 #include "G4Triton.hh" 73 #include "G4UnitsTable.hh" << 64 #include "G4He3.hh" >> 65 #include "G4Alpha.hh" >> 66 #include "Histo.hh" 74 #include "globals.hh" 67 #include "globals.hh" 75 68 76 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 77 70 78 HistoManager* HistoManager::fManager = nullptr << 71 HistoManager* HistoManager::fManager = 0; 79 72 80 //....oooOO0OOooo........oooOO0OOooo........oo 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 81 74 82 HistoManager* HistoManager::GetPointer() 75 HistoManager* HistoManager::GetPointer() 83 { 76 { 84 if (nullptr == fManager) { << 77 if(!fManager) { 85 static HistoManager manager; 78 static HistoManager manager; 86 fManager = &manager; 79 fManager = &manager; 87 } 80 } 88 return fManager; 81 return fManager; 89 } 82 } 90 83 91 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 92 85 93 HistoManager::HistoManager() 86 HistoManager::HistoManager() 94 : fPrimaryDef(nullptr), << 95 fRadius(10 * CLHEP::cm), << 96 fLength(300. * CLHEP::mm), << 97 fEdepMax(1.0 * CLHEP::GeV) << 98 { 87 { 99 fHisto = new Histo(); << 88 verbose= 0; 100 fHisto->SetVerbose(fVerbose); << 89 nSlices = 300; 101 fNeutron = G4Neutron::Neutron(); << 90 nBinsE = 100; >> 91 nHisto = 19; >> 92 length = 300.*mm; >> 93 histo = new Histo(verbose); >> 94 neutron = G4Neutron::Neutron(); 102 } 95 } 103 96 104 //....oooOO0OOooo........oooOO0OOooo........oo 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 105 98 106 HistoManager::~HistoManager() 99 HistoManager::~HistoManager() 107 { 100 { 108 delete fHisto; << 101 delete histo; 109 } 102 } 110 103 111 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 112 105 113 void HistoManager::BookHisto() << 106 void HistoManager::bookHisto() 114 { 107 { 115 fHistoBooked = true; << 108 histo->add1D("1", 116 fHisto->Add1D("1", "Energy deposition (MeV/m << 109 "Energy deposition (MeV/mm/event) in the target",nSlices,0.0,length/mm,MeV/mm); 117 MeV / mm); << 110 118 fHisto->Add1D("2", "Log10 Energy (MeV) of ga << 111 histo->add1D("2", 119 fHisto->Add1D("3", "Log10 Energy (MeV) of el << 112 "Log10 Energy (MeV) of gammas",nBinsE,-5.,5.,1.0); 120 fHisto->Add1D("4", "Log10 Energy (MeV) of po << 113 121 fHisto->Add1D("5", "Log10 Energy (MeV) of pr << 114 histo->add1D("3", 122 fHisto->Add1D("6", "Log10 Energy (MeV) of ne << 115 "Log10 Energy (MeV) of electrons",nBinsE,-5.,5.,1.0); 123 fHisto->Add1D("7", "Log10 Energy (MeV) of ch << 116 124 fHisto->Add1D("8", "Log10 Energy (MeV) of pi << 117 histo->add1D("4", 125 fHisto->Add1D("9", "Log10 Energy (MeV) of ch << 118 "Log10 Energy (MeV) of positrons",nBinsE,-5.,5.,1.0); 126 fHisto->Add1D("10", "Log10 Energy (MeV) of n << 119 127 fHisto->Add1D("11", "Log10 Energy (MeV) of d << 120 histo->add1D("5", 128 fHisto->Add1D("12", "Log10 Energy (MeV) of H << 121 "Log10 Energy (MeV) of protons",nBinsE,-5.,5.,1.0); 129 fHisto->Add1D("13", "Log10 Energy (MeV) of G << 122 130 fHisto->Add1D("14", "Log10 Energy (MeV) of m << 123 histo->add1D("6", 131 fHisto->Add1D("15", "log10 Energy (MeV) of s << 124 "Log10 Energy (MeV) of neutrons",nBinsE,-5.,5.,1.0); 132 fHisto->Add1D("16", "log10 Energy (MeV) of f << 125 133 fHisto->Add1D("17", "log10 Energy (MeV) of b << 126 histo->add1D("7", 134 fHisto->Add1D("18", "log10 Energy (MeV) of l << 127 "Log10 Energy (MeV) of charged pions",nBinsE,-4.,6.,1.0); 135 fHisto->Add1D("19", "log10 Energy (MeV) of l << 128 136 fHisto->Add1D("20", "Log10 Energy (MeV) of p << 129 histo->add1D("8", 137 fHisto->Add1D("21", "Log10 Energy (MeV) of p << 130 "Log10 Energy (MeV) of pi0",nBinsE,-4.,6.,1.0); 138 fHisto->Add1D("22", "Energy deposition in th << 131 139 1.0); << 132 histo->add1D("9", 140 fHisto->Add1D("23", "EM energy deposition in << 133 "Log10 Energy (MeV) of charged kaons",nBinsE,-4.,6.,1.0); 141 1.0); << 134 142 fHisto->Add1D("24", "Pion energy deposition << 135 histo->add1D("10", 143 1.1, 1.0); << 136 "Log10 Energy (MeV) of neutral kaons",nBinsE,-4.,6.,1.0); 144 fHisto->Add1D("25", "Proton energy depositio << 137 145 1.1, 1.0); << 138 histo->add1D("11", 146 fHisto->Add1D("26", "Energy (MeV) of pi+", f << 139 "Log10 Energy (MeV) of deuterons and tritons",nBinsE,-5.,5.,1.0); 147 fHisto->Add1D("27", "Energy (MeV) of pi-", f << 140 148 fHisto->Add1D("28", "Energy (MeV) of pi0", f << 141 histo->add1D("12", >> 142 "Log10 Energy (MeV) of He3 and alpha",nBinsE,-5.,5.,1.0); >> 143 >> 144 histo->add1D("13", >> 145 "Log10 Energy (MeV) of Generic Ions",nBinsE,-5.,5.,1.0); >> 146 >> 147 histo->add1D("14", >> 148 "Log10 Energy (MeV) of muons",nBinsE,-4.,6.,1.0); >> 149 >> 150 histo->add1D("15", >> 151 "log10 Energy (MeV) of side-leaked neutrons",nBinsE,-5.,5.,1.0); >> 152 >> 153 histo->add1D("16", >> 154 "log10 Energy (MeV) of forward-leaked neutrons",nBinsE,-5.,5.,1.0); >> 155 >> 156 histo->add1D("17", >> 157 "log10 Energy (MeV) of backward-leaked neutrons",nBinsE,-5.,5.,1.0); >> 158 >> 159 histo->add1D("18", >> 160 "log10 Energy (MeV) of leaking protons",nBinsE,-4.,6.,1.0); >> 161 >> 162 histo->add1D("19", >> 163 "log10 Energy (MeV) of leaking charged pions",nBinsE,-4.,6.,1.0); >> 164 149 } 165 } 150 166 151 //....oooOO0OOooo........oooOO0OOooo........oo 167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 152 168 153 void HistoManager::BeginOfRun() 169 void HistoManager::BeginOfRun() 154 { 170 { 155 fR2 = fRadius * fRadius; << 171 absZ0 = -0.5*length; 156 fAbsZ0 = 0.5 * fLength; << 172 n_evt = 0; 157 fNevt = 0; << 173 n_elec = 0; 158 fNelec = 0; << 174 n_posit = 0; 159 fNposit = 0; << 175 n_gam = 0; 160 fNgam = 0; << 176 n_step = 0; 161 fNstep = 0; << 177 n_prot_leak = 0; 162 fNprot_leak = 0; << 178 n_pion_leak = 0; 163 fNpiofNleak = 0; << 179 n_ions = 0; 164 fNions = 0; << 180 n_deut = 0; 165 fNdeut = 0; << 181 n_alpha = 0; 166 fNalpha = 0; << 182 n_kaons = 0; 167 fNkaons = 0; << 183 n_muons = 0; 168 fNmuons = 0; << 184 n_cpions = 0; 169 fNcpions = 0; << 185 n_pi0 = 0; 170 fNpi0 = 0; << 186 n_neutron = 0; 171 fNneutron = 0; << 187 n_proton = 0; 172 fNproton = 0; << 188 n_aproton = 0; 173 fNaproton = 0; << 189 n_neu_forw = 0; 174 fNneu_forw = 0; << 190 n_neu_leak = 0; 175 fNneu_leak = 0; << 191 n_neu_back = 0; 176 fNneu_back = 0; << 192 177 << 193 bookHisto(); 178 fEdepSum = 0.0; << 194 histo->book(); 179 fEdepSum2 = 0.0; << 195 180 << 196 if(verbose > 0) 181 if (!fHistoBooked) { << 197 G4cout << "HistoManager: Histograms are booked and run has been started" 182 BookHisto(); << 198 <<G4endl<<" BeginOfRun (After histo->book)"<< G4endl; 183 } << 184 fHisto->Book(); << 185 << 186 if (fVerbose > 0) { << 187 G4cout << "HistoManager: Histograms are bo << 188 } << 189 } 199 } 190 200 191 //....oooOO0OOooo........oooOO0OOooo........oo 201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 192 202 193 void HistoManager::EndOfRun() 203 void HistoManager::EndOfRun() 194 { 204 { >> 205 195 G4cout << "HistoManager: End of run actions 206 G4cout << "HistoManager: End of run actions are started" << G4endl; 196 207 197 // Average values 208 // Average values 198 G4cout << "================================= << 209 G4cout<<"========================================================"<<G4endl; 199 210 200 G4double x = (G4double)fNevt; << 211 G4double x = (G4double)n_evt; 201 if (fNevt > 0) { << 212 if(n_evt > 0) x = 1.0/x; 202 x = 1.0 / x; << 203 } << 204 213 205 G4double xe = x * (G4double)fNelec; << 214 G4double xe = x*(G4double)n_elec; 206 G4double xg = x * (G4double)fNgam; << 215 G4double xg = x*(G4double)n_gam; 207 G4double xp = x * (G4double)fNposit; << 216 G4double xp = x*(G4double)n_posit; 208 G4double xs = x * (G4double)fNstep; << 217 G4double xs = x*(G4double)n_step; 209 G4double xn = x * (G4double)fNneutron; << 218 G4double xn = x*(G4double)n_neutron; 210 G4double xpn = x * (G4double)fNproton; << 219 G4double xpn = x*(G4double)n_proton; 211 G4double xap = x * (G4double)fNaproton; << 220 G4double xap = x*(G4double)n_aproton; 212 G4double xnf = x * (G4double)fNneu_forw; << 221 G4double xnf = x*(G4double)n_neu_forw; 213 G4double xnb = x * (G4double)fNneu_leak; << 222 G4double xnb = x*(G4double)n_neu_leak; 214 G4double xnbw = x * (G4double)fNneu_back; << 223 G4double xnbw= x*(G4double)n_neu_back; 215 G4double xpl = x * (G4double)fNprot_leak; << 224 G4double xpl = x*(G4double)n_prot_leak; 216 G4double xal = x * (G4double)fNpiofNleak; << 225 G4double xal = x*(G4double)n_pion_leak; 217 G4double xpc = x * (G4double)fNcpions; << 226 G4double xpc = x*(G4double)n_cpions; 218 G4double xp0 = x * (G4double)fNpi0; << 227 G4double xp0 = x*(G4double)n_pi0; 219 G4double xpk = x * (G4double)fNkaons; << 228 G4double xpk = x*(G4double)n_kaons; 220 G4double xpm = x * (G4double)fNmuons; << 229 G4double xpm = x*(G4double)n_muons; 221 G4double xid = x * (G4double)fNdeut; << 230 G4double xid = x*(G4double)n_deut; 222 G4double xia = x * (G4double)fNalpha; << 231 G4double xia = x*(G4double)n_alpha; 223 G4double xio = x * (G4double)fNions; << 232 G4double xio = x*(G4double)n_ions; 224 << 233 225 fEdepSum *= x; << 234 G4cout << "Beam particle " 226 fEdepSum2 *= x; << 235 << primaryDef->GetParticleName() <<G4endl; 227 fEdepSum2 -= fEdepSum * fEdepSum; << 236 G4cout << "Beam Energy(MeV) " 228 if (fEdepSum2 > 0.0) { << 237 << primaryKineticEnergy/MeV <<G4endl; 229 fEdepSum2 = std::sqrt(fEdepSum2); << 238 G4cout << "Number of events " << n_evt <<G4endl; 230 } << 231 else { << 232 fEdepSum2 = 0.0; << 233 } << 234 << 235 G4cout << "Beam particle << 236 G4cout << "Beam Energy(MeV) << 237 G4cout << "Number of events << 238 G4cout << std::setprecision(4) << "Average e << 239 << " RMS(MeV) " << fEdepSum2 / MeV << 240 G4cout << std::setprecision(4) << "Average n 239 G4cout << std::setprecision(4) << "Average number of steps " << xs << G4endl; 241 G4cout << std::setprecision(4) << "Average n 240 G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl; 242 G4cout << std::setprecision(4) << "Average n 241 G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl; 243 G4cout << std::setprecision(4) << "Average n 242 G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl; 244 G4cout << std::setprecision(4) << "Average n 243 G4cout << std::setprecision(4) << "Average number of neutrons " << xn << G4endl; 245 G4cout << std::setprecision(4) << "Average n 244 G4cout << std::setprecision(4) << "Average number of protons " << xpn << G4endl; 246 G4cout << std::setprecision(4) << "Average n 245 G4cout << std::setprecision(4) << "Average number of antiprotons " << xap << G4endl; 247 G4cout << std::setprecision(4) << "Average n 246 G4cout << std::setprecision(4) << "Average number of pi+ & pi- " << xpc << G4endl; 248 G4cout << std::setprecision(4) << "Average n 247 G4cout << std::setprecision(4) << "Average number of pi0 " << xp0 << G4endl; 249 G4cout << std::setprecision(4) << "Average n 248 G4cout << std::setprecision(4) << "Average number of kaons " << xpk << G4endl; 250 G4cout << std::setprecision(4) << "Average n 249 G4cout << std::setprecision(4) << "Average number of muons " << xpm << G4endl; 251 G4cout << std::setprecision(4) << "Average n 250 G4cout << std::setprecision(4) << "Average number of deuterons+tritons " << xid << G4endl; 252 G4cout << std::setprecision(4) << "Average n 251 G4cout << std::setprecision(4) << "Average number of He3+alpha " << xia << G4endl; 253 G4cout << std::setprecision(4) << "Average n 252 G4cout << std::setprecision(4) << "Average number of ions " << xio << G4endl; 254 G4cout << std::setprecision(4) << "Average n 253 G4cout << std::setprecision(4) << "Average number of forward neutrons " << xnf << G4endl; 255 G4cout << std::setprecision(4) << "Average n 254 G4cout << std::setprecision(4) << "Average number of reflected neutrons " << xnb << G4endl; 256 G4cout << std::setprecision(4) << "Average n 255 G4cout << std::setprecision(4) << "Average number of leaked neutrons " << xnbw << G4endl; 257 G4cout << std::setprecision(4) << "Average n 256 G4cout << std::setprecision(4) << "Average number of proton leak " << xpl << G4endl; 258 G4cout << std::setprecision(4) << "Average n 257 G4cout << std::setprecision(4) << "Average number of pion leak " << xal << G4endl; 259 G4cout << "================================= << 258 G4cout<<"========================================================"<<G4endl; 260 G4cout << G4endl; << 259 G4cout<<G4endl; 261 260 262 // normalise histograms 261 // normalise histograms 263 for (G4int i = 0; i < fNHisto; i++) { << 262 for(G4int i=0; i<nHisto; i++) { 264 fHisto->ScaleH1(i, x); << 263 histo->scale(i,x); 265 } 264 } 266 265 267 fHisto->Save(); << 266 if(verbose > 1) histo->print(0); 268 } << 267 histo->save(); 269 << 268 histo->reset(); 270 //....oooOO0OOooo........oooOO0OOooo........oo << 271 << 272 void HistoManager::BeginOfEvent() << 273 { << 274 fEdepEvt = 0.0; << 275 fEdepEM = 0.0; << 276 fEdepPI = 0.0; << 277 fEdepP = 0.0; << 278 } << 279 << 280 //....oooOO0OOooo........oooOO0OOooo........oo << 281 << 282 void HistoManager::EndOfEvent() << 283 { << 284 fEdepSum += fEdepEvt; << 285 fEdepSum2 += fEdepEvt * fEdepEvt; << 286 fHisto->Fill(21, fEdepEvt / fPrimaryKineticE << 287 fHisto->Fill(22, fEdepEM / fPrimaryKineticEn << 288 fHisto->Fill(23, fEdepPI / fPrimaryKineticEn << 289 fHisto->Fill(24, fEdepP / fPrimaryKineticEne << 290 } 269 } 291 270 292 //....oooOO0OOooo........oooOO0OOooo........oo 271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 293 272 294 void HistoManager::ScoreNewTrack(const G4Track 273 void HistoManager::ScoreNewTrack(const G4Track* track) 295 { 274 { 296 const G4ParticleDefinition* pd = track->GetD 275 const G4ParticleDefinition* pd = track->GetDefinition(); 297 const G4String name = pd->GetParticleName(); << 276 G4String name = pd->GetParticleName(); 298 const G4double ee = track->GetKineticEnergy( << 277 G4double e = track->GetKineticEnergy(); 299 G4double e = ee; << 300 278 301 // Primary track 279 // Primary track 302 if (0 == track->GetParentID()) { << 280 if(0 == track->GetParentID()) { 303 fNevt++; << 281 304 fPrimaryKineticEnergy = e; << 282 n_evt++; 305 fPrimaryDef = pd; << 283 primaryKineticEnergy = e; >> 284 primaryDef = pd; 306 G4ThreeVector dir = track->GetMomentumDire 285 G4ThreeVector dir = track->GetMomentumDirection(); 307 if (1 < fVerbose) << 286 if(1 < verbose) 308 G4cout << "### Primary " << name << " ki << 287 G4cout << "### Primary " << name 309 << "; m(MeV)= " << pd->GetPDGMass << 288 << " kinE(MeV)= " << e/MeV 310 << "; dir= " << track->GetMoment << 289 << "; m(MeV)= " << pd->GetPDGMass()/MeV >> 290 << "; pos(mm)= " << track->GetPosition()/mm >> 291 << "; dir= " << track->GetMomentumDirection() >> 292 << G4endl; 311 293 312 // Secondary track 294 // Secondary track 313 } << 295 } else { 314 else { << 296 if(1 < verbose) 315 if (1 < fVerbose) << 297 G4cout << "=== Secondary " << name 316 G4cout << "=== Secondary " << name << " << 298 << " kinE(MeV)= " << e/MeV 317 << "; m(MeV)= " << pd->GetPDGMass << 299 << "; m(MeV)= " << pd->GetPDGMass()/MeV 318 << "; dir= " << track->GetMoment << 300 << "; pos(mm)= " << track->GetPosition()/mm 319 e = (e > 0.0) ? std::log10(e / MeV) : -100 << 301 << "; dir= " << track->GetMomentumDirection() 320 if (pd == G4Gamma::Gamma()) { << 302 << G4endl; 321 fNgam++; << 303 e = std::log10(e/MeV); 322 fHisto->Fill(1, e, 1.0); << 304 if(pd == G4Gamma::Gamma()) { 323 } << 305 n_gam++; 324 else if (pd == G4Electron::Electron()) { << 306 histo->fill(1,e,1.0); 325 fNelec++; << 307 } else if ( pd == G4Electron::Electron()) { 326 fHisto->Fill(2, e, 1.0); << 308 n_elec++; 327 } << 309 histo->fill(2,e,1.0); 328 else if (pd == G4Positron::Positron()) { << 310 } else if ( pd == G4Positron::Positron()) { 329 fNposit++; << 311 n_posit++; 330 fHisto->Fill(3, e, 1.0); << 312 histo->fill(3,e,1.0); 331 } << 313 } else if ( pd == G4Proton::Proton()) { 332 else if (pd == G4Proton::Proton()) { << 314 n_proton++; 333 fNproton++; << 315 histo->fill(4,e,1.0); 334 fHisto->Fill(4, e, 1.0); << 316 } else if ( pd == neutron) { 335 } << 317 n_neutron++; 336 else if (pd == fNeutron) { << 318 histo->fill(5,e,1.0); 337 fNneutron++; << 319 } else if ( pd == G4AntiProton::AntiProton()) { 338 fHisto->Fill(5, e, 1.0); << 320 n_aproton++; 339 } << 321 } else if ( pd == G4PionPlus::PionPlus() || pd == G4PionMinus::PionMinus()) { 340 else if (pd == G4AntiProton::AntiProton()) << 322 n_cpions++; 341 fNaproton++; << 323 histo->fill(6,e,1.0); 342 } << 324 } else if ( pd == G4PionZero::PionZero()) { 343 else if (pd == G4PionPlus::PionPlus()) { << 325 n_pi0++; 344 fNcpions++; << 326 histo->fill(7,e,1.0); 345 fHisto->Fill(6, e, 1.0); << 327 } else if ( pd == G4KaonPlus::KaonPlus() || pd == G4KaonMinus::KaonMinus()) { 346 fHisto->Fill(19, e, 1.0); << 328 n_kaons++; 347 fHisto->Fill(25, ee, 1.0); << 329 histo->fill(8,e,1.0); 348 } << 330 } else if ( pd == G4KaonZeroShort::KaonZeroShort() || pd == G4KaonZeroLong::KaonZeroLong()) { 349 else if (pd == G4PionMinus::PionMinus()) { << 331 n_kaons++; 350 fNcpions++; << 332 histo->fill(9,e,1.0); 351 fHisto->Fill(6, e, 1.0); << 333 } else if ( pd == G4Deuteron::Deuteron() || pd == G4Triton::Triton()) { 352 fHisto->Fill(20, e, 1.0); << 334 n_deut++; 353 fHisto->Fill(26, ee, 1.0); << 335 histo->fill(10,e,1.0); 354 } << 336 } else if ( pd == G4He3::He3() || pd == G4Alpha::Alpha()) { 355 else if (pd == G4PionZero::PionZero()) { << 337 n_alpha++; 356 fNpi0++; << 338 histo->fill(11,e,1.0); 357 fHisto->Fill(7, e, 1.0); << 339 } else if ( pd->GetParticleType() == "nucleus") { 358 fHisto->Fill(27, ee, 1.0); << 340 n_ions++; 359 } << 341 histo->fill(12,e,1.0); 360 else if (pd == G4KaonPlus::KaonPlus() || p << 342 } else if ( pd == G4MuonPlus::MuonPlus() || pd == G4MuonMinus::MuonMinus()) { 361 fNkaons++; << 343 n_muons++; 362 fHisto->Fill(8, e, 1.0); << 344 histo->fill(13,e,1.0); 363 } << 364 else if (pd == G4KaonZeroShort::KaonZeroSh << 365 fNkaons++; << 366 fHisto->Fill(9, e, 1.0); << 367 } << 368 else if (pd == G4Deuteron::Deuteron() || p << 369 fNdeut++; << 370 fHisto->Fill(10, e, 1.0); << 371 } << 372 else if (pd == G4He3::He3() || pd == G4Alp << 373 fNalpha++; << 374 fHisto->Fill(11, e, 1.0); << 375 } << 376 else if (pd->GetParticleType() == "nucleus << 377 fNions++; << 378 fHisto->Fill(12, e, 1.0); << 379 } << 380 else if (pd == G4MuonPlus::MuonPlus() || p << 381 fNmuons++; << 382 fHisto->Fill(13, e, 1.0); << 383 } 345 } 384 } 346 } 385 } 347 } 386 348 387 //....oooOO0OOooo........oooOO0OOooo........oo 349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 388 350 389 void HistoManager::AddTargetStep(const G4Step* 351 void HistoManager::AddTargetStep(const G4Step* step) 390 { 352 { 391 ++fNstep; << 353 n_step++; 392 const G4double fEdep = step->GetTotalEnergyD << 393 const G4Track* track = step->GetTrack(); 354 const G4Track* track = step->GetTrack(); 394 const G4ParticleDefinition* pd = track->GetD << 355 currentDef = track->GetDefinition(); 395 if (1 < fVerbose) { << 356 currentKinEnergy = track->GetKineticEnergy(); 396 G4cout << "TargetSD::ProcessHits: " << pd- << 397 << " E(MeV)=" << track->GetKineticE << 398 << " beta1= " << step->GetPreStepPo << 399 << " beta2= " << step->GetPostStepP << 400 << " weight= " << step->GetTrack()- << 401 << " t(ns)=" << track->GetGlobalTim << 402 } << 403 if (fEdep > 0.0) { << 404 G4ThreeVector pos = << 405 (step->GetPreStepPoint()->GetPosition() << 406 << 407 G4double z = pos.z() + fAbsZ0; << 408 << 409 // scoring << 410 fEdepEvt += fEdep; << 411 fHisto->Fill(0, z, fEdep); << 412 << 413 if (pd == G4Gamma::Gamma() || pd == G4Elec << 414 fEdepEM += fEdep; << 415 } << 416 else if (pd == G4PionPlus::PionPlus() || p << 417 fEdepPI += fEdep; << 418 } << 419 else if (pd == G4Proton::Proton() || pd == << 420 fEdepP += fEdep; << 421 } << 422 357 423 if (1 < fVerbose) { << 358 G4ThreeVector pos = 424 G4cout << "HistoManager::AddEnergy: e(ke << 359 (step->GetPreStepPoint()->GetPosition() + 425 << "; step(mm)= " << step->GetSte << 360 step->GetPostStepPoint()->GetPosition())*0.5; 426 << " E(MeV)= " << track->GetKinet << 361 427 } << 362 G4double z = pos.z() - absZ0; 428 } << 363 histo->fill(0,z,step->GetTotalEnergyDeposit()); >> 364 >> 365 if(1 < verbose) >> 366 G4cout << "HistoManager::AddEnergy: e(keV)= " << step->GetTotalEnergyDeposit()/keV >> 367 << "; z(mm)= " << z/mm >> 368 << "; step(mm)= " << step->GetStepLength()/mm >> 369 << " by " << currentDef->GetParticleName() >> 370 << " E(MeV)= " << currentKinEnergy/MeV >> 371 << G4endl; 429 } 372 } 430 373 431 //....oooOO0OOooo........oooOO0OOooo........oo 374 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 432 375 433 void HistoManager::AddLeakingParticle(const G4 376 void HistoManager::AddLeakingParticle(const G4Track* track) 434 { 377 { 435 const G4ParticleDefinition* pd = track->GetD << 378 const G4ParticleDefinition* pd = track->GetDefinition(); 436 const G4StepPoint* sp = track->GetStep()->Ge << 379 G4double e = std::log10(track->GetKineticEnergy()/MeV); 437 const G4double ekin = sp->GetKineticEnergy() << 438 G4double e = -100.; << 439 if (ekin > 0.0) { << 440 e = std::log10(ekin / CLHEP::MeV); << 441 } << 442 else { << 443 G4cout << "### HistoManager::AddLeakingPar << 444 << " Eprestep(MeV)=" << ekin << " << 445 return; << 446 } << 447 380 448 const G4ThreeVector& pos = sp->GetPosition() << 381 G4ThreeVector pos = track->GetPosition(); 449 const G4ThreeVector& dir = sp->GetMomentumDi << 382 G4ThreeVector dir = track->GetMomentumDirection(); 450 G4double x = pos.x(); 383 G4double x = pos.x(); 451 G4double y = pos.y(); 384 G4double y = pos.y(); 452 G4double z = pos.z(); 385 G4double z = pos.z(); 453 G4double vx = dir.x(); << 386 454 G4double vy = dir.y(); << 455 G4double vz = dir.z(); << 456 << 457 G4bool isLeaking = false; 387 G4bool isLeaking = false; 458 388 459 // Forward << 389 // Forward 460 const G4double del = 0.001 * CLHEP::mm; << 390 if(z > -absZ0 && dir.z() > 0.0) { 461 if (std::abs(z - fAbsZ0) < del && vz > 0.0) << 391 if(pd == neutron) { 462 isLeaking = true; << 392 n_neu_forw++; 463 if (pd == fNeutron) { << 393 histo->fill(15,e,1.0); 464 ++fNneu_forw; << 394 } else isLeaking = true; 465 fHisto->Fill(15, e, 1.0); << 466 } << 467 else << 468 isLeaking = true; << 469 395 470 // Backward 396 // Backward 471 } << 397 } else if (z < absZ0 && dir.z() < 0.0) { 472 else if (std::abs(z + fAbsZ0) < del && vz < << 398 if(pd == neutron) { 473 isLeaking = true; << 399 n_neu_leak++; 474 if (pd == fNeutron) { << 400 histo->fill(16,e,1.0); 475 ++fNneu_back; << 401 } else isLeaking = true; 476 fHisto->Fill(16, e, 1.0); << 477 } << 478 else << 479 isLeaking = true; << 480 402 481 // Side 403 // Side 482 } << 404 } else if (std::abs(z) <= -absZ0 && x*dir.x() + y*dir.y() > 0.0) { 483 else if (std::abs(z) <= fAbsZ0 + del && x * << 484 && std::abs(x * x + y * y - fR2) < << 485 { << 486 isLeaking = true; 405 isLeaking = true; 487 if (pd == fNeutron) { << 406 if(pd == neutron) { 488 ++fNneu_leak; << 407 n_neu_back++; 489 fHisto->Fill(14, e, 1.0); << 408 histo->fill(14,e,1.0); 490 } << 409 } else isLeaking = true; 491 else << 492 isLeaking = true; << 493 } 410 } 494 411 495 // protons and pions 412 // protons and pions 496 if (isLeaking) { << 413 if(isLeaking) { 497 if (pd == G4Proton::Proton()) { << 414 if(pd == G4Proton::Proton()) { 498 fHisto->Fill(17, e, 1.0); << 415 histo->fill(17,e,1.0); 499 ++fNprot_leak; << 416 n_prot_leak++; 500 } << 417 } else if (pd == G4PionPlus::PionPlus() || pd == G4PionMinus::PionMinus()) { 501 else if (pd == G4PionPlus::PionPlus() || p << 418 histo->fill(18,e,1.0); 502 fHisto->Fill(18, e, 1.0); << 419 n_pion_leak++; 503 ++fNpiofNleak; << 504 } 420 } 505 } 421 } 506 } 422 } 507 423 508 //....oooOO0OOooo........oooOO0OOooo........oo 424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 509 425 510 void HistoManager::SetVerbose(G4int val) << 426 void HistoManager::SetVerbose(G4int val) 511 { 427 { 512 fVerbose = val; << 428 verbose = val; 513 fHisto->SetVerbose(val); << 429 histo->setVerbose(val); 514 } 430 } 515 431 516 //....oooOO0OOooo........oooOO0OOooo........oo 432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 517 433 518 void HistoManager::Fill(G4int id, G4double x, << 519 { << 520 fHisto->Fill(id, x, w); << 521 } << 522 << 523 //....oooOO0OOooo........oooOO0OOooo........oo << 524 434