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