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 | // $Id: G4DNATest.cc,v 1.14 2006/06/29 19:43:57 gunter Exp $ |
---|
27 | // GEANT4 tag $Name: geant4-09-03-cand-01 $ |
---|
28 | |
---|
29 | #include "globals.hh" |
---|
30 | #include "G4ios.hh" |
---|
31 | #include <fstream> |
---|
32 | #include <iomanip> |
---|
33 | #include <memory> |
---|
34 | #include <cstdlib> |
---|
35 | |
---|
36 | #include "G4ParticleDefinition.hh" |
---|
37 | #include "G4ParticleTypes.hh" |
---|
38 | #include "G4ParticleTable.hh" |
---|
39 | #include "G4Material.hh" |
---|
40 | #include "G4MaterialTable.hh" |
---|
41 | #include "G4ProcessManager.hh" |
---|
42 | |
---|
43 | #include "G4EnergyLossTables.hh" |
---|
44 | #include "G4VParticleChange.hh" |
---|
45 | #include "G4ParticleChange.hh" |
---|
46 | #include "G4DynamicParticle.hh" |
---|
47 | #include "G4ForceCondition.hh" |
---|
48 | |
---|
49 | #include "G4LowEnergyBremsstrahlung.hh" |
---|
50 | #include "G4LowEnergyIonisation.hh" |
---|
51 | #include "G4eIonisation.hh" |
---|
52 | #include "G4MultipleScattering.hh" |
---|
53 | #include "G4eIonisation.hh" |
---|
54 | #include "G4eBremsstrahlung.hh" |
---|
55 | #include "G4eplusAnnihilation.hh" |
---|
56 | |
---|
57 | #include "G4ComptonScattering.hh" |
---|
58 | #include "G4PhotoElectricEffect.hh" |
---|
59 | |
---|
60 | #include "G4RunManager.hh" |
---|
61 | |
---|
62 | #include "G4Electron.hh" |
---|
63 | #include "G4Positron.hh" |
---|
64 | #include "G4Gamma.hh" |
---|
65 | |
---|
66 | #include "G4GRSVolume.hh" |
---|
67 | #include "G4Box.hh" |
---|
68 | #include "G4PVPlacement.hh" |
---|
69 | #include "G4Step.hh" |
---|
70 | #include "G4ProductionCutsTable.hh" |
---|
71 | #include "G4MaterialCutsCouple.hh" |
---|
72 | |
---|
73 | #include "G4UnitsTable.hh" |
---|
74 | |
---|
75 | #include "AIDA/IManagedObject.h" |
---|
76 | #include "AIDA/IAnalysisFactory.h" |
---|
77 | #include "AIDA/ITreeFactory.h" |
---|
78 | #include "AIDA/ITree.h" |
---|
79 | #include "AIDA/IHistogramFactory.h" |
---|
80 | #include "AIDA/IHistogram1D.h" |
---|
81 | #include "AIDA/IHistogram2D.h" |
---|
82 | #include "AIDA/IHistogram3D.h" |
---|
83 | #include "AIDA/ITupleFactory.h" |
---|
84 | #include "AIDA/ITuple.h" |
---|
85 | |
---|
86 | // DNA |
---|
87 | #include "G4DNAGenericIonsManager.hh" |
---|
88 | #include "G4DNAElectronElasticBrenner.hh" |
---|
89 | #include "G4DNAElectronElasticEmfietzoglou.hh" |
---|
90 | #include "G4DNAProtonExcitation.hh" |
---|
91 | #include "G4DNAHeliumExcitation.hh" |
---|
92 | #include "G4DNAAlphaPlusExcitation.hh" |
---|
93 | #include "G4DNAAlphaPlusPlusExcitation.hh" |
---|
94 | #include "G4DNAElectronExcitation.hh" |
---|
95 | #include "G4DNAProtonRuddIonization.hh" |
---|
96 | #include "G4DNAProtonChargeDecrease.hh" |
---|
97 | #include "G4DNAHydrogenChargeIncrease.hh" |
---|
98 | #include "G4DNAHydrogenRuddIonization.hh" |
---|
99 | #include "G4DNAProtonBornExcitation.hh" |
---|
100 | #include "G4DNAElectronBornExcitation.hh" |
---|
101 | #include "G4DNAAlphaPlusChargeDecrease.hh" |
---|
102 | #include "G4DNAAlphaPlusChargeIncrease.hh" |
---|
103 | #include "G4DNAAlphaPlusPlusChargeDecrease.hh" |
---|
104 | #include "G4DNAHeliumChargeIncrease.hh" |
---|
105 | |
---|
106 | //! \brief Options structure |
---|
107 | struct Options |
---|
108 | { |
---|
109 | //! \brief Mean free path test |
---|
110 | bool meanFreePathTest; |
---|
111 | //! \brief Post step do it test |
---|
112 | bool postStepDoItTest; |
---|
113 | //! \brief Post step do it test |
---|
114 | bool randomEnergy; |
---|
115 | //! \brief Output file name |
---|
116 | const char *outputFileName; |
---|
117 | //! \brief Material name |
---|
118 | const char *material; |
---|
119 | //! \brief Process name |
---|
120 | const char *process; |
---|
121 | //! \brief Particle name |
---|
122 | const char *particle; |
---|
123 | //! \brief Minimum energy |
---|
124 | G4double minEnergy; |
---|
125 | //! \brief Maximum energy |
---|
126 | G4double maxEnergy; |
---|
127 | //! \brief Number of energy step |
---|
128 | G4int nEnergySteps; |
---|
129 | //! \brief Number of interactions |
---|
130 | G4int nIterations; |
---|
131 | }; |
---|
132 | |
---|
133 | //! \brief Default output file name |
---|
134 | struct Options defaultOptions = { false, false, false, "G4DNATest.hbook", "Water", "G4DNAElectronElasticBrenner", "electron", 7.5*eV, 200*eV, 300, 1 }; |
---|
135 | |
---|
136 | //! \brief Creates some materials |
---|
137 | void CreateMaterials(void) |
---|
138 | { |
---|
139 | G4Element * H = new G4Element ("Hydrogen", "H", 1., 1.01*g/mole); |
---|
140 | G4Element * O = new G4Element ("Oxygen", "O", 8., 16.00*g/mole); |
---|
141 | G4Element * C = new G4Element ("Carbon", "C", 6., 12.00*g/mole); |
---|
142 | G4Element * Cs = new G4Element ("Cesium", "Cs", 55., 132.905*g/mole); |
---|
143 | G4Element * I = new G4Element ("Iodine", "I", 53., 126.9044*g/mole); |
---|
144 | |
---|
145 | G4Material * Si = new G4Material("Silicon", 14., 28.055*g/mole, 2.33*g/cm3); |
---|
146 | G4Material * Fe = new G4Material("Iron", 26., 55.85*g/mole, 7.87*g/cm3); |
---|
147 | G4Material * Cu = new G4Material("Copper", 29., 63.55*g/mole, 8.96*g/cm3); |
---|
148 | G4Material * W = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3); |
---|
149 | G4Material * Pb = new G4Material("Lead", 82., 207.19*g/mole, 11.35*g/cm3); |
---|
150 | G4Material * U = new G4Material("Uranium", 92., 238.03*g/mole, 18.95*g/cm3); |
---|
151 | G4Material * maO = new G4Material("Oxygen", 8., 16.00*g/mole, 1.1*g/cm3); |
---|
152 | G4Material * water = new G4Material ("Water", 1.*g/cm3, 2); |
---|
153 | water->AddElement(H, 2); |
---|
154 | water->AddElement(O, 1); |
---|
155 | |
---|
156 | G4Material* ethane = new G4Material ("Ethane", 0.4241*g/cm3, 2); |
---|
157 | ethane->AddElement(H, 6); |
---|
158 | ethane->AddElement(C, 2); |
---|
159 | |
---|
160 | G4Material* csI = new G4Material ("CsI", 4.53*g/cm3, 2); |
---|
161 | csI->AddElement(Cs, 1); |
---|
162 | csI->AddElement(I, 1); |
---|
163 | |
---|
164 | // This is needed to suppress some warnings. These lines can be deleted; |
---|
165 | Si->GetTemperature(); |
---|
166 | Fe->GetTemperature(); |
---|
167 | Cu->GetTemperature(); |
---|
168 | W->GetTemperature(); |
---|
169 | Pb->GetTemperature(); |
---|
170 | U->GetTemperature(); |
---|
171 | maO->GetTemperature(); |
---|
172 | water->GetTemperature(); |
---|
173 | ethane->GetTemperature(); |
---|
174 | csI->GetTemperature(); |
---|
175 | } |
---|
176 | |
---|
177 | //! \brief Process the options arguments |
---|
178 | //! \param argc Number of arguments |
---|
179 | //! \param argv Pointer to the arguments |
---|
180 | //! \param options Structure to fill-in |
---|
181 | void processOptions(int argc, char ** argv, struct Options * options) |
---|
182 | { |
---|
183 | options->meanFreePathTest = defaultOptions.meanFreePathTest; |
---|
184 | options->postStepDoItTest = defaultOptions.meanFreePathTest; |
---|
185 | options->randomEnergy = defaultOptions.randomEnergy; |
---|
186 | options->outputFileName = defaultOptions.outputFileName; |
---|
187 | options->material = defaultOptions.material; |
---|
188 | options->process = defaultOptions.process; |
---|
189 | options->particle = defaultOptions.particle; |
---|
190 | options->minEnergy = defaultOptions.minEnergy; |
---|
191 | options->maxEnergy = defaultOptions.maxEnergy; |
---|
192 | options->nEnergySteps = defaultOptions.nEnergySteps; |
---|
193 | options->nIterations = defaultOptions.nIterations; |
---|
194 | |
---|
195 | int i(1); |
---|
196 | |
---|
197 | while (i<argc) |
---|
198 | { |
---|
199 | if (argv[i][0]=='-' && argv[i][2]==0) |
---|
200 | { |
---|
201 | switch(argv[i][1]) |
---|
202 | { |
---|
203 | case 'h': |
---|
204 | case '?': |
---|
205 | G4cout << argv[0] << " [-h|-?] [-a] [-b] [-r] [-o <file name>] [-m <material name>] [-p <process name>] [-P <particle name>] [-e <min energy in eV>] [-E <max energy in eV>] [-s <energy steps>] [-n <iterations>] " << G4endl |
---|
206 | << G4endl |
---|
207 | << "-h|-? Shows this help" << G4endl |
---|
208 | << "-a Enables mean free path test" << G4endl |
---|
209 | << "-b Enables post step do it test" << G4endl |
---|
210 | << "-r Energy is choosen at random within the range" << G4endl |
---|
211 | << "-o <arg> Set the output file name (default: \"" << defaultOptions.outputFileName << "\")" << G4endl |
---|
212 | << "-m <arg> Set the material (default: \"" << defaultOptions.material << "\")" << G4endl |
---|
213 | << "-p <arg> Set the process (default: \"" << defaultOptions.process << "\")" << G4endl |
---|
214 | << "-P <arg> Set the incoming particle (default: \"" << defaultOptions.particle << "\")" << G4endl |
---|
215 | << "-e <arg> Set the low energy range in eV (default: " << defaultOptions.minEnergy/eV << " eV)" << G4endl |
---|
216 | << "-E <arg> Set the high energy range in eV (default: " << defaultOptions.maxEnergy/eV << " eV)" << G4endl |
---|
217 | << "-s <arg> Set the energy range step (default: " << defaultOptions.nEnergySteps << ")"<< G4endl |
---|
218 | << "-n <arg> Set the number of iterations for the post step do it (default: " << defaultOptions.nIterations << ")" << G4endl; |
---|
219 | exit(0); |
---|
220 | break; |
---|
221 | |
---|
222 | case 'a': |
---|
223 | options->meanFreePathTest=true; |
---|
224 | break; |
---|
225 | |
---|
226 | case 'b': |
---|
227 | options->postStepDoItTest=true; |
---|
228 | break; |
---|
229 | |
---|
230 | case 'r': |
---|
231 | options->randomEnergy=true; |
---|
232 | break; |
---|
233 | |
---|
234 | case 'o': |
---|
235 | i++; |
---|
236 | if (i<argc) |
---|
237 | { |
---|
238 | options->outputFileName = argv[i]; |
---|
239 | break; |
---|
240 | } |
---|
241 | |
---|
242 | case 'm': |
---|
243 | i++; |
---|
244 | if (i<argc) |
---|
245 | { |
---|
246 | options->material = argv[i]; |
---|
247 | break; |
---|
248 | } |
---|
249 | |
---|
250 | case 'p': |
---|
251 | i++; |
---|
252 | if (i<argc) |
---|
253 | { |
---|
254 | options->process = argv[i]; |
---|
255 | break; |
---|
256 | } |
---|
257 | |
---|
258 | case 'P': |
---|
259 | i++; |
---|
260 | if (i<argc) |
---|
261 | { |
---|
262 | options->particle = argv[i]; |
---|
263 | break; |
---|
264 | } |
---|
265 | |
---|
266 | case 'e': |
---|
267 | i++; |
---|
268 | if (i<argc) |
---|
269 | { |
---|
270 | options->minEnergy = std::atof(argv[i])*eV; |
---|
271 | if (options->minEnergy <= 0.) |
---|
272 | { |
---|
273 | G4cout << argv[0] << ": Energy must be > 0." << G4endl; |
---|
274 | exit(-1); |
---|
275 | } |
---|
276 | |
---|
277 | break; |
---|
278 | } |
---|
279 | |
---|
280 | case 'E': |
---|
281 | i++; |
---|
282 | if (i<argc) |
---|
283 | { |
---|
284 | options->maxEnergy = std::atof(argv[i])*eV; |
---|
285 | if (options->maxEnergy <= 0.) |
---|
286 | { |
---|
287 | G4cout << argv[0] << ": Energy must be > 0." << G4endl; |
---|
288 | exit(-1); |
---|
289 | } |
---|
290 | |
---|
291 | break; |
---|
292 | } |
---|
293 | |
---|
294 | case 's': |
---|
295 | i++; |
---|
296 | if (i<argc) |
---|
297 | { |
---|
298 | options->nEnergySteps = atoi(argv[i]); |
---|
299 | if (options->nEnergySteps <= 1) |
---|
300 | { |
---|
301 | G4cout << argv[0] << ": Expected at least two steps." << G4endl; |
---|
302 | exit(-1); |
---|
303 | } |
---|
304 | |
---|
305 | break; |
---|
306 | } |
---|
307 | |
---|
308 | G4cout << argv[0] << ": Expected one more parameter in " << argv[i] << " option. Use -h option for help." << G4endl; |
---|
309 | exit(-1); |
---|
310 | |
---|
311 | case 'n': |
---|
312 | i++; |
---|
313 | if (i<argc) |
---|
314 | { |
---|
315 | options->nIterations = atoi(argv[i]); |
---|
316 | if (options->nIterations <= 0) |
---|
317 | { |
---|
318 | G4cout << argv[0] << ": Expected at least one iteration." << G4endl; |
---|
319 | exit(-1); |
---|
320 | } |
---|
321 | |
---|
322 | break; |
---|
323 | } |
---|
324 | |
---|
325 | G4cout << argv[0] << ": Expected one more parameter in " << argv[i] << " option. Use -h option for help." << G4endl; |
---|
326 | exit(-1); |
---|
327 | |
---|
328 | default: |
---|
329 | G4cout << argv[0] << ": Unknown " << argv[i] << " option. Use -h option for help." << G4endl; |
---|
330 | exit(-1); |
---|
331 | } |
---|
332 | } |
---|
333 | else |
---|
334 | { |
---|
335 | G4cout << argv[0] << ": Bad arguments. Use -h option for help." << G4endl; |
---|
336 | exit(-1); |
---|
337 | } |
---|
338 | |
---|
339 | i++; |
---|
340 | } |
---|
341 | |
---|
342 | if (options->minEnergy >= options->maxEnergy) |
---|
343 | { |
---|
344 | G4cout << argv[0] << ": Mininum energy is higher than maximum energy" << G4endl; |
---|
345 | exit(-1); |
---|
346 | } |
---|
347 | |
---|
348 | G4cout << "Mean free path test: "; |
---|
349 | if (options->meanFreePathTest) |
---|
350 | G4cout << "On"; |
---|
351 | else |
---|
352 | G4cout << "Off"; |
---|
353 | G4cout << G4endl << "Post step do it test: "; |
---|
354 | if (options->postStepDoItTest) |
---|
355 | G4cout << "On"; |
---|
356 | else |
---|
357 | G4cout << "Off"; |
---|
358 | G4cout << G4endl << "Random energy generation: "; |
---|
359 | if (options->randomEnergy) |
---|
360 | G4cout << "On"; |
---|
361 | else |
---|
362 | G4cout << "Off"; |
---|
363 | G4cout << G4endl << "Output file: " << options->outputFileName << G4endl; |
---|
364 | G4cout << "Material: " << options->material << G4endl; |
---|
365 | G4cout << "Process: " << options->process << G4endl; |
---|
366 | G4cout << "Min energy: " << options->minEnergy/eV << " eV" << G4endl; |
---|
367 | G4cout << "Max energy: " << options->maxEnergy/eV << " eV" << G4endl; |
---|
368 | G4cout << "N energy steps: " << options->nEnergySteps << G4endl; |
---|
369 | G4cout << "N iterations: " << options->nIterations << G4endl; |
---|
370 | } |
---|
371 | |
---|
372 | //! \brief Return the selected material |
---|
373 | //! \param options Options for the material choice |
---|
374 | //! \return The material |
---|
375 | G4Material * GetSelectedMaterial(const struct Options & options) |
---|
376 | { |
---|
377 | const G4MaterialTable* theMaterialTable=G4Material::GetMaterialTable(); |
---|
378 | |
---|
379 | G4int i(G4Material::GetNumberOfMaterials()); |
---|
380 | |
---|
381 | while (i>0) |
---|
382 | { |
---|
383 | i--; |
---|
384 | |
---|
385 | if ((*theMaterialTable)[i]->GetName()==options.material) |
---|
386 | return (*theMaterialTable)[i]; |
---|
387 | } |
---|
388 | |
---|
389 | i=G4Material::GetNumberOfMaterials(); |
---|
390 | |
---|
391 | G4cout << "Available materials are: " << G4endl; |
---|
392 | while (i>0) |
---|
393 | { |
---|
394 | i--; |
---|
395 | G4cout << (*theMaterialTable)[i]->GetName(); |
---|
396 | |
---|
397 | if (i>0) |
---|
398 | G4cout << ", "; |
---|
399 | } |
---|
400 | |
---|
401 | G4cout << G4endl; |
---|
402 | |
---|
403 | exit(-2); |
---|
404 | return 0; |
---|
405 | } |
---|
406 | |
---|
407 | //! \brief Creates the geometry |
---|
408 | //! \param options Options for the material choice |
---|
409 | //! \return The world volume |
---|
410 | G4PVPlacement * CreateGeometry(const struct Options & options) |
---|
411 | { |
---|
412 | G4Box* theFrame = new G4Box ("Frame", 1*mm, 1*mm, 1*mm); |
---|
413 | |
---|
414 | G4LogicalVolume* logicalFrame = new G4LogicalVolume(theFrame, GetSelectedMaterial(options), "LFrame", 0, 0, 0); |
---|
415 | |
---|
416 | G4PVPlacement * placement = new G4PVPlacement(0, G4ThreeVector(), "PFrame", logicalFrame, 0, false, 0); |
---|
417 | |
---|
418 | G4cout << "[OK] Geometry built" << G4endl; |
---|
419 | return placement; |
---|
420 | } |
---|
421 | |
---|
422 | //! \brief Get process from options |
---|
423 | //! \param options Options for the process choice |
---|
424 | //! \return The choosen process |
---|
425 | G4VLowEnergyTestableDiscreteProcess * GetSelectedProcess(const struct Options & options) |
---|
426 | { |
---|
427 | static G4VLowEnergyTestableDiscreteProcess ** processes=0; |
---|
428 | if (!processes) |
---|
429 | { |
---|
430 | G4DNAGenericIonsManager * genericIonsManager; |
---|
431 | genericIonsManager=G4DNAGenericIonsManager::Instance(); |
---|
432 | |
---|
433 | processes=new G4VLowEnergyTestableDiscreteProcess * [18]; |
---|
434 | processes[0]=new G4DNAElectronElasticBrenner; |
---|
435 | processes[1]=new G4DNAElectronElasticEmfietzoglou; |
---|
436 | processes[2]=new G4DNAProtonExcitation; |
---|
437 | processes[3]=new G4DNAHeliumExcitation; |
---|
438 | processes[4]=new G4DNAAlphaPlusExcitation; |
---|
439 | processes[5]=new G4DNAAlphaPlusPlusExcitation; |
---|
440 | processes[6]=new G4DNAElectronExcitation; |
---|
441 | processes[7]=new G4DNAProtonRuddIonization; |
---|
442 | processes[8]=new G4DNAProtonChargeDecrease; |
---|
443 | processes[9]=new G4DNAHydrogenChargeIncrease; |
---|
444 | processes[10]=new G4DNAHydrogenRuddIonization; |
---|
445 | processes[11]=new G4DNAProtonBornExcitation; |
---|
446 | processes[12]=new G4DNAElectronBornExcitation; |
---|
447 | processes[13]=new G4DNAAlphaPlusChargeDecrease; |
---|
448 | processes[14]=new G4DNAAlphaPlusChargeIncrease; |
---|
449 | processes[15]=new G4DNAAlphaPlusPlusChargeDecrease; |
---|
450 | processes[16]=new G4DNAHeliumChargeIncrease; |
---|
451 | processes[17]=0; |
---|
452 | } |
---|
453 | |
---|
454 | unsigned long i(0); |
---|
455 | while (processes[i]) |
---|
456 | { |
---|
457 | if (processes[i]->GetProcessName()==options.process) |
---|
458 | return processes[i]; |
---|
459 | |
---|
460 | i++; |
---|
461 | } |
---|
462 | |
---|
463 | G4cout << "Available processes are: " << G4endl; |
---|
464 | i=0; |
---|
465 | while (processes[i]) |
---|
466 | { |
---|
467 | G4cout << processes[i]->GetProcessName(); |
---|
468 | i++; |
---|
469 | |
---|
470 | if (processes[i]) |
---|
471 | G4cout << ", "; |
---|
472 | } |
---|
473 | |
---|
474 | G4cout << G4endl; |
---|
475 | |
---|
476 | exit(-2); |
---|
477 | return 0; |
---|
478 | } |
---|
479 | |
---|
480 | G4ParticleDefinition * const * ParticleList() |
---|
481 | { |
---|
482 | static G4ParticleDefinition * * particles=0; |
---|
483 | |
---|
484 | if (!particles) |
---|
485 | { |
---|
486 | G4DNAGenericIonsManager * genericIonsManager; |
---|
487 | genericIonsManager=G4DNAGenericIonsManager::Instance(); |
---|
488 | |
---|
489 | particles=new G4ParticleDefinition * [10]; |
---|
490 | particles[0]=G4Electron::Electron(); |
---|
491 | particles[1]=G4Positron::Positron(); |
---|
492 | particles[2]=G4Gamma::Gamma(); |
---|
493 | particles[3]=G4AntiProton::AntiProton(); |
---|
494 | particles[4]=G4Proton::Proton(); |
---|
495 | particles[5]=genericIonsManager->GetIon("alpha++"); |
---|
496 | particles[6]=genericIonsManager->GetIon("alpha+"); |
---|
497 | particles[7]=genericIonsManager->GetIon("helium"); |
---|
498 | particles[8]=genericIonsManager->GetIon("hydrogen"); |
---|
499 | particles[9]=0; |
---|
500 | } |
---|
501 | |
---|
502 | return particles; |
---|
503 | } |
---|
504 | |
---|
505 | G4int GetParticleIndex(G4ParticleDefinition * const particle) |
---|
506 | { |
---|
507 | G4int i(0); |
---|
508 | G4ParticleDefinition * const * particles(ParticleList()); |
---|
509 | |
---|
510 | while (*particles) |
---|
511 | { |
---|
512 | if ((*particles)==particle) |
---|
513 | return i; |
---|
514 | |
---|
515 | i++; |
---|
516 | particles++; |
---|
517 | } |
---|
518 | |
---|
519 | return -1; |
---|
520 | } |
---|
521 | |
---|
522 | G4ParticleDefinition * GetSelectedParticle(const struct Options & options) |
---|
523 | { |
---|
524 | G4ParticleDefinition * const * particles(ParticleList()); |
---|
525 | while (*particles) |
---|
526 | { |
---|
527 | if ((*particles)->GetParticleName()==options.particle) |
---|
528 | return (*particles); |
---|
529 | |
---|
530 | particles++; |
---|
531 | } |
---|
532 | |
---|
533 | G4cout << "Available particles are: " << G4endl; |
---|
534 | particles=ParticleList(); |
---|
535 | |
---|
536 | while (*particles) |
---|
537 | { |
---|
538 | G4cout << (*particles)->GetParticleName(); |
---|
539 | |
---|
540 | particles++; |
---|
541 | |
---|
542 | if (*particles) |
---|
543 | G4cout << ", "; |
---|
544 | } |
---|
545 | |
---|
546 | G4cout << G4endl; |
---|
547 | |
---|
548 | exit(-2); |
---|
549 | return 0; |
---|
550 | } |
---|
551 | |
---|
552 | //! \brief Setup processes |
---|
553 | //! \param options Options for the process choice |
---|
554 | void SetPhysics(const struct Options & options, G4VPhysicalVolume * world) |
---|
555 | { |
---|
556 | G4ProductionCutsTable * cutsTable(G4ProductionCutsTable::GetProductionCutsTable()); |
---|
557 | G4ProductionCuts * cuts(cutsTable->GetDefaultProductionCuts()); |
---|
558 | |
---|
559 | G4ParticleDefinition * const * particles(ParticleList()); |
---|
560 | while (*particles) |
---|
561 | { |
---|
562 | cuts->SetProductionCut(1*micrometer, *particles); |
---|
563 | particles++; |
---|
564 | } |
---|
565 | |
---|
566 | cutsTable->UpdateCoupleTable(world); |
---|
567 | G4cout << "[OK] Cuts are defined " << G4endl; |
---|
568 | |
---|
569 | G4ParticleDefinition * particle(GetSelectedParticle(options)); |
---|
570 | G4VProcess * eProcess=GetSelectedProcess(options); |
---|
571 | if (! (eProcess->IsApplicable(*particle))) |
---|
572 | { |
---|
573 | G4cout<< "Process " << eProcess->GetProcessName() << " is not applicable to " << particle->GetParticleName() << G4endl; |
---|
574 | exit(0); |
---|
575 | return; |
---|
576 | } |
---|
577 | |
---|
578 | G4cout<< "[OK] Process " << eProcess->GetProcessName() << " is applicable to " << particle->GetParticleName() << G4endl; |
---|
579 | |
---|
580 | G4ProcessManager * gProcessManager(new G4ProcessManager(particle)); |
---|
581 | particle->SetProcessManager(gProcessManager); |
---|
582 | gProcessManager->AddDiscreteProcess(eProcess); |
---|
583 | |
---|
584 | G4cout << "[OK] Processes are defined " << G4endl; |
---|
585 | |
---|
586 | |
---|
587 | G4cout << "[OK] Building physics tables" << G4endl; |
---|
588 | eProcess->BuildPhysicsTable(* particle); |
---|
589 | |
---|
590 | G4cout << "[OK] Physics tables built" << G4endl; |
---|
591 | } |
---|
592 | |
---|
593 | //! \brief Generates the step |
---|
594 | //! \param options Options related to the track generation |
---|
595 | //! \return The generated track |
---|
596 | G4Step * GenerateStep(const struct Options & options) |
---|
597 | { |
---|
598 | G4ThreeVector momentumDirection; |
---|
599 | |
---|
600 | momentumDirection.setRThetaPhi(1., std::acos(2*G4UniformRand()-1.), twopi * G4UniformRand()); |
---|
601 | |
---|
602 | G4double lnEnergyMin=std::log(options.minEnergy); |
---|
603 | G4double lnEnergyMax=std::log(options.maxEnergy); |
---|
604 | G4DynamicParticle * dynamicElectron(new G4DynamicParticle(GetSelectedParticle(options), momentumDirection, std::exp(lnEnergyMin+(lnEnergyMax-lnEnergyMin)*G4UniformRand()))); |
---|
605 | |
---|
606 | G4Track * aTrack(new G4Track(dynamicElectron, 0., G4ThreeVector(0., 0., 0.))); |
---|
607 | |
---|
608 | G4Step* aStep(new G4Step()); |
---|
609 | aStep->SetTrack(aTrack); |
---|
610 | aTrack->SetStep(aStep); |
---|
611 | |
---|
612 | G4Material * material(GetSelectedMaterial(options)); |
---|
613 | G4ProductionCutsTable * cutsTable(G4ProductionCutsTable::GetProductionCutsTable()); |
---|
614 | const G4MaterialCutsCouple * theCouple(cutsTable->GetMaterialCutsCouple(material, cutsTable->GetDefaultProductionCuts())); |
---|
615 | |
---|
616 | G4StepPoint * aPoint(new G4StepPoint()); |
---|
617 | aPoint->SetPosition(G4ThreeVector(0., 0., 0.)); |
---|
618 | aPoint->SetMaterial(material); |
---|
619 | aPoint->SetMaterialCutsCouple(theCouple); |
---|
620 | aPoint->SetSafety(10000.*cm); |
---|
621 | |
---|
622 | aStep->SetPreStepPoint(aPoint); |
---|
623 | |
---|
624 | return aStep; |
---|
625 | } |
---|
626 | |
---|
627 | void ProgressBar(G4int remainingIterations) |
---|
628 | { |
---|
629 | static time_t startingTime; |
---|
630 | static time_t nextDumpTime; |
---|
631 | static G4int startingIteration(0); |
---|
632 | time_t now; |
---|
633 | |
---|
634 | if (remainingIterations==0) |
---|
635 | { |
---|
636 | startingIteration=0; |
---|
637 | } |
---|
638 | else if (startingIteration==0) |
---|
639 | { |
---|
640 | startingTime=time(0); |
---|
641 | nextDumpTime=startingTime+3; |
---|
642 | startingIteration=remainingIterations; |
---|
643 | } |
---|
644 | else |
---|
645 | { |
---|
646 | now=time(0); |
---|
647 | if (now>nextDumpTime) |
---|
648 | { |
---|
649 | nextDumpTime=now+10; |
---|
650 | G4double time; |
---|
651 | G4double perc; |
---|
652 | |
---|
653 | time=std::floor(static_cast<G4double>(now-startingTime)/static_cast<G4double>(startingIteration-remainingIterations)*static_cast<G4double>(remainingIterations)+0.5); |
---|
654 | perc=std::floor(static_cast<G4double>(remainingIterations)/static_cast<G4double>(startingIteration)*200.+.5)/2.; |
---|
655 | |
---|
656 | G4cout << " " << perc << " % Remaining time: " << time << " s \r"; |
---|
657 | G4cout.flush(); |
---|
658 | } |
---|
659 | } |
---|
660 | } |
---|
661 | |
---|
662 | //! \brief Test the mean free path table |
---|
663 | //! \param tupleFactory The tuple factory |
---|
664 | //! \param options Options related to the mean free path test |
---|
665 | void MeanFreePathTest(AIDA::ITupleFactory * tupleFactory, const struct Options & options) |
---|
666 | { |
---|
667 | AIDA::ITuple* iTuple = tupleFactory->create("1", "Mean Free Path Ntuple", "double k, log_k, mfp, log_mfp, cpu_time"); |
---|
668 | |
---|
669 | G4double energy(options.minEnergy); |
---|
670 | G4double stpEnergy(std::pow(options.maxEnergy/energy, 1./static_cast<G4double>(options.nEnergySteps-1))); |
---|
671 | G4int step(options.nEnergySteps); |
---|
672 | |
---|
673 | G4ForceCondition condition; |
---|
674 | G4VLowEnergyTestableDiscreteProcess * process(GetSelectedProcess(options)); |
---|
675 | |
---|
676 | G4double mfp; |
---|
677 | clock_t time; |
---|
678 | |
---|
679 | ProgressBar(0); |
---|
680 | while (step>0) |
---|
681 | { |
---|
682 | G4Step * aStep(GenerateStep(options)); |
---|
683 | G4Track * aTrack(aStep->GetTrack()); |
---|
684 | |
---|
685 | if (!options.randomEnergy) |
---|
686 | { |
---|
687 | aTrack->SetKineticEnergy(energy); |
---|
688 | energy*=stpEnergy; |
---|
689 | } |
---|
690 | ProgressBar(step); |
---|
691 | step--; |
---|
692 | |
---|
693 | time=clock(); |
---|
694 | mfp=process->DumpMeanFreePath(*aTrack, 1.*mm, &condition)/m; |
---|
695 | time=clock()-time; |
---|
696 | |
---|
697 | iTuple->fill(iTuple->findColumn("k"), aTrack->GetKineticEnergy()/eV); |
---|
698 | iTuple->fill(iTuple->findColumn("log_k"), std::log10(aTrack->GetKineticEnergy()/eV)); |
---|
699 | iTuple->fill(iTuple->findColumn("mfp"), mfp); |
---|
700 | iTuple->fill(iTuple->findColumn("log_mfp"), std::log10(mfp)); |
---|
701 | iTuple->fill(iTuple->findColumn("cpu_time"), static_cast<G4double>(time)/static_cast<G4double>(CLOCKS_PER_SEC)); |
---|
702 | iTuple->addRow(); |
---|
703 | |
---|
704 | delete aTrack; |
---|
705 | delete aStep; |
---|
706 | } |
---|
707 | } |
---|
708 | |
---|
709 | //! \brief Test the post step do it |
---|
710 | //! \param tupleFactory The tuple factory |
---|
711 | //! \param options Options related to the post step do it test |
---|
712 | void PostStepDoItTest(AIDA::ITupleFactory * tupleFactory, const struct Options & options) |
---|
713 | { |
---|
714 | AIDA::ITuple* iTuple = tupleFactory->create("2", "Post Step Do It Test", "double iteration, step, in_k, log_in_k, in_theta, in_phi, in_pol_theta, in_pol_phi, e_deposit, log_e_deposit, trk_status, out_k, log_out_k, out_theta, out_phi, out_pol_theta, out_pol_phi, n_secondaries, sec0_type, sec0_k, log_sec0_k, sec0_theta, sec0_phi, sec0_pol_theta, sec0_pol_phi, sec1_type, sec1_k, log_sec1_k, sec1_theta, sec1_phi, sec1_pol_theta, sec1_pol_phi, sec2_type, sec2_k, log_sec2_k, sec2_theta, sec2_phi, sec2_pol_theta, sec2_pol_phi, sec3_type, sec3_k, log_sec3_k, sec3_theta, sec3_phi, sec3_pol_theta, sec3_pol_phi, cpu_time"); |
---|
715 | |
---|
716 | G4double energy(options.minEnergy); |
---|
717 | G4double stpEnergy(std::pow(options.maxEnergy/energy, 1./static_cast<G4double>(options.nEnergySteps-1))); |
---|
718 | G4int step(options.nEnergySteps); |
---|
719 | |
---|
720 | G4VLowEnergyTestableDiscreteProcess * process(GetSelectedProcess(options)); |
---|
721 | clock_t time; |
---|
722 | |
---|
723 | char strBuffer[20]; |
---|
724 | |
---|
725 | ProgressBar(0); |
---|
726 | while (step>0) |
---|
727 | { |
---|
728 | G4Step * aStep(GenerateStep(options)); |
---|
729 | G4Track * aTrack(aStep->GetTrack()); |
---|
730 | const G4DynamicParticle * aParticle(aTrack->GetDynamicParticle()); |
---|
731 | G4ThreeVector vector; |
---|
732 | |
---|
733 | if (!options.randomEnergy) |
---|
734 | { |
---|
735 | aTrack->SetKineticEnergy(energy); |
---|
736 | energy*=stpEnergy; |
---|
737 | } |
---|
738 | ProgressBar(step); |
---|
739 | step--; |
---|
740 | |
---|
741 | G4int iteration(options.nIterations); |
---|
742 | |
---|
743 | while (iteration>0) |
---|
744 | { |
---|
745 | iteration--; |
---|
746 | |
---|
747 | aStep->SetStepLength(1*micrometer); |
---|
748 | |
---|
749 | iTuple->fill(iTuple->findColumn("iteration"), iteration); |
---|
750 | iTuple->fill(iTuple->findColumn("step"), aStep->GetStepLength()/m); |
---|
751 | |
---|
752 | iTuple->fill(iTuple->findColumn("in_k"), aParticle->GetKineticEnergy()/eV); |
---|
753 | iTuple->fill(iTuple->findColumn("log_in_k"), std::log10(aParticle->GetKineticEnergy()/eV)); |
---|
754 | vector=aParticle->GetMomentumDirection(); |
---|
755 | iTuple->fill(iTuple->findColumn("in_theta"), vector.theta()); |
---|
756 | iTuple->fill(iTuple->findColumn("in_phi"), vector.phi()); |
---|
757 | vector=aParticle->GetPolarization(); |
---|
758 | iTuple->fill(iTuple->findColumn("in_pol_theta"), vector.theta()); |
---|
759 | iTuple->fill(iTuple->findColumn("in_pol_phi"), vector.phi()); |
---|
760 | |
---|
761 | time=clock(); |
---|
762 | G4ParticleChange * particleChange(dynamic_cast<G4ParticleChange *>(process->PostStepDoIt(*aTrack, *aStep))); |
---|
763 | time=clock()-time; |
---|
764 | |
---|
765 | aTrack->SetKineticEnergy(particleChange->GetEnergy()); |
---|
766 | aTrack->SetMomentumDirection(*particleChange->GetMomentumDirection()); |
---|
767 | aTrack->SetPolarization(*particleChange->GetPolarization()); |
---|
768 | |
---|
769 | iTuple->fill(iTuple->findColumn("e_deposit"), particleChange->GetLocalEnergyDeposit()/eV); |
---|
770 | iTuple->fill(iTuple->findColumn("log_e_deposit"), std::log10(particleChange->GetLocalEnergyDeposit()/eV)); |
---|
771 | iTuple->fill(iTuple->findColumn("trk_status"), particleChange->GetTrackStatus()); |
---|
772 | |
---|
773 | iTuple->fill(iTuple->findColumn("out_k"), aParticle->GetKineticEnergy()/eV); |
---|
774 | iTuple->fill(iTuple->findColumn("log_out_k"), std::log10(aParticle->GetKineticEnergy()/eV)); |
---|
775 | vector=aParticle->GetMomentumDirection(); |
---|
776 | iTuple->fill(iTuple->findColumn("out_theta"), vector.theta()); |
---|
777 | iTuple->fill(iTuple->findColumn("out_phi"), vector.phi()); |
---|
778 | vector=aParticle->GetPolarization(); |
---|
779 | iTuple->fill(iTuple->findColumn("out_pol_theta"), vector.theta()); |
---|
780 | iTuple->fill(iTuple->findColumn("out_pol_phi"), vector.phi()); |
---|
781 | |
---|
782 | G4int n(particleChange->GetNumberOfSecondaries()); |
---|
783 | iTuple->fill(iTuple->findColumn("n_secondaries"), n); |
---|
784 | |
---|
785 | while (n>0) |
---|
786 | { |
---|
787 | n--; |
---|
788 | G4Track * aSecTrack(particleChange->GetSecondary(n)); |
---|
789 | const G4DynamicParticle * aSecParticle(aSecTrack->GetDynamicParticle()); |
---|
790 | |
---|
791 | sprintf(strBuffer, "sec%d_type", n); |
---|
792 | iTuple->fill(iTuple->findColumn(strBuffer), GetParticleIndex(aSecParticle->GetDefinition())); |
---|
793 | sprintf(strBuffer, "sec%d_k", n); |
---|
794 | iTuple->fill(iTuple->findColumn(strBuffer), aSecParticle->GetKineticEnergy()/eV); |
---|
795 | sprintf(strBuffer, "log_sec%d_k", n); |
---|
796 | iTuple->fill(iTuple->findColumn(strBuffer), std::log10(aSecParticle->GetKineticEnergy()/eV)); |
---|
797 | vector=aSecParticle->GetMomentumDirection(); |
---|
798 | sprintf(strBuffer, "sec%d_theta", n); |
---|
799 | iTuple->fill(iTuple->findColumn(strBuffer), vector.theta()); |
---|
800 | sprintf(strBuffer, "sec%d_phi", n); |
---|
801 | iTuple->fill(iTuple->findColumn(strBuffer), vector.phi()); |
---|
802 | vector=aSecParticle->GetPolarization(); |
---|
803 | sprintf(strBuffer, "sec%d_pol_theta", n); |
---|
804 | iTuple->fill(iTuple->findColumn(strBuffer), vector.theta()); |
---|
805 | sprintf(strBuffer, "sec%d_pol_phi", n); |
---|
806 | iTuple->fill(iTuple->findColumn(strBuffer), vector.phi()); |
---|
807 | |
---|
808 | delete aSecTrack; |
---|
809 | } |
---|
810 | |
---|
811 | iTuple->fill(iTuple->findColumn("cpu_time"), static_cast<G4double>(time)/static_cast<G4double>(CLOCKS_PER_SEC)); |
---|
812 | iTuple->addRow(); |
---|
813 | |
---|
814 | particleChange->Clear(); |
---|
815 | } |
---|
816 | |
---|
817 | delete aTrack; |
---|
818 | delete aStep; |
---|
819 | } |
---|
820 | } |
---|
821 | |
---|
822 | //! \brief Main function |
---|
823 | //! \param argc Number of arguments |
---|
824 | //! \param argv Pointer to the arguments |
---|
825 | //! \return The exit value |
---|
826 | int main(int argc, char ** argv) |
---|
827 | { |
---|
828 | struct Options options; |
---|
829 | processOptions(argc, argv, &options); |
---|
830 | |
---|
831 | CreateMaterials(); |
---|
832 | |
---|
833 | GetSelectedProcess(options); |
---|
834 | GetSelectedMaterial(options); |
---|
835 | |
---|
836 | G4RunManager* rm = new G4RunManager(); |
---|
837 | rm->GeometryHasBeenModified(); |
---|
838 | G4VPhysicalVolume * world(CreateGeometry(options)); |
---|
839 | rm->DefineWorldVolume(world); |
---|
840 | G4cout << "[OK] World is defined " << G4endl; |
---|
841 | |
---|
842 | SetPhysics(options, world); |
---|
843 | |
---|
844 | if (!(options.meanFreePathTest || options.postStepDoItTest)) |
---|
845 | { |
---|
846 | G4cout << "[OK] Program completed" << G4endl; |
---|
847 | return 0; |
---|
848 | } |
---|
849 | |
---|
850 | // HBOOK initialization |
---|
851 | AIDA::IAnalysisFactory * analysisFactory(AIDA_createAnalysisFactory()); |
---|
852 | AIDA::ITreeFactory * treeFactory(analysisFactory->createTreeFactory()); |
---|
853 | AIDA::ITree * tree(treeFactory->create(options.outputFileName, "hbook", false, true)); |
---|
854 | G4cout << "[OK] Tree store: " << tree->storeName() << G4endl; |
---|
855 | |
---|
856 | AIDA::ITupleFactory * tupleFactory(analysisFactory->createTupleFactory(*tree)); |
---|
857 | |
---|
858 | // Mean free path test |
---|
859 | if (options.meanFreePathTest) |
---|
860 | { |
---|
861 | G4cout << "[OK] Mean free path test started" << G4endl; |
---|
862 | MeanFreePathTest(tupleFactory, options); |
---|
863 | G4cout << "[OK] Mean free path test completed" << G4endl; |
---|
864 | } |
---|
865 | |
---|
866 | // Post step do it test |
---|
867 | if (options.postStepDoItTest) |
---|
868 | { |
---|
869 | G4cout << "[OK] Post step do it test started" << G4endl; |
---|
870 | PostStepDoItTest(tupleFactory, options); |
---|
871 | G4cout << "[OK] Post step do it test completed" << G4endl; |
---|
872 | } |
---|
873 | |
---|
874 | G4cout << "[OK] Storing analysis data" << G4endl; |
---|
875 | tree->commit(); |
---|
876 | tree->close(); |
---|
877 | |
---|
878 | G4cout << "[OK] Deleting analysis data" << G4endl; |
---|
879 | delete tupleFactory; |
---|
880 | delete tree; |
---|
881 | delete treeFactory; |
---|
882 | delete analysisFactory; |
---|
883 | |
---|
884 | G4cout << "[OK] Program completed" << G4endl; |
---|
885 | return 0; |
---|
886 | } |
---|