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