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/Hadr02/src/HistoManager.cc 26 /// \file hadronic/Hadr02/src/HistoManager.cc 27 /// \brief Implementation of the HistoManager 27 /// \brief Implementation of the HistoManager class 28 // 28 // 29 // 29 // 30 //-------------------------------------------- 30 //--------------------------------------------------------------------------- 31 // 31 // 32 // ClassName: HistoManager 32 // ClassName: HistoManager 33 // 33 // 34 // 34 // 35 // Author: V.Ivanchenko 30/01/01 35 // Author: V.Ivanchenko 30/01/01 36 // 36 // 37 // Modified: 37 // Modified: 38 // 04.06.2006 Adoptation of hadr01 (V.Ivanchen 38 // 04.06.2006 Adoptation of hadr01 (V.Ivanchenko) 39 // 16.11.2006 Add beamFlag (V.Ivanchenko) 39 // 16.11.2006 Add beamFlag (V.Ivanchenko) 40 // 16.10.2012 Renamed G4IonFTFPBinaryCascadePh 40 // 16.10.2012 Renamed G4IonFTFPBinaryCascadePhysics as G4IonPhysics (A.Ribon) 41 // 41 // 42 //-------------------------------------------- 42 //---------------------------------------------------------------------------- 43 // 43 // 44 44 45 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 46 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 47 47 48 #include "HistoManager.hh" 48 #include "HistoManager.hh" 49 << 49 #include "G4UnitsTable.hh" 50 #include "DetectorConstruction.hh" << 50 #include "G4Neutron.hh" 51 #include "Histo.hh" << 51 #include "G4Proton.hh" 52 #include "IonHIJINGPhysics.hh" << 53 #include "IonUrQMDPhysics.hh" << 54 << 55 #include "G4Alpha.hh" << 56 #include "G4AntiProton.hh" 52 #include "G4AntiProton.hh" 57 #include "G4BuilderType.hh" << 58 #include "G4Deuteron.hh" << 59 #include "G4Electron.hh" << 60 #include "G4Gamma.hh" 53 #include "G4Gamma.hh" 61 #include "G4He3.hh" << 54 #include "G4Electron.hh" 62 #include "G4KaonMinus.hh" << 55 #include "G4Positron.hh" 63 #include "G4KaonPlus.hh" << 64 #include "G4KaonZeroLong.hh" << 65 #include "G4KaonZeroShort.hh" << 66 #include "G4Material.hh" << 67 #include "G4MuonMinus.hh" << 68 #include "G4MuonPlus.hh" 56 #include "G4MuonPlus.hh" 69 #include "G4Neutron.hh" << 57 #include "G4MuonMinus.hh" 70 #include "G4NistManager.hh" << 71 #include "G4PionMinus.hh" << 72 #include "G4PionPlus.hh" 58 #include "G4PionPlus.hh" >> 59 #include "G4PionMinus.hh" 73 #include "G4PionZero.hh" 60 #include "G4PionZero.hh" 74 #include "G4Positron.hh" << 61 #include "G4KaonPlus.hh" 75 #include "G4Proton.hh" << 62 #include "G4KaonMinus.hh" 76 #include "G4RunManager.hh" << 63 #include "G4KaonZeroShort.hh" 77 #include "G4SystemOfUnits.hh" << 64 #include "G4KaonZeroLong.hh" >> 65 #include "G4Deuteron.hh" 78 #include "G4Triton.hh" 66 #include "G4Triton.hh" 79 #include "G4UnitsTable.hh" << 67 #include "G4He3.hh" >> 68 #include "G4Alpha.hh" >> 69 #include "Histo.hh" >> 70 #include "globals.hh" 80 #include "G4VModularPhysicsList.hh" 71 #include "G4VModularPhysicsList.hh" 81 #include "G4VPhysicsConstructor.hh" 72 #include "G4VPhysicsConstructor.hh" 82 #include "globals.hh" << 73 #include "G4RunManager.hh" >> 74 #include "DetectorConstruction.hh" >> 75 #include "G4NistManager.hh" >> 76 #include "G4Material.hh" >> 77 #include "G4BuilderType.hh" >> 78 #include "G4SystemOfUnits.hh" >> 79 >> 80 #include "IonUrQMDPhysics.hh" >> 81 #include "IonHIJINGPhysics.hh" 83 82 84 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 85 84 86 HistoManager* HistoManager::fManager = 0; 85 HistoManager* HistoManager::fManager = 0; 87 86 88 //....oooOO0OOooo........oooOO0OOooo........oo 87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 89 88 90 HistoManager* HistoManager::GetPointer() 89 HistoManager* HistoManager::GetPointer() 91 { 90 { 92 if (!fManager) { << 91 if(!fManager) { 93 static HistoManager manager; 92 static HistoManager manager; 94 fManager = &manager; 93 fManager = &manager; 95 } 94 } 96 return fManager; 95 return fManager; 97 } 96 } 98 97 99 //....oooOO0OOooo........oooOO0OOooo........oo 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 100 99 101 HistoManager::HistoManager() 100 HistoManager::HistoManager() 102 { 101 { 103 fVerbose = 0; << 102 fVerbose= 0; 104 fNSlices = 100; << 103 fNSlices = 100; 105 fNBinsE = 100; << 104 fNBinsE = 100; 106 fNHisto = 20; << 105 fNHisto = 20; 107 fLength = 300. * mm; << 106 fLength = 300.*mm; 108 fEdepMax = 1.0 * GeV; << 107 fEdepMax = 1.0*GeV; 109 108 110 fPrimaryDef = 0; << 109 fPrimaryDef= 0; 111 fPrimaryKineticEnergy = 0.0; 110 fPrimaryKineticEnergy = 0.0; 112 fMaterial = 0; << 111 fMaterial = 0; 113 fBeamFlag = true; << 112 fBeamFlag = true; 114 113 115 fHisto = new Histo(); << 114 fHisto = new Histo(); 116 fHisto->SetVerbose(fVerbose); 115 fHisto->SetVerbose(fVerbose); 117 fNeutron = G4Neutron::Neutron(); << 116 fNeutron = G4Neutron::Neutron(); 118 fPhysList = 0; << 117 fPhysList = 0; 119 fIonPhysics = 0; << 118 fIonPhysics= 0; 120 Initialise(); 119 Initialise(); 121 } 120 } 122 121 123 //....oooOO0OOooo........oooOO0OOooo........oo 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 124 123 125 void HistoManager::Initialise() 124 void HistoManager::Initialise() 126 { 125 { 127 fAbsZ0 = -0.5 * fLength; << 126 fAbsZ0 = -0.5*fLength; 128 fNevt = 0; << 127 fNevt = 0; 129 fNelec = 0; << 128 fNelec = 0; 130 fNposit = 0; << 129 fNposit = 0; 131 fNgam = 0; << 130 fNgam = 0; 132 fNstep = 0; << 131 fNstep = 0; 133 fNions = 0; << 132 fNions = 0; 134 fNdeut = 0; << 133 fNdeut = 0; 135 fNalpha = 0; << 134 fNalpha = 0; 136 fNkaons = 0; << 135 fNkaons = 0; 137 fNmuons = 0; << 136 fNmuons = 0; 138 fNcpions = 0; << 137 fNcpions = 0; 139 fNpi0 = 0; << 138 fNpi0 = 0; 140 fNneutron = 0; << 139 fNneutron = 0; 141 fNproton = 0; << 140 fNproton = 0; 142 fNaproton = 0; << 141 fNaproton = 0; 143 << 142 144 fEdepEvt = 0.0; << 143 fEdepEvt = 0.0; 145 fEdepSum = 0.0; << 144 fEdepSum = 0.0; 146 fEdepSum2 = 0.0; << 145 fEdepSum2 = 0.0; 147 } 146 } 148 147 149 //....oooOO0OOooo........oooOO0OOooo........oo 148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 150 149 151 HistoManager::~HistoManager() 150 HistoManager::~HistoManager() 152 { 151 { 153 delete fHisto; 152 delete fHisto; 154 } 153 } 155 154 156 //....oooOO0OOooo........oooOO0OOooo........oo 155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 157 156 158 void HistoManager::BookHisto() 157 void HistoManager::BookHisto() 159 { 158 { 160 fHisto->Add1D("0", "Energy deposition (MeV/m << 159 fHisto->Add1D("0","Energy deposition (MeV/mm/event) in the target", 161 MeV / mm); << 160 fNSlices,0.0,fLength/mm,MeV/mm); 162 fHisto->Add1D("1", "Log10 Energy (GeV) of ga << 161 fHisto->Add1D("1","Log10 Energy (GeV) of gammas",fNBinsE,-5.,5.,1.0); 163 fHisto->Add1D("2", "Log10 Energy (GeV) of el << 162 fHisto->Add1D("2","Log10 Energy (GeV) of electrons",fNBinsE,-5.,5.,1.0); 164 fHisto->Add1D("3", "Log10 Energy (GeV) of po << 163 fHisto->Add1D("3","Log10 Energy (GeV) of positrons",fNBinsE,-5.,5.,1.0); 165 fHisto->Add1D("4", "Log10 Energy (GeV) of pr << 164 fHisto->Add1D("4","Log10 Energy (GeV) of protons",fNBinsE,-5.,5.,1.0); 166 fHisto->Add1D("5", "Log10 Energy (GeV) of ne << 165 fHisto->Add1D("5","Log10 Energy (GeV) of neutrons",fNBinsE,-5.,5.,1.0); 167 fHisto->Add1D("6", "Log10 Energy (GeV) of ch << 166 fHisto->Add1D("6","Log10 Energy (GeV) of charged pions",fNBinsE,-4.,6.,1.0); 168 fHisto->Add1D("7", "Log10 Energy (GeV) of pi << 167 fHisto->Add1D("7","Log10 Energy (GeV) of pi0",fNBinsE,-4.,6.,1.0); 169 fHisto->Add1D("8", "Log10 Energy (GeV) of ch << 168 fHisto->Add1D("8","Log10 Energy (GeV) of charged kaons",fNBinsE,-4.,6.,1.0); 170 fHisto->Add1D("9", "Log10 Energy (GeV) of ne << 169 fHisto->Add1D("9","Log10 Energy (GeV) of neutral kaons",fNBinsE,-4.,6.,1.0); 171 fHisto->Add1D("10", "Log10 Energy (GeV) of d << 170 fHisto->Add1D("10","Log10 Energy (GeV) of deuterons and tritons", 172 fHisto->Add1D("11", "Log10 Energy (GeV) of H << 171 fNBinsE,-5.,5.,1.0); 173 fHisto->Add1D("12", "Log10 Energy (GeV) of G << 172 fHisto->Add1D("11","Log10 Energy (GeV) of He3 and alpha",fNBinsE,-5.,5.,1.0); 174 fHisto->Add1D("13", "Log10 Energy (GeV) of m << 173 fHisto->Add1D("12","Log10 Energy (GeV) of Generic Ions",fNBinsE,-5.,5.,1.0); 175 fHisto->Add1D("14", "Log10 Energy (GeV) of p << 174 fHisto->Add1D("13","Log10 Energy (GeV) of muons",fNBinsE,-4.,6.,1.0); 176 fHisto->Add1D("15", "Log10 Energy (GeV) of p << 175 fHisto->Add1D("14","Log10 Energy (GeV) of pi+",fNBinsE,-4.,6.,1.0); 177 fHisto->Add1D("16", "X Section (mb) of Secon << 176 fHisto->Add1D("15","Log10 Energy (GeV) of pi-",fNBinsE,-4.,6.,1.0); 178 1.0); << 177 fHisto->Add1D("16","X Section (mb) of Secondary Fragments Z with E>1 GeV (mb)" 179 fHisto->Add1D("17", "Secondary Fragment A E> << 178 ,25,0.5,25.5,1.0); 180 fHisto->Add1D("18", "Secondary Fragment Z E< << 179 fHisto->Add1D("17","Secondary Fragment A E>1 GeV",50,0.5,50.5,1.0); 181 fHisto->Add1D("19", "Secondary Fragment A E< << 180 fHisto->Add1D("18","Secondary Fragment Z E<1 GeV",25,0.5,25.5,1.0); 182 fHisto->Add1D("20", "X Section (mb) of Secon << 181 fHisto->Add1D("19","Secondary Fragment A E<1 GeV",50,0.5,50.5,1.0); 183 fHisto->Add1D("21", "Secondary Fragment A ", << 182 fHisto->Add1D("20","X Section (mb) of Secondary Fragments Z (mb) ", >> 183 25,0.5,25.5,1.0); >> 184 fHisto->Add1D("21","Secondary Fragment A ",50,0.5,50.5,1.0); 184 } 185 } 185 186 186 //....oooOO0OOooo........oooOO0OOooo........oo 187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 187 188 188 void HistoManager::BeginOfRun() 189 void HistoManager::BeginOfRun() 189 { 190 { 190 Initialise(); 191 Initialise(); 191 BookHisto(); 192 BookHisto(); 192 fHisto->Book(); 193 fHisto->Book(); 193 194 194 if (fVerbose > 0) { << 195 if(fVerbose > 0) { 195 G4cout << "HistoManager: Histograms are bo << 196 G4cout << "HistoManager: Histograms are booked and run has been started" 196 << " BeginOfRun (After fHisto->boo << 197 <<G4endl<<" BeginOfRun (After fHisto->book)"<< G4endl; 197 } 198 } 198 } 199 } 199 200 200 //....oooOO0OOooo........oooOO0OOooo........oo 201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 201 202 202 void HistoManager::EndOfRun() 203 void HistoManager::EndOfRun() 203 { 204 { >> 205 204 G4cout << "HistoManager: End of run actions 206 G4cout << "HistoManager: End of run actions are started" << G4endl; 205 207 206 // Average values 208 // Average values 207 G4cout << "================================= << 209 G4cout<<"========================================================"<<G4endl; 208 210 209 G4double x = (G4double)fNevt; 211 G4double x = (G4double)fNevt; 210 if (fNevt > 0) { << 212 if(fNevt > 0) { x = 1.0/x; } 211 x = 1.0 / x; << 212 } << 213 213 214 G4double xe = x * fNelec; << 214 G4double xe = x*fNelec; 215 G4double xg = x * fNgam; << 215 G4double xg = x*fNgam; 216 G4double xp = x * fNposit; << 216 G4double xp = x*fNposit; 217 G4double xs = x * fNstep; << 217 G4double xs = x*fNstep; 218 G4double xn = x * fNneutron; << 218 G4double xn = x*fNneutron; 219 G4double xpn = x * fNproton; << 219 G4double xpn = x*fNproton; 220 G4double xap = x * fNaproton; << 220 G4double xap = x*fNaproton; 221 G4double xpc = x * fNcpions; << 221 G4double xpc = x*fNcpions; 222 G4double xp0 = x * fNpi0; << 222 G4double xp0 = x*fNpi0; 223 G4double xpk = x * fNkaons; << 223 G4double xpk = x*fNkaons; 224 G4double xpm = x * fNmuons; << 224 G4double xpm = x*fNmuons; 225 G4double xid = x * fNdeut; << 225 G4double xid = x*fNdeut; 226 G4double xia = x * fNalpha; << 226 G4double xia = x*fNalpha; 227 G4double xio = x * fNions; << 227 G4double xio = x*fNions; 228 228 229 fEdepSum *= x; << 229 fEdepSum *= x; 230 fEdepSum2 *= x; 230 fEdepSum2 *= x; 231 fEdepSum2 -= fEdepSum * fEdepSum; << 231 fEdepSum2 -= fEdepSum*fEdepSum; 232 if (fEdepSum2 > 0.0) << 232 if(fEdepSum2 > 0.0) fEdepSum2 = std::sqrt(fEdepSum2); 233 fEdepSum2 = std::sqrt(fEdepSum2); << 233 else fEdepSum2 = 0.0; 234 else << 234 235 fEdepSum2 = 0.0; << 235 G4cout << "Beam particle " 236 << 236 << fPrimaryDef->GetParticleName() <<G4endl; 237 G4cout << "Beam particle << 237 G4cout << "Beam Energy(GeV) " 238 G4cout << "Beam Energy(GeV) << 238 << fPrimaryKineticEnergy/GeV <<G4endl; 239 G4cout << "Number of events << 239 G4cout << "Number of events " 240 G4cout << std::setprecision(4) << "Average e << 240 << fNevt <<G4endl; 241 << " RMS(GeV) " << fEdepSum2 / GeV << 241 G4cout << std::setprecision(4) << "Average energy deposit (GeV) " 242 G4cout << std::setprecision(4) << "Average n << 242 << fEdepSum/GeV 243 G4cout << std::setprecision(4) << "Average n << 243 << " RMS(GeV) " << fEdepSum2/GeV << G4endl; 244 G4cout << std::setprecision(4) << "Average n << 244 G4cout << std::setprecision(4) << "Average number of steps " 245 G4cout << std::setprecision(4) << "Average n << 245 << xs << G4endl; 246 G4cout << std::setprecision(4) << "Average n << 246 G4cout << std::setprecision(4) << "Average number of gamma " 247 G4cout << std::setprecision(4) << "Average n << 247 << xg << G4endl; 248 G4cout << std::setprecision(4) << "Average n << 248 G4cout << std::setprecision(4) << "Average number of e- " 249 G4cout << std::setprecision(4) << "Average n << 249 << xe << G4endl; 250 G4cout << std::setprecision(4) << "Average n << 250 G4cout << std::setprecision(4) << "Average number of e+ " 251 G4cout << std::setprecision(4) << "Average n << 251 << xp << G4endl; 252 G4cout << std::setprecision(4) << "Average n << 252 G4cout << std::setprecision(4) << "Average number of neutrons " 253 G4cout << std::setprecision(4) << "Average n << 253 << xn << G4endl; 254 G4cout << std::setprecision(4) << "Average n << 254 G4cout << std::setprecision(4) << "Average number of protons " 255 G4cout << std::setprecision(4) << "Average n << 255 << xpn << G4endl; 256 G4cout << "================================= << 256 G4cout << std::setprecision(4) << "Average number of antiprotons " 257 G4cout << G4endl; << 257 << xap << G4endl; >> 258 G4cout << std::setprecision(4) << "Average number of pi+ & pi- " >> 259 << xpc << G4endl; >> 260 G4cout << std::setprecision(4) << "Average number of pi0 " >> 261 << xp0 << G4endl; >> 262 G4cout << std::setprecision(4) << "Average number of kaons " >> 263 << xpk << G4endl; >> 264 G4cout << std::setprecision(4) << "Average number of muons " >> 265 << xpm << G4endl; >> 266 G4cout << std::setprecision(4) << "Average number of deuterons+tritons " >> 267 << xid << G4endl; >> 268 G4cout << std::setprecision(4) << "Average number of He3+alpha " >> 269 << xia << G4endl; >> 270 G4cout << std::setprecision(4) << "Average number of ions " >> 271 << xio << G4endl; >> 272 G4cout<<"========================================================"<<G4endl; >> 273 G4cout<<G4endl; 258 274 259 // normalise histograms 275 // normalise histograms 260 for (G4int i = 0; i < fNHisto; i++) { << 276 for(G4int i=0; i<fNHisto; i++) { 261 fHisto->ScaleH1(i, x); << 277 fHisto->ScaleH1(i,x); 262 } 278 } 263 // will work only for pure material - 1 elem 279 // will work only for pure material - 1 element 264 // G4cout << fMaterial << G4endl; << 280 // G4cout << fMaterial << G4endl; 265 G4double F = 1000 / (fLength * barn * fMater << 281 G4double F= 1000/(fLength*barn*fMaterial->GetTotNbOfAtomsPerVolume()); 266 if (F > 0.0) { << 282 if(F > 0.0) { 267 fHisto->ScaleH1(16, F); << 283 fHisto->ScaleH1(16, F); 268 fHisto->ScaleH1(20, F); << 284 fHisto->ScaleH1(20, F); 269 } 285 } 270 286 271 fHisto->Save(); 287 fHisto->Save(); 272 } 288 } 273 289 274 //....oooOO0OOooo........oooOO0OOooo........oo 290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 275 291 276 void HistoManager::BeginOfEvent() 292 void HistoManager::BeginOfEvent() 277 { 293 { 278 ++fNevt; 294 ++fNevt; 279 fEdepEvt = 0.0; 295 fEdepEvt = 0.0; 280 } 296 } 281 << 297 282 //....oooOO0OOooo........oooOO0OOooo........oo 298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 283 299 284 void HistoManager::EndOfEvent() 300 void HistoManager::EndOfEvent() 285 { 301 { 286 fEdepSum += fEdepEvt; << 302 fEdepSum += fEdepEvt; 287 fEdepSum2 += fEdepEvt * fEdepEvt; << 303 fEdepSum2 += fEdepEvt*fEdepEvt; 288 } 304 } 289 305 290 //....oooOO0OOooo........oooOO0OOooo........oo 306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 291 307 292 void HistoManager::ScoreNewTrack(const G4Track 308 void HistoManager::ScoreNewTrack(const G4Track* track) 293 { 309 { 294 const G4ParticleDefinition* pd = track->GetD 310 const G4ParticleDefinition* pd = track->GetDefinition(); 295 G4String name = pd->GetParticleName(); 311 G4String name = pd->GetParticleName(); 296 G4double e = track->GetKineticEnergy(); 312 G4double e = track->GetKineticEnergy(); 297 313 298 // Primary track 314 // Primary track 299 if (0 == track->GetParentID()) { << 315 if(0 == track->GetParentID()) { >> 316 300 fPrimaryKineticEnergy = e; 317 fPrimaryKineticEnergy = e; 301 fPrimaryDef = pd; 318 fPrimaryDef = pd; 302 G4ThreeVector dir = track->GetMomentumDire 319 G4ThreeVector dir = track->GetMomentumDirection(); 303 if (1 < fVerbose) << 320 if(1 < fVerbose) 304 G4cout << "### Primary " << name << " ki << 321 G4cout << "### Primary " << name 305 << "; m(GeV)= " << pd->GetPDGMass << 322 << " kinE(GeV)= " << e/GeV 306 << "; dir= " << track->GetMoment << 323 << "; m(GeV)= " << pd->GetPDGMass()/GeV 307 << 324 << "; pos(mm)= " << track->GetPosition()/mm >> 325 << "; dir= " << track->GetMomentumDirection() >> 326 << G4endl; >> 327 308 // Secondary track 328 // Secondary track 309 } << 329 } else { 310 else { << 330 if(1 < fVerbose) { 311 if (1 < fVerbose) { << 331 G4cout << "=== Secondary " << name 312 G4cout << "=== Secondary " << name << " << 332 << " kinE(GeV)= " << e/GeV 313 << "; m(GeV)= " << pd->GetPDGMass << 333 << "; m(GeV)= " << pd->GetPDGMass()/GeV 314 << "; dir= " << track->GetMoment << 334 << "; pos(mm)= " << track->GetPosition()/mm 315 } << 335 << "; dir= " << track->GetMomentumDirection() 316 e = std::log10(e / GeV); << 336 << G4endl; 317 if (pd == G4Gamma::Gamma()) { << 337 } >> 338 e = std::log10(e/GeV); >> 339 if(pd == G4Gamma::Gamma()) { 318 fNgam++; 340 fNgam++; 319 fHisto->Fill(1, e, 1.0); << 341 fHisto->Fill(1,e,1.0); 320 } << 342 } else if ( pd == G4Electron::Electron()) { 321 else if (pd == G4Electron::Electron()) { << 322 fNelec++; 343 fNelec++; 323 fHisto->Fill(2, e, 1.0); << 344 fHisto->Fill(2,e,1.0); 324 } << 345 } else if ( pd == G4Positron::Positron()) { 325 else if (pd == G4Positron::Positron()) { << 326 fNposit++; 346 fNposit++; 327 fHisto->Fill(3, e, 1.0); << 347 fHisto->Fill(3,e,1.0); 328 } << 348 } else if ( pd == G4Proton::Proton()) { 329 else if (pd == G4Proton::Proton()) { << 330 fNproton++; 349 fNproton++; 331 fHisto->Fill(4, e, 1.0); << 350 fHisto->Fill(4,e,1.0); 332 } << 351 } else if ( pd == fNeutron) { 333 else if (pd == fNeutron) { << 334 fNneutron++; 352 fNneutron++; 335 fHisto->Fill(5, e, 1.0); << 353 fHisto->Fill(5,e,1.0); 336 } << 354 } else if ( pd == G4AntiProton::AntiProton()) { 337 else if (pd == G4AntiProton::AntiProton()) << 338 fNaproton++; 355 fNaproton++; 339 } << 356 } else if ( pd == G4PionPlus::PionPlus() ) { 340 else if (pd == G4PionPlus::PionPlus()) { << 341 fNcpions++; 357 fNcpions++; 342 fHisto->Fill(6, e, 1.0); << 358 fHisto->Fill(6,e,1.0); 343 fHisto->Fill(14, e, 1.0); << 359 fHisto->Fill(14,e,1.0); 344 } << 360 345 else if (pd == G4PionMinus::PionMinus()) { << 361 } else if ( pd == G4PionMinus::PionMinus()) { 346 fNcpions++; 362 fNcpions++; 347 fHisto->Fill(6, e, 1.0); << 363 fHisto->Fill(6,e,1.0); 348 fHisto->Fill(15, e, 1.0); << 364 fHisto->Fill(15,e,1.0); 349 } << 365 350 else if (pd == G4PionZero::PionZero()) { << 366 } else if ( pd == G4PionZero::PionZero()) { 351 fNpi0++; 367 fNpi0++; 352 fHisto->Fill(7, e, 1.0); << 368 fHisto->Fill(7,e,1.0); 353 } << 369 } else if ( pd == G4KaonPlus::KaonPlus() || 354 else if (pd == G4KaonPlus::KaonPlus() || p << 370 pd == G4KaonMinus::KaonMinus()) { 355 fNkaons++; 371 fNkaons++; 356 fHisto->Fill(8, e, 1.0); << 372 fHisto->Fill(8,e,1.0); 357 } << 373 } else if ( pd == G4KaonZeroShort::KaonZeroShort() || 358 else if (pd == G4KaonZeroShort::KaonZeroSh << 374 pd == G4KaonZeroLong::KaonZeroLong()) { 359 fNkaons++; 375 fNkaons++; 360 fHisto->Fill(9, e, 1.0); << 376 fHisto->Fill(9,e,1.0); 361 } << 377 } else if ( pd == G4Deuteron::Deuteron() || pd == G4Triton::Triton()) { 362 else if (pd == G4Deuteron::Deuteron() || p << 363 fNdeut++; 378 fNdeut++; 364 fHisto->Fill(10, e, 1.0); << 379 fHisto->Fill(10,e,1.0); 365 } << 380 } else if ( pd == G4He3::He3() || pd == G4Alpha::Alpha()) { 366 else if (pd == G4He3::He3() || pd == G4Alp << 367 fNalpha++; 381 fNalpha++; 368 fHisto->Fill(11, e, 1.0); << 382 fHisto->Fill(11,e,1.0); 369 } << 383 } else if ( pd->GetParticleType() == "nucleus") { 370 else if (pd->GetParticleType() == "nucleus << 371 fNions++; 384 fNions++; 372 fHisto->Fill(12, e, 1.0); << 385 fHisto->Fill(12,e,1.0); 373 G4double Z = pd->GetPDGCharge() / eplus; << 386 G4double Z = pd->GetPDGCharge()/eplus; 374 G4double A = (G4double)pd->GetBaryonNumb 387 G4double A = (G4double)pd->GetBaryonNumber(); 375 if (e > 0.0) { << 388 if(e > 0.0) { 376 fHisto->Fill(16, Z, 1.0); << 389 fHisto->Fill(16,Z,1.0); 377 fHisto->Fill(17, A, 1.0); << 390 fHisto->Fill(17,A,1.0); >> 391 } else { >> 392 fHisto->Fill(18,Z,1.0); >> 393 fHisto->Fill(19,A,1.0); 378 } 394 } 379 else { << 395 fHisto->Fill(20,Z,1.0); 380 fHisto->Fill(18, Z, 1.0); << 396 fHisto->Fill(21,A,1.0); 381 fHisto->Fill(19, A, 1.0); << 397 } else if ( pd == G4MuonPlus::MuonPlus() || 382 } << 398 pd == G4MuonMinus::MuonMinus()) { 383 fHisto->Fill(20, Z, 1.0); << 384 fHisto->Fill(21, A, 1.0); << 385 } << 386 else if (pd == G4MuonPlus::MuonPlus() || p << 387 fNmuons++; 399 fNmuons++; 388 fHisto->Fill(13, e, 1.0); << 400 fHisto->Fill(13,e,1.0); 389 } 401 } 390 } 402 } 391 } 403 } 392 404 393 //....oooOO0OOooo........oooOO0OOooo........oo 405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 394 406 395 void HistoManager::AddTargetStep(const G4Step* 407 void HistoManager::AddTargetStep(const G4Step* step) 396 { 408 { 397 ++fNstep; 409 ++fNstep; 398 G4double edep = step->GetTotalEnergyDeposit( 410 G4double edep = step->GetTotalEnergyDeposit(); 399 411 400 if (edep > 0.0) { << 412 if(edep > 0.0) { 401 G4ThreeVector pos = << 413 G4ThreeVector pos = 402 (step->GetPreStepPoint()->GetPosition() << 414 (step->GetPreStepPoint()->GetPosition() + >> 415 step->GetPostStepPoint()->GetPosition())*0.5; 403 416 404 G4double z = pos.z() - fAbsZ0; 417 G4double z = pos.z() - fAbsZ0; 405 418 406 // scoring 419 // scoring 407 fEdepEvt += edep; 420 fEdepEvt += edep; 408 fHisto->Fill(0, z, edep); << 421 fHisto->Fill(0,z,edep); 409 << 422 410 if (1 < fVerbose) { << 423 if(1 < fVerbose) { 411 const G4Track* track = step->GetTrack(); 424 const G4Track* track = step->GetTrack(); 412 G4cout << "HistoManager::AddEnergy: e(ke << 425 G4cout << "HistoManager::AddEnergy: e(keV)= " << edep/keV 413 << "; step(mm)= " << step->GetSte << 426 << "; z(mm)= " << z/mm 414 << track->GetDefinition()->GetPar << 427 << "; step(mm)= " << step->GetStepLength()/mm 415 << " E(GeV)= " << track->GetKinet << 428 << " by " << track->GetDefinition()->GetParticleName() >> 429 << " E(GeV)= " << track->GetKineticEnergy()/GeV >> 430 << G4endl; 416 } 431 } 417 } 432 } 418 } 433 } 419 434 420 //....oooOO0OOooo........oooOO0OOooo........oo 435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 421 436 422 void HistoManager::SetVerbose(G4int val) << 437 void HistoManager::SetVerbose(G4int val) 423 { 438 { 424 fVerbose = val; << 439 fVerbose = val; 425 fHisto->SetVerbose(val); 440 fHisto->SetVerbose(val); 426 } 441 } 427 442 428 //....oooOO0OOooo........oooOO0OOooo........oo 443 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 429 444 430 void HistoManager::Fill(G4int id, G4double x, 445 void HistoManager::Fill(G4int id, G4double x, G4double w) 431 { 446 { 432 fHisto->Fill(id, x, w); 447 fHisto->Fill(id, x, w); 433 } 448 } 434 449 435 //....oooOO0OOooo........oooOO0OOooo........oo 450 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 436 451 437 void HistoManager::SetIonPhysics(const G4Strin 452 void HistoManager::SetIonPhysics(const G4String& nam) 438 { 453 { 439 if (fIonPhysics) { << 454 if(fIonPhysics) { 440 G4cout << "### HistoManager WARNING: Ion P << 455 G4cout << "### HistoManager WARNING: Ion Physics is already defined: <" 441 << "> is ignored!" << G4endl; << 456 << nam << "> is ignored!" << G4endl; 442 } << 457 } else if(nam == "HIJING") { 443 else if (nam == "HIJING") { << 444 #ifdef G4_USE_HIJING 458 #ifdef G4_USE_HIJING 445 fIonPhysics = new IonHIJINGPhysics(); 459 fIonPhysics = new IonHIJINGPhysics(); 446 fPhysList->ReplacePhysics(fIonPhysics); 460 fPhysList->ReplacePhysics(fIonPhysics); 447 G4RunManager::GetRunManager()->PhysicsHasB 461 G4RunManager::GetRunManager()->PhysicsHasBeenModified(); 448 G4cout << "### SetIonPhysics: Ion Physics << 462 G4cout << "### SetIonPhysics: Ion Physics FTFP/Binary is added" >> 463 << G4endl; 449 #else 464 #else 450 G4cout << "Error: Ion Physics HIJING is re << 465 G4cout << "Error: Ion Physics HIJING is requested but is not available" >> 466 <<G4endl; 451 #endif 467 #endif 452 } << 468 } else if(nam == "UrQMD") { 453 else if (nam == "UrQMD") { << 454 #ifdef G4_USE_URQMD 469 #ifdef G4_USE_URQMD 455 fIonPhysics = new IonUrQMDPhysics(1); 470 fIonPhysics = new IonUrQMDPhysics(1); 456 fPhysList->ReplacePhysics(fIonPhysics); 471 fPhysList->ReplacePhysics(fIonPhysics); 457 G4RunManager::GetRunManager()->PhysicsHasB 472 G4RunManager::GetRunManager()->PhysicsHasBeenModified(); 458 G4cout << "### SetIonPhysics: Ion Physics << 473 G4cout << "### SetIonPhysics: Ion Physics UrQMD is added" >> 474 << G4endl; 459 #else 475 #else 460 G4cout << "Error: Ion Physics UrQMD is req << 476 G4cout << "Error: Ion Physics UrQMD is requested but is not available" >> 477 <<G4endl; 461 #endif 478 #endif 462 } << 479 } else { 463 else { << 480 G4cout << "### HistoManager WARNING: Ion Physics <" 464 G4cout << "### HistoManager WARNING: Ion P << 481 << nam << "> is unknown!" << G4endl; 465 } 482 } 466 } 483 } 467 484 468 //....oooOO0OOooo........oooOO0OOooo........oo 485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 486 469 487