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