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