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