source: trunk/source/processes/electromagnetic/lowenergy/test/G4CrossSectionHandlerTest.cc@ 1350

Last change on this file since 1350 was 1350, checked in by garnier, 15 years ago

update to last version 4.9.4
File size: 12.9 KB
RevLine 
[1350]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: G4CrossSectionHandlerTest.cc,v 1.4 2006/06/29 19:43:53 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// -------------------------------------------------------------------
31// GEANT 4 class file --- Copyright CERN 1998
32// CERN Geneva Switzerland
33//
34//
35// File name: G4DataHandlerTest
36//
37// Author: Maria Grazia Pia
38//
39// Creation date: 3 August 2001
40//
41// Modifications:
42//
43// -------------------------------------------------------------------
44
45#include "CLHEP/Hist/TupleManager.h"
46#include "CLHEP/Hist/HBookFile.h"
47#include "CLHEP/Hist/Histogram.h"
48#include "CLHEP/Hist/Tuple.h"
49
50#include "globals.hh"
51#include "G4ios.hh"
52#include <fstream>
53#include <iomanip>
54
55#include "G4CrossSectionHandler.hh"
56#include "G4VEMDataSet.hh"
57#include "G4VDataSetAlgorithm.hh"
58#include "G4LogLogInterpolation.hh"
59#include "G4Material.hh"
60#include "G4MaterialTable.hh"
61
62HepTupleManager* hbookManager;
63
64int main()
65{
66 G4cout.setf( ios::scientific, ios::floatfield );
67
68
69 // ---- HBOOK initialization
70
71
72 hbookManager = new HBookFile("datahandler.hbook", 58);
73 assert (hbookManager != 0);
74
75 // ---- primary ntuple ------
76 HepTuple* ntuple = hbookManager->ntuple("CrossSectionNtuple");
77 assert (ntuple != 0);
78
79 // Create some materials
80 G4Material* Be = new G4Material("Beryllium", 4., 9.01*g/mole, 1.848*g/cm3);
81 G4Material* Graphite = new G4Material("Graphite",6., 12.00*g/mole, 2.265*g/cm3 );
82 G4Material* Al = new G4Material("Aluminium", 13., 26.98*g/mole, 2.7 *g/cm3);
83 G4Material* Si = new G4Material("Silicon", 14., 28.055*g/mole, 2.33*g/cm3);
84 G4Material* LAr = new G4Material("LArgon", 18., 39.95*g/mole, 1.393*g/cm3);
85 G4Material* Fe = new G4Material("Iron", 26., 55.85*g/mole, 7.87*g/cm3);
86 G4Material* Cu = new G4Material("Copper", 29., 63.55*g/mole, 8.96*g/cm3);
87 G4Material* W = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3);
88 G4Material* Pb = new G4Material("Lead", 82., 207.19*g/mole, 11.35*g/cm3);
89 G4Material* U = new G4Material("Uranium", 92., 238.03*g/mole, 18.95*g/cm3);
90
91 G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
92 G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
93 G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
94 G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
95 G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
96
97 G4Material* maO = new G4Material("Oxygen", 8., 16.00*g/mole, 1.1*g/cm3);
98
99 G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
100 water->AddElement(H,2);
101 water->AddElement(O,1);
102
103 G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
104 ethane->AddElement(H,6);
105 ethane->AddElement(C,2);
106
107 G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
108 csi->AddElement(Cs,1);
109 csi->AddElement(I,1);
110
111 static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
112
113 G4int nMaterials = G4Material::GetNumberOfMaterials();
114 G4cout << "The MaterialTable contains "
115 << nMaterials
116 << " materials "
117 << G4endl;
118
119 G4cout << "G4DataSetManager test: dump LLNL data sets" << G4endl;
120
121 // G4cout << "Enter file name " << G4endl;
122
123 G4VDataSetAlgorithm* interpolation = new G4LogLogInterpolation();
124 G4cout << "Interpolation created" << G4endl;
125
126 G4CrossSectionHandler* manager = new G4CrossSectionHandler(interpolation,10*eV);
127 G4cout << "CrossSectionHandler created" << G4endl;
128
129 G4cout << " Load Atom (1) or Shell (2) data?" << G4endl;
130 G4int type;
131 G4cin >> type;
132 G4String fileName;
133 if (type == 1)
134 {
135 // G4cin >> fileName;
136 fileName = "brem/br-cs-";
137 G4String name(fileName);
138 manager->LoadData(fileName);
139 G4cout << "Data Loaded" << G4endl;
140 }
141 if (type == 2)
142 {
143 // G4cin >> fileName;
144 fileName = "phot/pe-ss-cs-";
145 G4String name(fileName);
146 manager->LoadShellData(fileName);
147 G4cout << "Data Loaded" << G4endl;
148 }
149
150 G4cout << "Data loaded; print (1) or continue (2)? " << G4endl;
151 G4int printData;
152 G4cin >> printData;
153 if (printData == 1) manager->PrintData();
154
155 G4cout << "Enter Z " << G4endl;
156 G4int Z;
157 G4cin >> Z;
158
159 G4cout << "Enter energy" << G4endl;
160 G4double e;
161 G4cin >> e;
162 G4double sigma = manager->FindValue(Z,e) / barn;
163 G4cout << "Value retrieved by manager = " << sigma << G4endl;
164
165 G4cout << "Available materials are: " << G4endl;
166 for (G4int mat = 0; mat < nMaterials; mat++)
167 {
168 G4cout << mat << ") "
169 << (*theMaterialTable)[mat]->GetName()
170 << G4endl;
171 }
172
173 G4int materialId;
174 G4cout << "Which material? " << G4endl;
175 G4cin >> materialId;
176 G4Material* material = (*theMaterialTable)[materialId] ;
177
178 G4double materialSigma = manager->ValueForMaterial(material,e);
179 G4cout << "Material value calculated by manager = " << materialSigma << G4endl;
180
181 G4int selectedZ = manager->SelectRandomAtom(material,e);
182 G4cout << "Select random atom: Z = " << selectedZ << G4endl;
183
184 G4int selectedShell = manager->SelectRandomShell(Z,e);
185 G4cout << "Select random shell: shell = " << selectedShell << G4endl;
186
187 G4cout << "Now clearing" << G4endl;
188 manager->Clear();
189
190 G4cout << "Now re-loading" << G4endl;
191 if (type == 1) manager->LoadData(fileName);
192 if (type == 2) manager->LoadShellData(fileName);
193 G4cout << "End of re-loading" << G4endl;
194
195 G4VEMDataSet* meanFreePathTable = manager->BuildMeanFreePathForMaterials();
196 size_t nRows = meanFreePathTable->NumberOfComponents();
197 meanFreePathTable->PrintData();
198 G4cout << "MeanFreePathTable has " << nRows << " components" << G4endl;
199
200 const G4VEMDataSet* materialSet = meanFreePathTable->GetComponent[materialId];
201
202 G4int i;
203 for (i=2; i<10; i++)
204 {
205 G4double e = (10. * i) * eV;
206 G4double value = manager->FindValue(Z,e) / barn;
207 G4double valueM = manager->ValueForMaterial(material,e);
208 G4double valueT = materialSet->FindValue(e);
209 if (valueT > 0.) valueT = 1. / valueT;
210 ntuple->column("e",e);
211 ntuple->column("sigmac",valueM);
212 ntuple->column("sigmat",valueT);
213 ntuple->dumpData();
214
215 // G4cout << "e = " << e/MeV << " --- Cross section = " << value << " b" << G4endl;
216 }
217 for (i=1; i<10; i++)
218 {
219 G4double e = (100. * i) * eV;
220 G4double value = manager->FindValue(Z,e) / barn;
221 // G4cout << "e = " << e/MeV << " --- Cross section = " << value << " b" << G4endl;
222 G4double valueM = manager->ValueForMaterial(material,e);
223 G4double valueT = materialSet->FindValue(e);
224 if (valueT > 0.) valueT = 1. / valueT;
225 ntuple->column("e",e);
226 ntuple->column("sigmac",valueM);
227 ntuple->column("sigmat",valueT);
228 ntuple->dumpData();
229 }
230 for (i=1; i<10; i++)
231 {
232 G4double e = (1. * i) * keV;
233 G4double value = manager->FindValue(Z,e) / barn;
234 G4double valueM = manager->ValueForMaterial(material,e);
235 G4double valueT = materialSet->FindValue(e);
236 if (valueT > 0.) valueT = 1. / valueT;
237 ntuple->column("e",e);
238 ntuple->column("sigmac",valueM);
239 ntuple->column("sigmat",valueT);
240 ntuple->dumpData();
241 // G4cout << "e = " << e/MeV << " --- Cross section = " << value << " b" << G4endl;
242 }
243 for (i=1; i<10; i++)
244 {
245 G4double e = (10. * i) * keV;
246 G4double value = manager->FindValue(Z,e) / barn;
247 G4double valueM = manager->ValueForMaterial(material,e);
248 G4double valueT = materialSet->FindValue(e);
249 if (valueT > 0.) valueT = 1. / valueT;
250 ntuple->column("e",e);
251 ntuple->column("sigmac",valueM);
252 ntuple->column("sigmat",valueT);
253 ntuple->dumpData();
254
255 // G4cout << "e = " << e/MeV << " --- Cross section = " << value << " b" << G4endl;
256 }
257 for (i=1; i<10; i++)
258 {
259 G4double e = (100. * i) * keV;
260 G4double value = manager->FindValue(Z,e) / barn;
261 G4double valueM = manager->ValueForMaterial(material,e);
262 G4double valueT = materialSet->FindValue(e);
263 if (valueT > 0.) valueT = 1. / valueT;
264 ntuple->column("e",e);
265 ntuple->column("sigmac",valueM);
266 ntuple->column("sigmat",valueT);
267 ntuple->dumpData();
268
269 // G4cout << "e = " << e/MeV << " --- Cross section = " << value << " b" << G4endl;
270 }
271 for (i=1; i<10; i++)
272 {
273 G4double e = (1. * i) * MeV;
274 G4double value = manager->FindValue(Z,e) / barn;
275 G4double valueM = manager->ValueForMaterial(material,e);
276 G4double valueT = materialSet->FindValue(e);
277 if (valueT > 0.) valueT = 1. / valueT;
278 ntuple->column("e",e);
279 ntuple->column("sigmac",valueM);
280 ntuple->column("sigmat",valueT);
281 ntuple->dumpData();
282
283 // G4cout << "e = " << e/MeV << " --- Cross section = " << value << " b" << G4endl;
284 }
285 for (i=1; i<10; i++)
286 {
287 G4double e = (10. * i) * MeV;
288 G4double value = manager->FindValue(Z,e) / barn;
289 G4double valueM = manager->ValueForMaterial(material,e);
290 G4double valueT = materialSet->FindValue(e);
291 if (valueT > 0.) valueT = 1. / valueT;
292 ntuple->column("e",e);
293 ntuple->column("sigmac",valueM);
294 ntuple->column("sigmat",valueT);
295 ntuple->dumpData();
296 // G4cout << "e = " << e/MeV << " --- Cross section = " << value << " b" << G4endl;
297 }
298 for (i=1; i<10; i++)
299 {
300 G4double e = (100. * i) * MeV;
301 G4double value = manager->FindValue(Z,e) / barn;
302 G4double valueM = manager->ValueForMaterial(material,e);
303 G4double valueT = materialSet->FindValue(e);
304 if (valueT > 0.) valueT = 1. / valueT;
305 ntuple->column("e",e);
306 ntuple->column("sigmac",valueM);
307 ntuple->column("sigmat",valueT);
308 ntuple->dumpData();
309
310 // G4cout << "e = " << e/MeV << " --- Cross section = " << value << " b" << G4endl;
311 }
312 for (i=1; i<10; i++)
313 {
314 G4double e = (1. * i) * GeV;
315 G4double value = manager->FindValue(Z,e) / barn;
316 G4double valueM = manager->ValueForMaterial(material,e);
317 G4double valueT = materialSet->FindValue(e);
318 if (valueT > 0.) valueT = 1. / valueT;
319 ntuple->column("e",e);
320 ntuple->column("sigmac",valueM);
321 ntuple->column("sigmat",valueT);
322 ntuple->dumpData();
323
324 // G4cout << "e = " << e/MeV << " --- Cross section = " << value << " b" << G4endl;
325 }
326 for (i=1; i<10; i++)
327 {
328 G4double e = (10. * i) * GeV;
329 G4double value = manager->FindValue(Z,e) / barn;
330 G4double valueM = manager->ValueForMaterial(material,e);
331 G4double valueT = materialSet->FindValue(e);
332 if (valueT > 0.) valueT = 1. / valueT;
333 ntuple->column("e",e);
334 ntuple->column("sigmac",valueM);
335 ntuple->column("sigmat",valueT);
336 ntuple->dumpData();
337 // G4cout << "e = " << e/MeV << " --- Cross section = " << value << " b" << G4endl;
338 }
339
340 G4double eFin = 100 * GeV;
341 G4double value = manager->FindValue(Z,eFin) / barn;
342 G4double valueM = manager->ValueForMaterial(material,eFin);
343 G4double valueT = materialSet->FindValue(eFin);
344 if (valueT > 0.) valueT = 1. / valueT;
345 ntuple->column("e",eFin);
346 ntuple->column("sigmac",valueM);
347 ntuple->column("sigmat",valueT);
348 ntuple->dumpData();
349
350 // G4cout << "e = " << e/MeV << " --- Cross section = " << value << " b" << G4endl;
351
352 hbookManager->write();
353 delete hbookManager;
354
355 delete meanFreePathTable;
356 delete manager;
357
358 cout << "END OF THE MAIN PROGRAM" << G4endl;
359}
360
361
362
363
364
365
366
367
Note: See TracBrowser for help on using the repository browser.