1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // $Id: G4hTestStoppingPower.cc,v 1.18 2006/06/29 19:44:46 gunter Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-03-cand-01 $ |
---|
29 | // |
---|
30 | // ------------------------------------------------------------------- |
---|
31 | // GEANT 4 class file --- Copyright CERN 1998 |
---|
32 | // CERN Geneva Switzerland |
---|
33 | // |
---|
34 | // |
---|
35 | // File name: G4hTestStoppingPower.cc |
---|
36 | // |
---|
37 | // Author: Vladimir Ivanchenko |
---|
38 | // |
---|
39 | // Creation date: 28 July 2000 |
---|
40 | // |
---|
41 | // Modifications: |
---|
42 | // |
---|
43 | // ------------------------------------------------------------------- |
---|
44 | |
---|
45 | #include "globals.hh" |
---|
46 | #include "G4ios.hh" |
---|
47 | #include "G4UnitsTable.hh" |
---|
48 | #include <fstream> |
---|
49 | #include <iomanip> |
---|
50 | |
---|
51 | #include "G4Material.hh" |
---|
52 | #include "G4Element.hh" |
---|
53 | #include "G4ProcessManager.hh" |
---|
54 | #include "G4hLowEnergyIonisation.hh" |
---|
55 | #include "G4hBetheBlochModel.hh" |
---|
56 | #include "G4hParametrisedLossModel.hh" |
---|
57 | #include "G4hNuclearStoppingModel.hh" |
---|
58 | #include "G4QAOLowEnergyLoss.hh" |
---|
59 | #include "G4VParticleChange.hh" |
---|
60 | #include "G4ParticleChange.hh" |
---|
61 | #include "G4DynamicParticle.hh" |
---|
62 | |
---|
63 | #include "G4hZiegler1977p.hh" |
---|
64 | #include "G4hZiegler1985p.hh" |
---|
65 | #include "G4hZiegler1977He.hh" |
---|
66 | #include "G4hICRU49p.hh" |
---|
67 | #include "G4hICRU49He.hh" |
---|
68 | #include "G4hSRIM2000p.hh" |
---|
69 | #include "G4hSRIM2003p.hh" |
---|
70 | |
---|
71 | #include "G4hZiegler1977Nuclear.hh" |
---|
72 | #include "G4hZiegler1985Nuclear.hh" |
---|
73 | //#include "G4hMollereNuclear.hh" // exist no more |
---|
74 | |
---|
75 | |
---|
76 | #include "G4Box.hh" |
---|
77 | #include "G4PVPlacement.hh" |
---|
78 | |
---|
79 | #include "G4Step.hh" |
---|
80 | #include "G4GRSVolume.hh" |
---|
81 | |
---|
82 | #include "G4ParticleTypes.hh" |
---|
83 | #include "G4ParticleTable.hh" |
---|
84 | #include "G4ParticleDefinition.hh" |
---|
85 | #include "G4ParticleWithCuts.hh" |
---|
86 | #include "G4Ions.hh" |
---|
87 | |
---|
88 | // New Histogramming (from AIDA and Anaphe): |
---|
89 | #include <memory> // for the auto_ptr(T> |
---|
90 | |
---|
91 | #include "AIDA/AIDA.h" |
---|
92 | |
---|
93 | |
---|
94 | #include "hTest/include/G4IonC12.hh" |
---|
95 | #include "hTest/include/G4IonAr40.hh" |
---|
96 | |
---|
97 | #include "G4Timer.hh" |
---|
98 | |
---|
99 | int main() |
---|
100 | { |
---|
101 | // ---- HBOOK initialization |
---|
102 | |
---|
103 | G4String hFile = "htest.hbook"; |
---|
104 | |
---|
105 | //--------- Materials definition --------- |
---|
106 | |
---|
107 | G4Material* Be = new G4Material("Beryllium", 4., 9.01*g/mole, 1.848*g/cm3); |
---|
108 | G4Material* Graphite = new G4Material("Graphite",6.,12.0*g/mole,2.265*g/cm3); |
---|
109 | G4Material* Al = new G4Material("Aluminium", 13., 26.98*g/mole, 2.7 *g/cm3); |
---|
110 | G4Material* Si = new G4Material("Silicon", 14., 28.055*g/mole, 2.33*g/cm3); |
---|
111 | G4Material* LAr = new G4Material("LArgon", 18., 39.95*g/mole, 1.393*g/cm3); |
---|
112 | G4Material* Ti = new G4Material("Titan", 22., 47.867*g/mole, 4.54*g/cm3); |
---|
113 | G4Material* Fe = new G4Material("Iron", 26., 55.85*g/mole, 7.87*g/cm3); |
---|
114 | G4Material* Cu = new G4Material("Copper", 29., 63.55*g/mole, 8.96*g/cm3); |
---|
115 | G4Material* W = new G4Material("Tungsten",74., 183.85*g/mole, 19.30*g/cm3); |
---|
116 | G4Material* Pb = new G4Material("Lead", 82., 207.19*g/mole, 11.35*g/cm3); |
---|
117 | G4Material* U = new G4Material("Uranium", 92., 238.03*g/mole, 18.95*g/cm3); |
---|
118 | |
---|
119 | G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole); |
---|
120 | G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole); |
---|
121 | G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole); |
---|
122 | G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole); |
---|
123 | G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole); |
---|
124 | |
---|
125 | // G4Material* water = new G4Material ("Water" ,"H_2O", 1.*g/cm3, 2); |
---|
126 | G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2); |
---|
127 | water->AddElement(H,2); |
---|
128 | water->AddElement(O,1); |
---|
129 | |
---|
130 | // G4Material* ethane = new G4Material ("Ethane" ,"C_2H_6", 0.4241*g/cm3, 2); |
---|
131 | G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2); |
---|
132 | ethane->AddElement(H,6); |
---|
133 | ethane->AddElement(C,2); |
---|
134 | |
---|
135 | // G4Material* csi = new G4Material ("CsI" , "CsI", 4.53*g/cm3, 2); |
---|
136 | G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2); |
---|
137 | csi->AddElement(Cs,1); |
---|
138 | csi->AddElement(I,1); |
---|
139 | |
---|
140 | G4Material* el[93]; |
---|
141 | G4double z; |
---|
142 | G4int j; |
---|
143 | |
---|
144 | for ( j=1; j<93; j++) |
---|
145 | { |
---|
146 | z = G4double(j); |
---|
147 | el[j] = new G4Material ("atom", z , 2.0*z*g/mole, 10.0 *g/cm3); |
---|
148 | } |
---|
149 | |
---|
150 | const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); |
---|
151 | |
---|
152 | // create table |
---|
153 | G4int numOfMaterials = G4Material::GetNumberOfMaterials(); |
---|
154 | |
---|
155 | G4double dimx = 1*mm, dimy = 1*mm, dimz = 1*mm; |
---|
156 | G4int imat = 0; |
---|
157 | G4cout << "Available materials are: " << G4endl; |
---|
158 | for (imat = 0; imat < numOfMaterials; imat++) { |
---|
159 | G4cout << imat << ") " << (*theMaterialTable)[imat]->GetName() << G4endl; |
---|
160 | } |
---|
161 | |
---|
162 | |
---|
163 | // Geometry definitions |
---|
164 | G4Box* theFrame = new G4Box ("Frame",dimx, dimy, dimz); |
---|
165 | |
---|
166 | G4LogicalVolume* LogicalFrame = new G4LogicalVolume(theFrame, |
---|
167 | (*theMaterialTable)[0], |
---|
168 | "LFrame", 0, 0, 0); |
---|
169 | |
---|
170 | G4PVPlacement* PhysicalFrame = new G4PVPlacement(0,G4ThreeVector(), |
---|
171 | "PFrame",LogicalFrame,0,false,0); |
---|
172 | |
---|
173 | G4cout << PhysicalFrame << G4endl; |
---|
174 | |
---|
175 | //--------- Particle definition --------- |
---|
176 | |
---|
177 | G4ParticleDefinition* gamma = G4Gamma::Gamma(); |
---|
178 | G4ParticleDefinition* electron = G4Electron::Electron(); |
---|
179 | G4ParticleDefinition* proton = G4Proton::Proton(); |
---|
180 | G4ParticleDefinition* antiproton =G4AntiProton::AntiProton(); |
---|
181 | G4ParticleDefinition* deuteron = G4Deuteron::Deuteron(); |
---|
182 | G4ParticleDefinition* alpha = G4Alpha::Alpha(); |
---|
183 | |
---|
184 | G4Ions* iC12 = new G4Ions::G4Ions( |
---|
185 | "IonC12", 11.14945*GeV, 0.0*MeV, +6.0*eplus, |
---|
186 | 0, +1, 0, |
---|
187 | 0, 0, 0, |
---|
188 | "static_nucleus", 0, +12, 0, |
---|
189 | true, -1.0, NULL); |
---|
190 | |
---|
191 | G4Ions* iAr40 = new G4Ions::G4Ions( |
---|
192 | "IonAr40", 37.291*GeV, 0.0*MeV, +18.0*eplus, |
---|
193 | 0, +1, 0, |
---|
194 | 0, 0, 0, |
---|
195 | "static_nucleus", 0, +40, 0, |
---|
196 | true, -1.0, NULL); |
---|
197 | |
---|
198 | G4Ions* iFe56 = new G4Ions::G4Ions( |
---|
199 | "IonFe56", 52.0308*GeV, 0.0*MeV, +26.0*eplus, |
---|
200 | 0, +1, 0, |
---|
201 | 0, 0, 0, |
---|
202 | "static_nucleus", 0, +56, 0, |
---|
203 | true, -1.0, NULL); |
---|
204 | |
---|
205 | G4ParticleDefinition* ionC12 = iC12->IonsDefinition(); |
---|
206 | G4ParticleDefinition* ionAr40 = iAr40->IonsDefinition(); |
---|
207 | G4ParticleDefinition* ionFe56 = iFe56->IonsDefinition(); |
---|
208 | |
---|
209 | G4ParticleDefinition* part[7]; |
---|
210 | part[0] = proton; |
---|
211 | part[1] = antiproton; |
---|
212 | part[2] = deuteron; |
---|
213 | part[3] = alpha; |
---|
214 | part[4] = ionC12; |
---|
215 | part[5] = ionAr40; |
---|
216 | part[6] = ionFe56; |
---|
217 | |
---|
218 | G4double ecut = 1000.0*mm; |
---|
219 | G4double pcut = 0.0001*mm; |
---|
220 | gamma->SetCuts(ecut); |
---|
221 | electron->SetCuts(ecut); |
---|
222 | G4cout << "Cuts are following: cutElectron = " << ecut |
---|
223 | << " mm; cutProton = " << pcut << " mm" << G4endl; |
---|
224 | |
---|
225 | //--------- Ionisation processes definition and build physics table -- |
---|
226 | |
---|
227 | // Define models for parametrisation of electronic energy losses |
---|
228 | // G4VLowEnergyModel* theBetheBlochModel = |
---|
229 | // new G4hBetheBlochModel("Bethe-Bloch") ; |
---|
230 | // theProtonModel = new G4hParametrisedLossModel(theProtonTable) ; |
---|
231 | // theAntiProtonModel = new G4QAOLowEnergyLoss(theAntiProtonTable) ; |
---|
232 | |
---|
233 | G4hNuclearStoppingModel* theNuclearStoppingModel = |
---|
234 | // new G4hNuclearStoppingModel("Ziegler1977") ; |
---|
235 | // new G4hNuclearStoppingModel("Ziegler1985") ; |
---|
236 | new G4hNuclearStoppingModel("ICRU_R49") ; |
---|
237 | |
---|
238 | G4VLowEnergyModel* theIonEffChargeModel = |
---|
239 | new G4hIonEffChargeSquare("Ziegler1988") ; |
---|
240 | |
---|
241 | |
---|
242 | G4hLowEnergyIonisation* hIon[8]; |
---|
243 | G4ProcessManager* theProcessManager[8]; |
---|
244 | G4int i; |
---|
245 | |
---|
246 | G4cout << "Define processes!" << G4endl; |
---|
247 | |
---|
248 | for( i=0; i<7; i++) { |
---|
249 | G4cout << "Ionisation process for particle " << i << " " << part[i]->GetParticleName() << G4endl; |
---|
250 | theProcessManager[i] = new G4ProcessManager(part[i]); |
---|
251 | part[i]->SetProcessManager(theProcessManager[i]); |
---|
252 | hIon[i] = new G4hLowEnergyIonisation(); |
---|
253 | hIon[i]->SetEnlossFluc(false) ; |
---|
254 | |
---|
255 | // hIon[i]->SetBarkasOff(); |
---|
256 | // hIon[i]->SetNuclearStoppingOff(); |
---|
257 | // hIon[i]->SetStoppingPowerTableName("ICRU_R49p"); |
---|
258 | |
---|
259 | theProcessManager[i]->AddProcess(hIon[i]); |
---|
260 | hIon[i]->BuildPhysicsTable(*part[i]); |
---|
261 | } |
---|
262 | |
---|
263 | |
---|
264 | |
---|
265 | //----------- Histogram ------------------------------------- |
---|
266 | |
---|
267 | G4cout << "Fill Hbook!" << G4endl; |
---|
268 | |
---|
269 | // Creating the analysis factory |
---|
270 | std::auto_ptr< AIDA::IAnalysisFactory > af( AIDA_createAnalysisFactory() ); |
---|
271 | |
---|
272 | // Creating the tree factory |
---|
273 | std::auto_ptr< AIDA::ITreeFactory > tf( af->createTreeFactory() ); |
---|
274 | |
---|
275 | // Creating a tree mapped to a new hbook file. |
---|
276 | std::auto_ptr< AIDA::ITree > tree( tf->create( hFile,"hbook", false,false) ); |
---|
277 | std::cout << "Tree store : " << tree->storeName() << std::endl; |
---|
278 | |
---|
279 | // Creating a histogram factory, whose histograms will be handled by the tree |
---|
280 | std::auto_ptr< AIDA::IHistogramFactory > hf( af->createHistogramFactory( *tree ) ); |
---|
281 | |
---|
282 | // G4Material* material ; |
---|
283 | |
---|
284 | G4double minE = 1.0*eV, maxE = 10000.0*MeV, s; |
---|
285 | const G4int num = 200; |
---|
286 | G4double tkin = 0.0; |
---|
287 | s = (std::log10(maxE)-std::log10(minE))/num; |
---|
288 | |
---|
289 | AIDA::IHistogram1D* h[71] ; |
---|
290 | |
---|
291 | // Test on Stopping Powers for all elements |
---|
292 | |
---|
293 | h[1] = hf->createHistogram1D("1","p 40 keV (keV*cm2/10^15!atoms) Ziegler1977p" |
---|
294 | ,92,0.5,92.5) ; |
---|
295 | h[2] = hf->createHistogram1D("2","p 100 keV (keV*cm2/10^15!atoms) Ziegler1977p" |
---|
296 | ,92,0.5,92.5) ; |
---|
297 | h[3] = hf->createHistogram1D("3","p 400 keV (keV*cm2/10^15!atoms) Ziegler1977p" |
---|
298 | ,92,0.5,92.5) ; |
---|
299 | |
---|
300 | h[4] = hf->createHistogram1D("4","p 1 MeV (keV*cm2/10^15!atoms) Ziegler1977p" |
---|
301 | ,92,0.5,92.5) ; |
---|
302 | h[5] = hf->createHistogram1D("5","p 4 MeV (keV*cm2/10^15!atoms) Ziegler1977p" |
---|
303 | ,92,0.5,92.5) ; |
---|
304 | |
---|
305 | h[6] = hf->createHistogram1D("6","p 40 keV (keV*cm2/10^15!atoms) ICRU_49p" |
---|
306 | ,92,0.5,92.5) ; |
---|
307 | h[7] = hf->createHistogram1D("7","p 100 keV (keV*cm2/10^15!atoms) ICRU_49p" |
---|
308 | ,92,0.5,92.5) ; |
---|
309 | h[8] = hf->createHistogram1D("8","p 400 keV (keV*cm2/10^15!atoms) ICRU_49p" |
---|
310 | ,92,0.5,92.5) ; |
---|
311 | h[9] = hf->createHistogram1D("9","p 1 MeV (keV*cm2/10^15!atoms) ICRU_49p" |
---|
312 | ,92,0.5,92.5) ; |
---|
313 | h[10] = hf->createHistogram1D("10","p 4 MeV (keV*cm2/10^15!atoms) ICRU_49p" |
---|
314 | ,92,0.5,92.5) ; |
---|
315 | |
---|
316 | h[11] = hf->createHistogram1D("11","p 40 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
317 | ,92,0.5,92.5) ; |
---|
318 | h[12] = hf->createHistogram1D("12","p 100 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
319 | ,92,0.5,92.5) ; |
---|
320 | h[13] = hf->createHistogram1D("13","p 400 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
321 | ,92,0.5,92.5) ; |
---|
322 | h[14] = hf->createHistogram1D("14","p 1 MeV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
323 | ,92,0.5,92.5) ; |
---|
324 | h[15] = hf->createHistogram1D("15","p 4 MeV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
325 | ,92,0.5,92.5) ; |
---|
326 | |
---|
327 | h[16] = hf->createHistogram1D("16","He 10 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
328 | ,92,0.5,92.5) ; |
---|
329 | h[17] = hf->createHistogram1D("17","He 40 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
330 | ,92,0.5,92.5) ; |
---|
331 | h[18] = hf->createHistogram1D("18","He 100 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
332 | ,92,0.5,92.5) ; |
---|
333 | |
---|
334 | h[19] = hf->createHistogram1D("19","He 400 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
335 | ,92,0.5,92.5) ; |
---|
336 | h[20] = hf->createHistogram1D("20","He 1 MeV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
337 | ,92,0.5,92.5) ; |
---|
338 | h[21] = hf->createHistogram1D("21","He 4 MeV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
339 | ,92,0.5,92.5) ; |
---|
340 | h[22] = hf->createHistogram1D("22","He 10 keV (keV*cm2/10^15!atoms) ICRU49He" |
---|
341 | ,92,0.5,92.5) ; |
---|
342 | h[23] = hf->createHistogram1D("23","He 40 keV (keV*cm2/10^15!atoms) ICRU49He" |
---|
343 | ,92,0.5,92.5) ; |
---|
344 | h[24] = hf->createHistogram1D("24","He 100 keV (keV*cm2/10^15!atoms) ICRU49He" |
---|
345 | ,92,0.5,92.5) ; |
---|
346 | h[25] = hf->createHistogram1D("25","He 400 keV (keV*cm2/10^15!atoms) ICRU49He" |
---|
347 | ,92,0.5,92.5) ; |
---|
348 | h[26] = hf->createHistogram1D("26","He 1 MeV (keV*cm2/10^15!atoms) ICRU49He" |
---|
349 | ,92,0.5,92.5) ; |
---|
350 | h[27] = hf->createHistogram1D("27","He 4 MeV (keV*cm2/10^15!atoms) ICRU49He" |
---|
351 | ,92,0.5,92.5) ; |
---|
352 | |
---|
353 | h[28]= hf->createHistogram1D("28","p in C (MeV/mm) ICRU49p" |
---|
354 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
355 | h[29]= hf->createHistogram1D("29","p in Al (MeV/mm) ICRU49p" |
---|
356 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
357 | h[30]= hf->createHistogram1D("30","p in Si (MeV/mm) ICRU49p" |
---|
358 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
359 | h[31]= hf->createHistogram1D("31","p in Cu (MeV/mm) ICRU49p" |
---|
360 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
361 | h[32]= hf->createHistogram1D("32","p in Fe (MeV/mm) ICRU49p" |
---|
362 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
363 | h[33]= hf->createHistogram1D("33","p in Pb (MeV/mm) ICRU49p" |
---|
364 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
365 | h[34]= hf->createHistogram1D("34","p in C2H6 (MeV/mm) ICRU49p" |
---|
366 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
367 | h[35]= hf->createHistogram1D("35","p in H2O (MeV/mm) ICRU49p" |
---|
368 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
369 | h[36]= hf->createHistogram1D("36","p in lAr (MeV/mm) ICRU49p" |
---|
370 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
371 | h[37]= hf->createHistogram1D("37","p in CsI (MeV/mm) ICRU49p" |
---|
372 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
373 | |
---|
374 | |
---|
375 | h[38]= hf->createHistogram1D("38","p in C (MeV/mm)Ziegler1985p" |
---|
376 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
377 | h[39]= hf->createHistogram1D("39","p in Al (MeV/mm)Ziegler1985p" |
---|
378 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
379 | h[40]= hf->createHistogram1D("40","p in Si (MeV/mm)Ziegler1985p" |
---|
380 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
381 | h[41]= hf->createHistogram1D("41","p in Cu (MeV/mm)Ziegler1985p" |
---|
382 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
383 | h[42]= hf->createHistogram1D("42","p in Fe (MeV/mm)Ziegler1985p" |
---|
384 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
385 | h[43]= hf->createHistogram1D("43","p in Pb (MeV/mm)Ziegler1985p" |
---|
386 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
387 | h[44]= hf->createHistogram1D("44","p in C2H6 (MeV/mm)Ziegler1985p" |
---|
388 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
389 | h[45]= hf->createHistogram1D("45","p in H2O (MeV/mm)Ziegler1985p" |
---|
390 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
391 | h[46]= hf->createHistogram1D("46","p in lAr (MeV/mm)Ziegler1985p" |
---|
392 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
393 | h[47]= hf->createHistogram1D("47","p in CsI (MeV/mm)Ziegler1985p" |
---|
394 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
395 | |
---|
396 | |
---|
397 | h[48]= hf->createHistogram1D("48","He effective charge for Cu" |
---|
398 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
399 | h[49]= hf->createHistogram1D("49","C12 effective charge in Cu" |
---|
400 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
401 | |
---|
402 | h[50]= hf->createHistogram1D("50","He in Al (MeV/(mg/cm2)) ICRU49p" |
---|
403 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
404 | h[51]= hf->createHistogram1D("51","C12 in Al (MeV/(mg/cm2)) ICRU49p" |
---|
405 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
406 | h[52]= hf->createHistogram1D("52","Ar40 in Al (MeV/(mg/cm2)) ICRU49p" |
---|
407 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
408 | |
---|
409 | // Table with the data |
---|
410 | h[53] = hf->createHistogram1D("53","Data p 40 keV (keV*cm2/10^15!atoms) Ziegler1977p" |
---|
411 | ,92,0.5,92.5) ; |
---|
412 | h[54] = hf->createHistogram1D("54","Data He 40 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
413 | ,92,0.5,92.5) ; |
---|
414 | |
---|
415 | h[55] = hf->createHistogram1D("55","p 40 keV (keV*cm2/10^15!atoms) Ziegler1985p" |
---|
416 | ,92,0.5,92.5) ; |
---|
417 | h[56] = hf->createHistogram1D("56","p 100 keV (keV*cm2/10^15!atoms) Ziegler1985p" |
---|
418 | ,92,0.5,92.5) ; |
---|
419 | h[57] = hf->createHistogram1D("57","p 400 keV (keV*cm2/10^15!atoms) Ziegler1985p" |
---|
420 | ,92,0.5,92.5) ; |
---|
421 | h[58] = hf->createHistogram1D("58","p 1 MeV (keV*cm2/10^15!atoms) Ziegler1985p" |
---|
422 | ,92,0.5,92.5) ; |
---|
423 | h[59] = hf->createHistogram1D("59","p 4 MeV (keV*cm2/10^15!atoms) Ziegler1985p" |
---|
424 | ,92,0.5,92.5) ; |
---|
425 | // Histo for Antiproton |
---|
426 | |
---|
427 | h[60]= hf->createHistogram1D("60","pbar in C (MeV/mm) QAOLoss" |
---|
428 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
429 | h[61]= hf->createHistogram1D("61","pbar in Al (MeV/mm) QAOLoss" |
---|
430 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
431 | h[62]= hf->createHistogram1D("62","pbar in Si (MeV/mm) QAOLoss" |
---|
432 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
433 | h[63]= hf->createHistogram1D("63","pbar in Cu (MeV/mm) QAOLoss" |
---|
434 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
435 | h[64]= hf->createHistogram1D("64","pbar in Fe (MeV/mm) QAOLoss" |
---|
436 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
437 | h[65]= hf->createHistogram1D("65","pbar in Pb (MeV/mm) QAOLoss" |
---|
438 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
439 | h[66]= hf->createHistogram1D("66","pbar in C2H6 (MeV/mm) QAOLoss" |
---|
440 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
441 | h[67]= hf->createHistogram1D("67","pbar in H2O (MeV/mm) QAOLoss" |
---|
442 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
443 | h[68]= hf->createHistogram1D("68","pbar in lAr (MeV/mm) QAOLoss" |
---|
444 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
445 | h[69]= hf->createHistogram1D("69","pbar in CsI (MeV/mm) QAOLoss" |
---|
446 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
447 | |
---|
448 | h[70] = hf->createHistogram1D("70","p 6.5 MeV (keV*cm2/10^15!atoms) Ziegler77p" |
---|
449 | ,92,0.5,92.5) ; |
---|
450 | |
---|
451 | |
---|
452 | // G4VhElectronicStoppingPower* Z77p = new G4hZiegler1977p() ; |
---|
453 | // G4VhElectronicStoppingPower* Z85p = new G4hZiegler1985p() ; |
---|
454 | G4VhElectronicStoppingPower* Z77p = new G4hSRIM2000p() ; |
---|
455 | G4VhElectronicStoppingPower* Z85p = new G4hSRIM2003p() ; |
---|
456 | G4VhElectronicStoppingPower* Z77He = new G4hZiegler1977He() ; |
---|
457 | G4VhElectronicStoppingPower* I49p = new G4hICRU49p() ; |
---|
458 | G4VhElectronicStoppingPower* I49He = new G4hICRU49He() ; |
---|
459 | |
---|
460 | G4double de; |
---|
461 | |
---|
462 | // Ziegler1977p |
---|
463 | G4double tau = 40.0*keV ; |
---|
464 | G4double q2; |
---|
465 | G4double rateMass=4.0026/1.007276; |
---|
466 | |
---|
467 | for ( j=1; j<93; j++) |
---|
468 | { |
---|
469 | z = G4double(j); |
---|
470 | q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin); |
---|
471 | de = Z77p->ElectronicStoppingPower(z, tau) ; |
---|
472 | h[1]->fill(z,de) ; |
---|
473 | de = Z85p->ElectronicStoppingPower(z, tau) ; |
---|
474 | h[55]->fill(z,de) ; |
---|
475 | de = I49p->ElectronicStoppingPower(z, tau) ; |
---|
476 | h[6]->fill(z,de) ; |
---|
477 | de = Z77He->ElectronicStoppingPower(z, tau) ; |
---|
478 | h[11]->fill(z,de) ; |
---|
479 | de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
480 | h[17]->fill(z,de) ; |
---|
481 | de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
482 | h[23]->fill(z,de) ; |
---|
483 | } |
---|
484 | tau = 100.0*keV ; |
---|
485 | for ( j=1; j<93; j++) |
---|
486 | { |
---|
487 | z = G4double(j); |
---|
488 | de = Z77p->ElectronicStoppingPower(z, tau) ; |
---|
489 | h[2]->fill(z,de) ; |
---|
490 | de = Z85p->ElectronicStoppingPower(z, tau) ; |
---|
491 | h[56]->fill(z,de) ; |
---|
492 | de = I49p->ElectronicStoppingPower(z, tau) ; |
---|
493 | h[7]->fill(z,de) ; |
---|
494 | de = Z77He->ElectronicStoppingPower(z, tau) ; |
---|
495 | h[12]->fill(z,de) ; |
---|
496 | q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin); |
---|
497 | de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
498 | h[18]->fill(z,de) ; |
---|
499 | de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
500 | h[24]->fill(z,de) ; |
---|
501 | } |
---|
502 | |
---|
503 | tau = 400.0*keV ; |
---|
504 | for ( j=1; j<93; j++) |
---|
505 | { |
---|
506 | z = G4double(j); |
---|
507 | de = Z77p->ElectronicStoppingPower(z, tau) ; |
---|
508 | h[3]->fill(z,de) ; |
---|
509 | de = Z85p->ElectronicStoppingPower(z, tau) ; |
---|
510 | h[57]->fill(z,de) ; |
---|
511 | de = I49p->ElectronicStoppingPower(z, tau) ; |
---|
512 | h[8]->fill(z,de) ; |
---|
513 | de = Z77He->ElectronicStoppingPower(z, tau) ; |
---|
514 | h[13]->fill(z,de) ; |
---|
515 | q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin); |
---|
516 | de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
517 | h[19]->fill(z,de) ; |
---|
518 | de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
519 | h[25]->fill(z,de) ; |
---|
520 | } |
---|
521 | tau = 1000.0*keV ; |
---|
522 | for ( j=1; j<93; j++) |
---|
523 | { |
---|
524 | z = G4double(j); |
---|
525 | de = Z77p->ElectronicStoppingPower(z, tau) ; |
---|
526 | h[4]->fill(z,de) ; |
---|
527 | de = Z85p->ElectronicStoppingPower(z, tau) ; |
---|
528 | h[58]->fill(z,de) ; |
---|
529 | de = I49p->ElectronicStoppingPower(z, tau) ; |
---|
530 | h[9]->fill(z,de) ; |
---|
531 | de = Z77He->ElectronicStoppingPower(z, tau) ; |
---|
532 | h[14]->fill(z,de) ; |
---|
533 | q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin); |
---|
534 | de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
535 | h[20]->fill(z,de) ; |
---|
536 | de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
537 | h[26]->fill(z,de) ; |
---|
538 | } |
---|
539 | tau = 4000*keV ; |
---|
540 | for ( j=1; j<93; j++) |
---|
541 | { |
---|
542 | z = G4double(j); |
---|
543 | de = Z77p->ElectronicStoppingPower(z, tau) ; |
---|
544 | h[5]->fill(z,de) ; |
---|
545 | de = Z85p->ElectronicStoppingPower(z, tau) ; |
---|
546 | h[59]->fill(z,de) ; |
---|
547 | de = I49p->ElectronicStoppingPower(z, tau) ; |
---|
548 | h[10]->fill(z,de) ; |
---|
549 | de = Z77He->ElectronicStoppingPower(z, tau) ; |
---|
550 | h[15]->fill(z,de) ; |
---|
551 | q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin); |
---|
552 | de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
553 | h[21]->fill(z,de) ; |
---|
554 | de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
555 | h[27]->fill(z,de) ; |
---|
556 | } |
---|
557 | tau = 6.5*MeV ; |
---|
558 | for ( j=1; j<93; j++) |
---|
559 | { |
---|
560 | z = G4double(j); |
---|
561 | de = Z77p->ElectronicStoppingPower(z, tau) ; |
---|
562 | h[70]->fill(z,de) ; |
---|
563 | } |
---|
564 | tau = 10.0*keV ; |
---|
565 | for ( j=1; j<93; j++) |
---|
566 | { |
---|
567 | z = G4double(j); |
---|
568 | q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin); |
---|
569 | de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
570 | h[16]->fill(z,de) ; |
---|
571 | de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
572 | h[22]->fill(z,de) ; |
---|
573 | } |
---|
574 | |
---|
575 | for ( j=1; j<93; j++) |
---|
576 | { |
---|
577 | static G4double p40[92] = { |
---|
578 | 6.22, 6.55, 7.61, 10.4, 12.9, 13.9, 15.9, 14.6, 11.7, 11.1, |
---|
579 | 14.2, 20.6, 20.7, 22.3, 18.1, 19.3, 27.5, 30.5, 27.9, 29.8, |
---|
580 | 28.3, 26.7, 25.1, 22.5, 19.8, 20.1, 18.0, 20.1, 20.4, 23.4, |
---|
581 | 27.5, 29.0, 29.3, 31.0, 31.1, 34.8, 31.5, 35.0, 35.5, 37.3, |
---|
582 | 38.3, 35.9, 38.0, 34.4, 33.4, 29.8, 30.9, 32.7, 34.9, 36.0, |
---|
583 | 40.9, 39.1, 43.1, 45.8, 40.8, 44.1, 44.9, 42.0, 41.0, 40.0, |
---|
584 | 39.0, 38.0, 37.1, 38.2, 35.3, 31.5, 29.9, 29.1, 28.3, 27.5, |
---|
585 | 28.1, 28.9, 27.3, 26.3, 29.9, 29.1, 28.4, 25.8, 28.0, 24.9, |
---|
586 | 27.3, 30.7, 34.3, 35.4, 35.6, 35.5, 39.8, 42.9, 43.7, 44.1, |
---|
587 | 42.4, 41.8} ; |
---|
588 | h[53]->fill(double(j),p40[j-1]) ; |
---|
589 | |
---|
590 | static G4double he40[92] = { |
---|
591 | 11.8, 16.7, 21.9, 24.1, 34.7, 36.1, 45.5, 46.1, 45.8, 44.9, |
---|
592 | 46.9, 51.7, 54.9, 63.2, 63.8, 67.4, 79.8, 80.8, 78.0, 83.2, |
---|
593 | 85.4, 84.2, 90.6, 85.7, 84.1, 86.4, 82.0, 77.7, 74.4, 75.1, |
---|
594 | 75.6, 80.7, 84.9, 87.0, 88.6, 103.0, 103.0, 113.0, 116.0, 123.0, |
---|
595 | 120.0, 115.0, 118.0, 115.0, 113.0, 112.0, 111.0, 110.0, 116.0, 117.0, |
---|
596 | 119.0, 123.0, 122.0, 139.0, 138.0, 143.0, 152.0, 141.0, 139.0, 137.0, |
---|
597 | 135.0, 134.0, 130.0, 137.0, 129.0, 128.0, 124.0, 125.0, 120.0, 118.0, |
---|
598 | 119.0, 120.0, 121.0, 122.0, 125.0, 127.0, 128.0, 120.0, 125.0, 128.0, |
---|
599 | 136.0, 138.0, 144.0, 148.0, 151.0, 162.0, 166.0, 177.0, 179.0, 183.0, |
---|
600 | 177.0, 176.0 } ; |
---|
601 | h[54]->fill(double(j),he40[j-1]) ; |
---|
602 | } |
---|
603 | |
---|
604 | G4cout << "Stopping Power histograms are filled!" << G4endl; |
---|
605 | |
---|
606 | // dedx |
---|
607 | for (j = 0 ; j < num-1 ; j++) { |
---|
608 | tkin = std::pow(10.0,(std::log10(minE) + (G4double(j)+0.5)*s)); |
---|
609 | de = hIon[0]->ComputeDEDX(part[0],Graphite,tkin) ; |
---|
610 | h[28]->fill(std::log10(tkin),de) ; |
---|
611 | de = hIon[0]->ComputeDEDX(part[0],Al,tkin) ; |
---|
612 | h[29]->fill(std::log10(tkin),de) ; |
---|
613 | de = hIon[0]->ComputeDEDX(part[0],Si,tkin) ; |
---|
614 | h[30]->fill(std::log10(tkin),de) ; |
---|
615 | de = hIon[0]->ComputeDEDX(part[0],Cu,tkin) ; |
---|
616 | h[31]->fill(std::log10(tkin),de) ; |
---|
617 | de = hIon[0]->ComputeDEDX(part[0],Fe,tkin) ; |
---|
618 | h[32]->fill(std::log10(tkin),de) ; |
---|
619 | de = hIon[0]->ComputeDEDX(part[0],Pb,tkin) ; |
---|
620 | h[33]->fill(std::log10(tkin),de) ; |
---|
621 | de = hIon[0]->ComputeDEDX(part[0],ethane,tkin) ; |
---|
622 | // G4cout << "ethane: E = " << tkin << "; dedx = " << de << G4endl; |
---|
623 | h[34]->fill(std::log10(tkin),de) ; |
---|
624 | de = hIon[0]->ComputeDEDX(part[0],water,tkin) ; |
---|
625 | // de += theNuclearStoppingModel->TheValue(part[1],water,tkin) ; |
---|
626 | // G4cout << "water : E = " << tkin << "; dedx = " << de << G4endl; |
---|
627 | h[35]->fill(std::log10(tkin),de) ; |
---|
628 | de = hIon[0]->ComputeDEDX(part[0],LAr,tkin) ; |
---|
629 | h[36]->fill(std::log10(tkin),de) ; |
---|
630 | de = hIon[0]->ComputeDEDX(part[0],csi,tkin) ; |
---|
631 | h[37]->fill(std::log10(tkin),de) ; |
---|
632 | } |
---|
633 | |
---|
634 | G4cout << "Proton's dEdx histograms are filled!" << G4endl; |
---|
635 | |
---|
636 | for (j = 0 ; j < num-1 ; j++) { |
---|
637 | tkin = std::pow(10.0,(std::log10(minE) + (G4double(j)+0.5)*s)); |
---|
638 | de = hIon[0]->ComputeDEDX(part[1],Graphite,tkin) ; |
---|
639 | h[60]->fill(std::log10(tkin),de) ; |
---|
640 | de = hIon[0]->ComputeDEDX(part[1],Al,tkin) ; |
---|
641 | h[61]->fill(std::log10(tkin),de) ; |
---|
642 | de = hIon[0]->ComputeDEDX(part[1],Si,tkin) ; |
---|
643 | h[62]->fill(std::log10(tkin),de) ; |
---|
644 | de = hIon[0]->ComputeDEDX(part[1],Cu,tkin) ; |
---|
645 | h[63]->fill(std::log10(tkin),de) ; |
---|
646 | de = hIon[0]->ComputeDEDX(part[1],Fe,tkin) ; |
---|
647 | h[64]->fill(std::log10(tkin),de) ; |
---|
648 | de = hIon[0]->ComputeDEDX(part[1],Pb,tkin) ; |
---|
649 | h[65]->fill(std::log10(tkin),de) ; |
---|
650 | de = hIon[0]->ComputeDEDX(part[1],ethane,tkin) ; |
---|
651 | // G4cout << "ethane: E = " << tkin << "; dedx = " << de << G4endl; |
---|
652 | h[66]->fill(std::log10(tkin),de) ; |
---|
653 | de = hIon[0]->ComputeDEDX(part[1],water,tkin) ; |
---|
654 | // G4cout << "water : E = " << tkin << "; dedx = " << de << G4endl; |
---|
655 | h[67]->fill(std::log10(tkin),de) ; |
---|
656 | de = hIon[0]->ComputeDEDX(part[1],LAr,tkin) ; |
---|
657 | h[68]->fill(std::log10(tkin),de) ; |
---|
658 | de = hIon[0]->ComputeDEDX(part[1],csi,tkin) ; |
---|
659 | h[69]->fill(std::log10(tkin),de) ; |
---|
660 | } |
---|
661 | |
---|
662 | G4cout << "AntiProton's dEdx histograms are filled!" << G4endl; |
---|
663 | |
---|
664 | G4double mProt = part[0]->GetPDGMass()*1.007276; |
---|
665 | G4double fact = cm/(2700.0*MeV) ; // to MeV/mg/cm^2 |
---|
666 | |
---|
667 | for (j = 0 ; j < num-1 ; j++) { |
---|
668 | tkin = std::pow(10.0,(std::log10(minE) + (G4double(j)+0.5)*s)); |
---|
669 | de = theIonEffChargeModel->TheValue(part[3],Cu,tkin) ; |
---|
670 | // G4cout << "E = " << tkin << "; dedx = " << de << G4endl; |
---|
671 | h[48]->fill(std::log10(tkin),de) ; |
---|
672 | de = theIonEffChargeModel->TheValue(part[4],Cu,tkin) ; |
---|
673 | h[49]->fill(std::log10(tkin),de) ; |
---|
674 | G4double tRed = tkin * (part[3]->GetPDGMass())/mProt ; |
---|
675 | de = hIon[3]->ComputeDEDX(part[3],Al,tRed) ; |
---|
676 | de += theNuclearStoppingModel->TheValue(part[3],Al,tRed) ; |
---|
677 | h[50]->fill(std::log10(tkin),de*fact) ; |
---|
678 | tRed = tkin * (part[4]->GetPDGMass())/mProt ; |
---|
679 | de = hIon[4]->ComputeDEDX(part[4],Al,tRed) ; |
---|
680 | //de += theNuclearStoppingModel->TheValue(part[4],Al,tRed) ; |
---|
681 | h[51]->fill(std::log10(tkin),de*fact) ; |
---|
682 | tRed = tkin * (part[5]->GetPDGMass())/mProt ; |
---|
683 | de = hIon[5]->ComputeDEDX(part[5],Al,tRed) ; |
---|
684 | //de += theNuclearStoppingModel->TheValue(part[5],Al,tRed) ; |
---|
685 | h[52]->fill(std::log10(tkin),de*fact) ; |
---|
686 | } |
---|
687 | |
---|
688 | G4cout << "Ions's dEdx histograms are filled!" << G4endl; |
---|
689 | |
---|
690 | theProcessManager[7] = new G4ProcessManager(part[0]); |
---|
691 | part[0]->SetProcessManager(theProcessManager[7]); |
---|
692 | hIon[7] = new G4hLowEnergyIonisation(); |
---|
693 | hIon[7]->SetElectronicStoppingPowerModel(part[0],"Ziegler1985p"); |
---|
694 | hIon[7]->SetEnlossFluc(false) ; |
---|
695 | //hIon[7]->SetNuclearStoppingOff(); |
---|
696 | // hIon[7]->SetBarkasOff(); |
---|
697 | theProcessManager[7]->AddProcess(hIon[7]); |
---|
698 | hIon[7]->BuildPhysicsTable(*part[0]); |
---|
699 | |
---|
700 | G4cout << "Ziegler's dEdx histograms will be filled" << G4endl; |
---|
701 | |
---|
702 | for (j = 0 ; j < num-1 ; j++) { |
---|
703 | tkin = std::pow(10.0,(std::log10(minE) + (G4double(j)+0.5)*s)); |
---|
704 | de = hIon[7]->ComputeDEDX(part[0],Graphite,tkin) ; |
---|
705 | h[38]->fill(std::log10(tkin),de) ; |
---|
706 | de = hIon[7]->ComputeDEDX(part[0],Al,tkin) ; |
---|
707 | h[39]->fill(std::log10(tkin),de) ; |
---|
708 | de = hIon[7]->ComputeDEDX(part[0],Si,tkin) ; |
---|
709 | h[40]->fill(std::log10(tkin),de) ; |
---|
710 | de = hIon[7]->ComputeDEDX(part[0],Cu,tkin) ; |
---|
711 | h[41]->fill(std::log10(tkin),de) ; |
---|
712 | de = hIon[7]->ComputeDEDX(part[0],Fe,tkin) ; |
---|
713 | h[42]->fill(std::log10(tkin),de) ; |
---|
714 | de = hIon[7]->ComputeDEDX(part[0],Pb,tkin) ; |
---|
715 | h[43]->fill(std::log10(tkin),de) ; |
---|
716 | de = hIon[7]->ComputeDEDX(part[0],ethane,tkin) ; |
---|
717 | // G4cout << "ethane: E = " << tkin << "; dedx = " << de << G4endl; |
---|
718 | h[44]->fill(std::log10(tkin),de) ; |
---|
719 | de = hIon[7]->ComputeDEDX(part[0],water,tkin) ; |
---|
720 | // G4cout << "water: E = " << tkin << "; dedx = " << de << G4endl; |
---|
721 | h[45]->fill(std::log10(tkin),de) ; |
---|
722 | de = hIon[7]->ComputeDEDX(part[0],LAr,tkin) ; |
---|
723 | h[46]->fill(std::log10(tkin),de) ; |
---|
724 | de = hIon[7]->ComputeDEDX(part[0],csi,tkin) ; |
---|
725 | h[47]->fill(std::log10(tkin),de) ; |
---|
726 | } |
---|
727 | |
---|
728 | //---------------------- Fill Ntuple ------------------------ |
---|
729 | /* |
---|
730 | G4cout << "Fill Ntuple!" << G4endl; |
---|
731 | |
---|
732 | // ---- primary ntuple ------ |
---|
733 | HepTuple* ntuple = hbookManager->ntuple("dEdx ntuple"); |
---|
734 | assert (ntuple != 0); |
---|
735 | |
---|
736 | for( i = 0; i < num; i++) { |
---|
737 | tkin = std::pow(10,(std::log10(minE) + i*s)); |
---|
738 | |
---|
739 | for ( G4int j = 0 ; j < numOfMaterials; j++ ) { |
---|
740 | |
---|
741 | material = (*theMaterialTable)[j] ; |
---|
742 | // get elements in the actual material, |
---|
743 | const G4ElementVector* theElementVector = material->GetElementVector() ; |
---|
744 | const G4double* theAtomicNumDensityVector = |
---|
745 | material->GetAtomicNumDensityVector() ; |
---|
746 | const G4int NumberOfElements = material->GetNumberOfElements() ; |
---|
747 | |
---|
748 | // loop for the elements in the material |
---|
749 | // to find out average values Z, vF, lF |
---|
750 | G4double z = 0.0, norm = 0.0 ; |
---|
751 | |
---|
752 | if( 1 == NumberOfElements ) { |
---|
753 | z = material->GetZ() ; |
---|
754 | |
---|
755 | } else { |
---|
756 | for (G4int iel=0; iel<NumberOfElements; iel++) { |
---|
757 | const G4Element* element = (*theElementVector)(iel) ; |
---|
758 | G4double z2 = element->GetZ() ; |
---|
759 | const G4double weight = theAtomicNumDensityVector[iel] ; |
---|
760 | norm += weight ; |
---|
761 | z += z2 * weight ; |
---|
762 | } |
---|
763 | z /= norm ; |
---|
764 | } |
---|
765 | |
---|
766 | for (G4int k=0 ; k<7; k++) { |
---|
767 | G4double dedx = hIon[k]->ComputeDEDX(part[k],material,tkin) ; |
---|
768 | G4double q = theIonEffChargeModel->TheValue(part[k],material,tkin); |
---|
769 | ntuple->column("part",k); |
---|
770 | ntuple->column("mat",j); |
---|
771 | ntuple->column("ie",i); |
---|
772 | ntuple->column("zeff",z); |
---|
773 | ntuple->column("tkin",tkin/MeV); |
---|
774 | ntuple->column("mass",(part[k]->GetPDGMass())/MeV); |
---|
775 | ntuple->column("char",(part[k]->GetPDGCharge())/eplus); |
---|
776 | ntuple->column("cha2",q); |
---|
777 | ntuple->column("dedx",dedx*mm/MeV); |
---|
778 | ntuple->dumpData(); |
---|
779 | } |
---|
780 | } |
---|
781 | } |
---|
782 | */ |
---|
783 | //----------- End of work ------------------------------------- |
---|
784 | |
---|
785 | std::cout << "Committing..." << std::endl; |
---|
786 | tree->commit(); |
---|
787 | std::cout << "Closing the tree..." << std::endl; |
---|
788 | tree->close(); |
---|
789 | |
---|
790 | G4cout << "Ntuple and Hbook are saved" << G4endl; |
---|
791 | |
---|
792 | // delete materials and elements |
---|
793 | delete Be; |
---|
794 | delete Graphite; |
---|
795 | delete Al; |
---|
796 | delete Si; |
---|
797 | delete LAr; |
---|
798 | delete Fe; |
---|
799 | delete Cu; |
---|
800 | delete W; |
---|
801 | delete Pb; |
---|
802 | delete U; |
---|
803 | delete H; |
---|
804 | delete C; |
---|
805 | delete Cs; |
---|
806 | delete I; |
---|
807 | delete Ti; |
---|
808 | delete O; |
---|
809 | delete water; |
---|
810 | delete ethane; |
---|
811 | delete csi; |
---|
812 | G4cout << "Materials are deleted" << G4endl; |
---|
813 | delete theIonEffChargeModel; |
---|
814 | delete Z77p; |
---|
815 | delete Z77He; |
---|
816 | delete I49p; |
---|
817 | delete I49He; |
---|
818 | |
---|
819 | cout<<"END OF THE MAIN PROGRAM"<<G4endl; |
---|
820 | } |
---|