Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/eRosita/physics/src/G4RDGenerator2BN.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/eRosita/physics/src/G4RDGenerator2BN.cc (Version 11.3.0) and /examples/advanced/eRosita/physics/src/G4RDGenerator2BN.cc (Version 9.6)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27 // -------------------------------------------     27 // -------------------------------------------------------------------
 28 //                                                 28 //
 29 // GEANT4 Class file                               29 // GEANT4 Class file
 30 //                                                 30 //
 31 //                                                 31 //
 32 // File name:     G4RDGenerator2BN                 32 // File name:     G4RDGenerator2BN
 33 //                                                 33 //
 34 // Author:        Andreia Trindade (andreia@li     34 // Author:        Andreia Trindade (andreia@lip.pt)
 35 //                Pedro Rodrigues  (psilva@lip     35 //                Pedro Rodrigues  (psilva@lip.pt)
 36 //                Luis Peralta     (luis@lip.p     36 //                Luis Peralta     (luis@lip.pt)
 37 //                                                 37 //
 38 // Creation date: 21 June 2003                     38 // Creation date: 21 June 2003
 39 //                                                 39 //
 40 // Modifications:                                  40 // Modifications: 
 41 // 21 Jun 2003                                     41 // 21 Jun 2003                                 First implementation acording with new design
 42 // 05 Nov 2003  MGP                                42 // 05 Nov 2003  MGP                            Fixed compilation warning
 43 // 14 Mar 2004                                     43 // 14 Mar 2004                                 Code optimization
 44 //                                                 44 //
 45 // Class Description:                              45 // Class Description: 
 46 //                                                 46 //
 47 // Concrete base class for Bremsstrahlung Angu     47 // Concrete base class for Bremsstrahlung Angular Distribution Generation - 2BN Distribution
 48 //                                                 48 //
 49 // Class Description: End                          49 // Class Description: End 
 50 //                                                 50 //
 51 // -------------------------------------------     51 // -------------------------------------------------------------------
 52 //                                                 52 //
 53 //                                                 53 //    
 54                                                    54 
 55 #include "G4RDGenerator2BN.hh"                     55 #include "G4RDGenerator2BN.hh"
 56 #include "Randomize.hh"                            56 #include "Randomize.hh"
 57 #include "G4PhysicalConstants.hh"                  57 #include "G4PhysicalConstants.hh"
 58 #include "G4SystemOfUnits.hh"                      58 #include "G4SystemOfUnits.hh"
 59 //                                                 59 //    
 60                                                    60 
 61 G4double G4RDGenerator2BN::Atab[320] =             61 G4double G4RDGenerator2BN::Atab[320] =
 62          { 0.0264767, 0.0260006, 0.0257112, 0.     62          { 0.0264767, 0.0260006, 0.0257112, 0.0252475, 0.024792, 0.0243443, 0.0240726, 0.0236367,
 63            0.023209,  0.0227892, 0.0225362, 0.     63            0.023209,  0.0227892, 0.0225362, 0.0221284, 0.0217282,0.0214894, 0.0211005, 0.0207189,
 64            0.0204935, 0.0201227, 0.0197588, 0.     64            0.0204935, 0.0201227, 0.0197588, 0.019546,  0.0191923,0.0188454, 0.0186445, 0.0183072,
 65            0.0179763, 0.0177866, 0.0174649, 0.     65            0.0179763, 0.0177866, 0.0174649, 0.0172828, 0.0169702,0.0166634, 0.0164915, 0.0161933,
 66            0.0160283, 0.0157384, 0.0155801, 0.     66            0.0160283, 0.0157384, 0.0155801, 0.0152981, 0.0151463,0.0148721, 0.0147263, 0.0144598,
 67            0.01432,   0.0140607, 0.0139267, 0.     67            0.01432,   0.0140607, 0.0139267, 0.0136744, 0.0135459,0.0133005, 0.0131773, 0.0129385,
 68            0.0128205, 0.0125881, 0.012475,  0.     68            0.0128205, 0.0125881, 0.012475,  0.0122488, 0.0121406,0.0119204, 0.0118167, 0.0117158,
 69            0.0115032, 0.0114067, 0.0111995, 0.     69            0.0115032, 0.0114067, 0.0111995, 0.0111072, 0.0110175,0.0108173, 0.0107316, 0.0105365,
 70            0.0104547, 0.0102646, 0.0101865, 0.     70            0.0104547, 0.0102646, 0.0101865, 0.010111,  0.00992684,0.0098548,0.00967532,0.00960671,
 71            0.00943171,0.00936643,0.00930328,0.     71            0.00943171,0.00936643,0.00930328,0.0091337, 0.00907372,0.00890831,0.00885141,0.00869003,
 72            0.00863611,0.00858428,0.00842757,0.     72            0.00863611,0.00858428,0.00842757,0.00837854,0.0082256,0.00817931,0.00803,0.00798639,
 73            0.00784058,0.00779958,0.00776046,0.     73            0.00784058,0.00779958,0.00776046,0.00761866,0.00758201,0.00744346,0.00740928,0.00727384,
 74            0.00724201,0.00710969,0.00708004,0.     74            0.00724201,0.00710969,0.00708004,0.00695074,0.00692333,0.00679688,0.00677166,0.00664801,
 75            0.00662484,0.00650396,0.00648286,0.     75            0.00662484,0.00650396,0.00648286,0.00636458,0.00634545,0.00622977,0.00621258,0.00609936,
 76            0.00608412,0.00597331,0.00595991,0.     76            0.00608412,0.00597331,0.00595991,0.00585143,0.00583988,0.0057337,0.0057239,0.00561991,
 77            0.0056119, 0.00551005,0.00550377,0.     77            0.0056119, 0.00551005,0.00550377,0.00540399,0.00539938,0.00530162,0.00529872,0.00520292,
 78            0.0051091, 0.00510777,0.00501582,0.     78            0.0051091, 0.00510777,0.00501582,0.00501608,0.00492594,0.00492781,0.00483942,0.0048429,
 79            0.00475622,0.00476127,0.00467625,0.     79            0.00475622,0.00476127,0.00467625,0.00468287,0.00459947,0.00451785,0.00452581,0.00444573,
 80            0.00445522,0.00437664,0.00438768,0.     80            0.00445522,0.00437664,0.00438768,0.00431057,0.00432316,0.00424745,0.0042616,0.00418726,
 81            0.004203,  0.00413,   0.00405869,0.     81            0.004203,  0.00413,   0.00405869,0.00407563,0.00400561,0.00402414,0.00395536,0.00397553,
 82            0.00390795,0.00392975,0.00386339,0.     82            0.00390795,0.00392975,0.00386339,0.00379862,0.00382167,0.00375805,0.00378276,0.00372031,
 83            0.00374678,0.00368538,0.00371363,0.     83            0.00374678,0.00368538,0.00371363,0.00365335,0.00359463,0.0036242, 0.00356653,0.003598,
 84            0.00354139,0.00357481,0.00351921,0.     84            0.00354139,0.00357481,0.00351921,0.00355464,0.00350005,0.0034471,0.00348403,0.00343208,
 85            0.0034712, 0.00342026,0.00346165,0.     85            0.0034712, 0.00342026,0.00346165,0.00341172,0.00345548,0.00340657,0.00335944,0.00340491,
 86            0.00335885,0.00340692,0.00336191,0.     86            0.00335885,0.00340692,0.00336191,0.00341273,0.00336879,0.00342249,0.00337962,0.00333889,
 87            0.00339463,0.00335506,0.00341401,0.     87            0.00339463,0.00335506,0.00341401,0.00337558,0.00343797,0.00340067,0.00336584,0.00343059,
 88            0.0033969, 0.00346557,0.00343302,0.     88            0.0033969, 0.00346557,0.00343302,0.00350594,0.00347448,0.00344563,0.00352163,0.00349383,
 89            0.00357485,0.00354807,0.00352395,0.     89            0.00357485,0.00354807,0.00352395,0.00360885,0.00358571,0.00367661,0.00365446,0.00375194,
 90            0.00373078,0.00371234,0.00381532,0.     90            0.00373078,0.00371234,0.00381532,0.00379787,0.00390882,0.00389241,0.00387881,0.00399675,
 91            0.00398425,0.00411183,0.00410042,0.     91            0.00398425,0.00411183,0.00410042,0.00409197,0.00422843,0.00422123,0.00436974,0.0043637,
 92            0.00436082,0.00452075,0.00451934,0.     92            0.00436082,0.00452075,0.00451934,0.00452125,0.00469406,0.00469756,0.00488741,0.00489221,
 93            0.00490102,0.00510782,0.00511801,0.     93            0.00490102,0.00510782,0.00511801,0.00513271,0.0053589, 0.00537524,0.00562643,0.00564452,
 94            0.0056677, 0.00594482,0.00596999,0.     94            0.0056677, 0.00594482,0.00596999,0.0059999, 0.00630758,0.00634014,0.00637849,0.00672136,
 95            0.00676236,0.00680914,0.00719407,0.     95            0.00676236,0.00680914,0.00719407,0.0072439, 0.00730063,0.0077349, 0.00779494,0.00786293,
 96            0.00835577,0.0084276, 0.00850759,0.     96            0.00835577,0.0084276, 0.00850759,0.00907162,0.00915592,0.00924925,0.00935226,0.00999779,
 97            0.0101059, 0.0102249, 0.0109763, 0.     97            0.0101059, 0.0102249, 0.0109763, 0.0111003, 0.0112367, 0.0113862, 0.0122637, 0.0124196,
 98            0.0125898, 0.0136311, 0.0138081, 0.     98            0.0125898, 0.0136311, 0.0138081, 0.0140011, 0.0142112, 0.0154536, 0.0156723, 0.0159099,
 99            0.016168,  0.0176664, 0.0179339, 0.     99            0.016168,  0.0176664, 0.0179339, 0.0182246, 0.0185396, 0.020381,  0.0207026, 0.0210558,
100            0.0214374, 0.0237377, 0.0241275, 0.    100            0.0214374, 0.0237377, 0.0241275, 0.0245528, 0.0250106, 0.0255038, 0.0284158, 0.0289213,
101            0.0294621, 0.0300526, 0.0338619, 0.    101            0.0294621, 0.0300526, 0.0338619, 0.0344537, 0.0351108, 0.0358099, 0.036554,  0.0416399
102          };                                       102          };
103                                                   103 
104 G4double G4RDGenerator2BN::ctab[320] =            104 G4double G4RDGenerator2BN::ctab[320] =
105     { 0.482253, 0.482253, 0.489021, 0.489021,     105     { 0.482253, 0.482253, 0.489021, 0.489021, 0.489021, 0.489021, 0.495933,
106       0.495933, 0.495933, 0.495933, 0.502993,     106       0.495933, 0.495933, 0.495933, 0.502993, 0.502993, 0.502993, 0.510204,
107       0.510204, 0.510204, 0.517572, 0.517572,     107       0.510204, 0.510204, 0.517572, 0.517572, 0.517572, 0.5251,   0.5251,
108       0.5251,   0.532793, 0.532793, 0.532793,     108       0.5251,   0.532793, 0.532793, 0.532793, 0.540657, 0.540657, 0.548697,
109       0.548697, 0.548697, 0.556917, 0.556917,     109       0.548697, 0.548697, 0.556917, 0.556917, 0.565323, 0.565323, 0.573921,
110       0.573921, 0.582717, 0.582717, 0.591716,     110       0.573921, 0.582717, 0.582717, 0.591716, 0.591716, 0.600925, 0.600925,
111       0.610352, 0.610352, 0.620001, 0.620001,     111       0.610352, 0.610352, 0.620001, 0.620001, 0.629882, 0.629882, 0.64,
112       0.64,     0.650364, 0.650364, 0.660982,     112       0.64,     0.650364, 0.650364, 0.660982, 0.660982, 0.671862, 0.683013,
113       0.683013, 0.694444, 0.694444, 0.706165,     113       0.683013, 0.694444, 0.694444, 0.706165, 0.718184, 0.718184, 0.730514,
114       0.730514, 0.743163, 0.743163, 0.756144,     114       0.730514, 0.743163, 0.743163, 0.756144, 0.769468, 0.769468, 0.783147,
115       0.783147, 0.797194, 0.797194, 0.811622,     115       0.783147, 0.797194, 0.797194, 0.811622, 0.826446, 0.826446, 0.84168,
116       0.84168,  0.857339, 0.857339, 0.873439,     116       0.84168,  0.857339, 0.857339, 0.873439, 0.889996, 0.889996, 0.907029,
117       0.907029, 0.924556, 0.924556, 0.942596,     117       0.907029, 0.924556, 0.924556, 0.942596, 0.942596, 0.961169, 0.980296,
118       0.980296, 1,        1,        1.0203,       118       0.980296, 1,        1,        1.0203,   1.0203,   1.04123,  1.04123,
119       1.06281,  1.06281,  1.08507,  1.08507,      119       1.06281,  1.06281,  1.08507,  1.08507,  1.10803,  1.10803,  1.13173,
120       1.13173,  1.1562,   1.1562,   1.18147,      120       1.13173,  1.1562,   1.1562,   1.18147,  1.18147,  1.20758,  1.20758,
121       1.23457,  1.23457,  1.26247,  1.26247,      121       1.23457,  1.23457,  1.26247,  1.26247,  1.29132,  1.29132,  1.32118,
122       1.32118,  1.35208,  1.35208,  1.38408,      122       1.32118,  1.35208,  1.35208,  1.38408,  1.38408,  1.41723,  1.41723,
123       1.45159,  1.45159,  1.45159,  1.48721,      123       1.45159,  1.45159,  1.45159,  1.48721,  1.48721,  1.52416,  1.52416,
124       1.5625,   1.5625,   1.60231,  1.60231,      124       1.5625,   1.5625,   1.60231,  1.60231,  1.64366,  1.64366,  1.68663,
125       1.68663,  1.68663,  1.7313,   1.7313,       125       1.68663,  1.68663,  1.7313,   1.7313,   1.77778,  1.77778,  1.82615,
126       1.82615,  1.87652,  1.87652,  1.92901,      126       1.82615,  1.87652,  1.87652,  1.92901,  1.92901,  1.98373,  1.98373,
127       1.98373,  2.04082,  2.04082,  2.1004,       127       1.98373,  2.04082,  2.04082,  2.1004,   2.1004,   2.16263,  2.16263,
128       2.22767,  2.22767,  2.22767,  2.29568,      128       2.22767,  2.22767,  2.22767,  2.29568,  2.29568,  2.36686,  2.36686,
129       2.44141,  2.44141,  2.51953,  2.51953,      129       2.44141,  2.44141,  2.51953,  2.51953,  2.51953,  2.60146,  2.60146,
130       2.68745,  2.68745,  2.77778,  2.77778,      130       2.68745,  2.68745,  2.77778,  2.77778,  2.87274,  2.87274,  2.87274,
131       2.97265,  2.97265,  3.07787,  3.07787,      131       2.97265,  2.97265,  3.07787,  3.07787,  3.18878,  3.18878,  3.30579,
132       3.30579,  3.30579,  3.42936,  3.42936,      132       3.30579,  3.30579,  3.42936,  3.42936,  3.55999,  3.55999,  3.69822,
133       3.69822,  3.84468,  3.84468,  3.84468,      133       3.69822,  3.84468,  3.84468,  3.84468,  4,        4,        4.16493,
134       4.16493,  4.34028,  4.34028,  4.34028,      134       4.16493,  4.34028,  4.34028,  4.34028,  4.52694,  4.52694,  4.7259,
135       4.7259,   4.93827,  4.93827,  4.93827,      135       4.7259,   4.93827,  4.93827,  4.93827,  5.16529,  5.16529,  5.40833,
136       5.40833,  5.40833,  5.66893,  5.66893,      136       5.40833,  5.40833,  5.66893,  5.66893,  5.94884,  5.94884,  6.25,
137       6.25,     6.25,     6.57462,  6.57462,      137       6.25,     6.25,     6.57462,  6.57462,  6.92521,  6.92521,  6.92521,
138       7.3046,   7.3046,   7.71605,  7.71605,      138       7.3046,   7.3046,   7.71605,  7.71605,  7.71605,  8.16327,  8.16327,
139       8.65052,  8.65052,  8.65052,  9.18274,      139       8.65052,  8.65052,  8.65052,  9.18274,  9.18274,  9.18274,  9.76562,
140       9.76562,  10.4058,  10.4058,  10.4058,      140       9.76562,  10.4058,  10.4058,  10.4058,  11.1111,  11.1111,  11.1111,
141       11.8906,  11.8906,  12.7551,  12.7551,      141       11.8906,  11.8906,  12.7551,  12.7551,  12.7551,  13.7174,  13.7174,
142       13.7174,  14.7929,  14.7929,  14.7929,      142       13.7174,  14.7929,  14.7929,  14.7929,  16,       16,       16,
143       17.3611,  17.3611,  17.3611,  18.9036,      143       17.3611,  17.3611,  17.3611,  18.9036,  18.9036,  18.9036,  20.6612,
144       20.6612,  20.6612,  22.6757,  22.6757,      144       20.6612,  20.6612,  22.6757,  22.6757,  22.6757,  22.6757,  25,
145       25,       25,       27.7008,  27.7008,      145       25,       25,       27.7008,  27.7008,  27.7008,  27.7008,  30.8642,
146       30.8642,  30.8642,  34.6021,  34.6021,      146       30.8642,  30.8642,  34.6021,  34.6021,  34.6021,  34.6021,  39.0625,
147       39.0625,  39.0625,  39.0625,  44.4444,      147       39.0625,  39.0625,  39.0625,  44.4444,  44.4444,  44.4444,  44.4444,
148       51.0204,  51.0204,  51.0204,  51.0204,      148       51.0204,  51.0204,  51.0204,  51.0204,  59.1716,  59.1716,  59.1716,
149       59.1716,  59.1716,  69.4444,  69.4444,      149       59.1716,  59.1716,  69.4444,  69.4444,  69.4444,  69.4444,  82.6446,
150       82.6446,  82.6446,  82.6446,  82.6446,      150       82.6446,  82.6446,  82.6446,  82.6446,  100
151     };                                            151     };
152                                                   152 
153                                                   153 
154 G4RDGenerator2BN::G4RDGenerator2BN(const G4Str    154 G4RDGenerator2BN::G4RDGenerator2BN(const G4String& name):G4RDVBremAngularDistribution(name)
155 {                                                 155 {
156   b = 1.2;                                        156   b = 1.2;
157   index_min = -300;                               157   index_min = -300;
158   index_max = 20;                                 158   index_max = 20;
159                                                   159 
160   // Set parameters minimum limits Ekelectron     160   // Set parameters minimum limits Ekelectron = 250 eV and kphoton = 100 eV
161   kmin = 100*eV;                                  161   kmin = 100*eV;
162   Ekmin = 250*eV;                                 162   Ekmin = 250*eV;
163   kcut = 100*eV;                                  163   kcut = 100*eV;
164                                                   164 
165   // Increment Theta value for surface interpo    165   // Increment Theta value for surface interpolation
166   dtheta = 0.1*rad;                               166   dtheta = 0.1*rad;
167                                                   167 
168   // Construct Majorant Surface to 2BN cross-s    168   // Construct Majorant Surface to 2BN cross-section
169   // (we dont need this. Pre-calculated values    169   // (we dont need this. Pre-calculated values are used instead due to performance issues
170   // ConstructMajorantSurface();                  170   // ConstructMajorantSurface();
171 }                                                 171 }
172                                                   172 
173 //                                                173 //    
174                                                   174 
175 G4RDGenerator2BN::~G4RDGenerator2BN()             175 G4RDGenerator2BN::~G4RDGenerator2BN() 
176 {;}                                               176 {;}
177                                                   177 
178 //                                                178 //
179                                                   179 
180 G4double G4RDGenerator2BN::PolarAngle(const G4    180 G4double G4RDGenerator2BN::PolarAngle(const G4double initial_energy,
181             const G4double final_energy,          181             const G4double final_energy,
182             const G4int) // Z                     182             const G4int) // Z
183 {                                                 183 {
184                                                   184 
185   G4double theta = 0;                             185   G4double theta = 0;
186                                                   186 
187   G4double k = initial_energy - final_energy;     187   G4double k = initial_energy - final_energy;
188                                                   188 
189   theta = Generate2BN(initial_energy, k);         189   theta = Generate2BN(initial_energy, k);
190                                                   190 
191   return theta;                                   191   return theta;
192 }                                                 192 }
193 //                                                193 //
194                                                   194 
195 G4double G4RDGenerator2BN::CalculateFkt(G4doub    195 G4double G4RDGenerator2BN::CalculateFkt(G4double k, G4double theta, G4double A, G4double c) const
196 {                                                 196 {
197   G4double Fkt_value = 0;                         197   G4double Fkt_value = 0;
198   Fkt_value = A*std::pow(k,-b)*theta/(1+c*thet    198   Fkt_value = A*std::pow(k,-b)*theta/(1+c*theta*theta);
199   return Fkt_value;                               199   return Fkt_value;
200 }                                                 200 }
201                                                   201 
202 G4double G4RDGenerator2BN::Calculatedsdkdt(G4d    202 G4double G4RDGenerator2BN::Calculatedsdkdt(G4double kout, G4double theta, G4double Eel) const
203 {                                                 203 {
204                                                   204 
205   G4double dsdkdt_value = 0.;                     205   G4double dsdkdt_value = 0.;
206   G4double Z = 1;                                 206   G4double Z = 1;
207   // classic radius (in cm)                       207   // classic radius (in cm)
208   G4double r0 = 2.82E-13;                         208   G4double r0 = 2.82E-13;
209   //squared classic radius (in barn)              209   //squared classic radius (in barn)
210   G4double r02 = r0*r0*1.0E+24;                   210   G4double r02 = r0*r0*1.0E+24;
211                                                   211 
212   // Photon energy cannot be greater than elec    212   // Photon energy cannot be greater than electron kinetic energy
213   if(kout > (Eel-electron_mass_c2)){              213   if(kout > (Eel-electron_mass_c2)){
214     dsdkdt_value = 0;                             214     dsdkdt_value = 0;
215     return dsdkdt_value;                          215     return dsdkdt_value;
216   }                                               216   }
217                                                   217 
218      G4double E0 = Eel/electron_mass_c2;          218      G4double E0 = Eel/electron_mass_c2;
219      G4double k =  kout/electron_mass_c2;         219      G4double k =  kout/electron_mass_c2;
220      G4double E =  E0 - k;                        220      G4double E =  E0 - k;
221                                                   221 
222      if(E <= 1*MeV ){                             222      if(E <= 1*MeV ){                                 // Kinematic limit at 1 MeV !!
223        dsdkdt_value = 0;                          223        dsdkdt_value = 0;
224        return dsdkdt_value;                       224        return dsdkdt_value;
225      }                                            225      }
226                                                   226 
227                                                   227 
228      G4double p0 = std::sqrt(E0*E0-1);            228      G4double p0 = std::sqrt(E0*E0-1);
229      G4double p = std::sqrt(E*E-1);               229      G4double p = std::sqrt(E*E-1);
230      G4double L = std::log((E*E0-1+p*p0)/(E*E0    230      G4double L = std::log((E*E0-1+p*p0)/(E*E0-1-p*p0));
231      G4double delta0 = E0 - p0*std::cos(theta)    231      G4double delta0 = E0 - p0*std::cos(theta);
232      G4double epsilon = std::log((E+p)/(E-p));    232      G4double epsilon = std::log((E+p)/(E-p));
233      G4double Z2 = Z*Z;                           233      G4double Z2 = Z*Z;
234      G4double sintheta2 = std::sin(theta)*std:    234      G4double sintheta2 = std::sin(theta)*std::sin(theta);
235      G4double E02 = E0*E0;                        235      G4double E02 = E0*E0;
236      G4double E2 = E*E;                           236      G4double E2 = E*E;
237      G4double p02 = E0*E0-1;                      237      G4double p02 = E0*E0-1;
238      G4double k2 = k*k;                           238      G4double k2 = k*k;
239      G4double delta02 = delta0*delta0;            239      G4double delta02 = delta0*delta0;
240      G4double delta04 =  delta02* delta02;        240      G4double delta04 =  delta02* delta02;
241      G4double Q = std::sqrt(p02+k2-2*k*p0*std:    241      G4double Q = std::sqrt(p02+k2-2*k*p0*std::cos(theta));
242      G4double Q2 = Q*Q;                           242      G4double Q2 = Q*Q;
243      G4double epsilonQ = std::log((Q+p)/(Q-p))    243      G4double epsilonQ = std::log((Q+p)/(Q-p));
244                                                   244 
245                                                   245 
246      dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1    246      dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) * (p/p0) *
247        ( (8 * (sintheta2*(2*E02+1))/(p02*delta    247        ( (8 * (sintheta2*(2*E02+1))/(p02*delta04)) -
248          ((2*(5*E02+2*E*E0+3))/(p02 * delta02)    248          ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) -
249          ((2*(p02-k2))/((Q2*delta02))) +          249          ((2*(p02-k2))/((Q2*delta02))) +
250          ((4*E)/(p02*delta0)) +                   250          ((4*E)/(p02*delta0)) +
251          (L/(p*p0))*(                             251          (L/(p*p0))*(
252                  ((4*E0*sintheta2*(3*k-p02*E))    252                  ((4*E0*sintheta2*(3*k-p02*E))/(p02*delta04)) +
253                  ((4*E02*(E02+E2))/(p02*delta0    253                  ((4*E02*(E02+E2))/(p02*delta02)) +
254                  ((2-2*(7*E02-3*E*E0+E2))/(p02    254                  ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) +
255                  (2*k*(E02+E*E0-1))/((p02*delt    255                  (2*k*(E02+E*E0-1))/((p02*delta0))
256                  ) -                              256                  ) -
257          ((4*epsilon)/(p*delta0)) +               257          ((4*epsilon)/(p*delta0)) +
258          ((epsilonQ)/(p*Q))*                      258          ((epsilonQ)/(p*Q))*
259          (4/delta02-(6*k/delta0)-(2*k*(p02-k2)    259          (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/(Q2*delta0))
260          );                                       260          );
261                                                   261 
262                                                   262 
263                                                   263 
264      dsdkdt_value = dsdkdt_value*std::sin(thet    264      dsdkdt_value = dsdkdt_value*std::sin(theta);
265      return dsdkdt_value;                         265      return dsdkdt_value;
266 }                                                 266 }
267                                                   267 
268 void G4RDGenerator2BN::ConstructMajorantSurfac    268 void G4RDGenerator2BN::ConstructMajorantSurface()
269 {                                                 269 {
270                                                   270 
271   G4double Eel;                                   271   G4double Eel;
272   G4int vmax;                                     272   G4int vmax;
273   G4double Ek;                                    273   G4double Ek;
274   G4double k, theta, k0, theta0, thetamax;        274   G4double k, theta, k0, theta0, thetamax;
275   G4double ds, df, dsmax, dk, dt;                 275   G4double ds, df, dsmax, dk, dt;
276   G4double ratmin;                                276   G4double ratmin;
277   G4double ratio = 0;                             277   G4double ratio = 0;
278   G4double Vds, Vdf;                              278   G4double Vds, Vdf;
279   G4double A, c;                                  279   G4double A, c;
280                                                   280 
281   G4cout << "**** Constructing Majorant Surfac    281   G4cout << "**** Constructing Majorant Surface for 2BN Distribution ****" << G4endl;
282                                                   282 
283   if(kcut > kmin) kmin = kcut;                    283   if(kcut > kmin) kmin = kcut;
284                                                   284 
285   G4int i = 0;                                    285   G4int i = 0;
286   // Loop on energy index                         286   // Loop on energy index
287   for(G4int index = index_min; index < index_m    287   for(G4int index = index_min; index < index_max; index++){
288                                                   288 
289   G4double fraction = index/100.;                 289   G4double fraction = index/100.;
290   Ek = std::pow(10.,fraction);                    290   Ek = std::pow(10.,fraction);
291   Eel = Ek + electron_mass_c2;                    291   Eel = Ek + electron_mass_c2;
292                                                   292 
293   // find x-section maximum at k=kmin             293   // find x-section maximum at k=kmin
294   dsmax = 0.;                                     294   dsmax = 0.;
295   thetamax = 0.;                                  295   thetamax = 0.;
296                                                   296 
297   for(theta = 0.; theta < pi; theta = theta +     297   for(theta = 0.; theta < pi; theta = theta + dtheta){
298                                                   298 
299     ds = Calculatedsdkdt(kmin, theta, Eel);       299     ds = Calculatedsdkdt(kmin, theta, Eel);
300                                                   300 
301     if(ds > dsmax){                               301     if(ds > dsmax){
302       dsmax = ds;                                 302       dsmax = ds;
303       thetamax = theta;                           303       thetamax = theta;
304     }                                             304     }
305   }                                               305   }
306                                                   306 
307   // Compute surface paremeters at kmin           307   // Compute surface paremeters at kmin
308   if(Ek < kmin || thetamax == 0){                 308   if(Ek < kmin || thetamax == 0){
309     c = 0;                                        309     c = 0;
310     A = 0;                                        310     A = 0;
311   }else{                                          311   }else{
312     c = 1/(thetamax*thetamax);                    312     c = 1/(thetamax*thetamax);
313     A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b    313     A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b));
314   }                                               314   }
315                                                   315 
316   // look for correction factor to normalizati    316   // look for correction factor to normalization at kmin 
317   ratmin = 1.;                                    317   ratmin = 1.;
318                                                   318 
319   // Volume under surfaces                        319   // Volume under surfaces
320   Vds = 0.;                                       320   Vds = 0.;
321   Vdf = 0.;                                       321   Vdf = 0.;
322   k0 = 0.;                                        322   k0 = 0.;
323   theta0 = 0.;                                    323   theta0 = 0.;
324                                                   324 
325   vmax = G4int(100.*std::log10(Ek/kmin));         325   vmax = G4int(100.*std::log10(Ek/kmin));
326                                                   326 
327   for(G4int v = 0; v < vmax; v++){                327   for(G4int v = 0; v < vmax; v++){
328     G4double fraction = (v/100.);                 328     G4double fraction = (v/100.);
329     k = std::pow(10.,fraction)*kmin;              329     k = std::pow(10.,fraction)*kmin;
330                                                   330 
331     for(theta = 0.; theta < pi; theta = theta     331     for(theta = 0.; theta < pi; theta = theta + dtheta){
332       dk = k - k0;                                332       dk = k - k0;
333       dt = theta - theta0;                        333       dt = theta - theta0;
334       ds = Calculatedsdkdt(k,theta, Eel);         334       ds = Calculatedsdkdt(k,theta, Eel);
335       Vds = Vds + ds*dk*dt;                       335       Vds = Vds + ds*dk*dt;
336       df = CalculateFkt(k,theta, A, c);           336       df = CalculateFkt(k,theta, A, c);
337       Vdf = Vdf + df*dk*dt;                       337       Vdf = Vdf + df*dk*dt;
338                                                   338 
339       if(ds != 0.){                               339       if(ds != 0.){
340   if(df != 0.) ratio = df/ds;                     340   if(df != 0.) ratio = df/ds;
341       }                                           341       }
342                                                   342 
343       if(ratio < ratmin && ratio != 0.){          343       if(ratio < ratmin && ratio != 0.){
344   ratmin = ratio;                                 344   ratmin = ratio;
345       }                                           345       }
346     }                                             346     }
347   }                                               347   }
348                                                   348 
349                                                   349 
350   // sampling efficiency                          350   // sampling efficiency
351   Atab[i] = A/ratmin * 1.04;                      351   Atab[i] = A/ratmin * 1.04;
352   ctab[i] = c;                                    352   ctab[i] = c;
353                                                   353 
354 //  G4cout << Ek << " " << i << " " << index <    354 //  G4cout << Ek << " " << i << " " << index << " " << Atab[i] << " " << ctab[i] << " " << G4endl;
355   i++;                                            355   i++;
356   }                                               356   }
357                                                   357 
358 }                                                 358 }
359                                                   359 
360 G4double G4RDGenerator2BN::Generate2BN(G4doubl    360 G4double G4RDGenerator2BN::Generate2BN(G4double Ek, G4double k) const 
361 {                                                 361 {
362                                                   362 
363   G4double Eel;                                   363   G4double Eel;
364   G4double t;                                     364   G4double t;
365   G4double cte2;                                  365   G4double cte2;
366   G4double y, u;                                  366   G4double y, u;
367   G4double fk, ft;                                367   G4double fk, ft;
368   G4double ds;                                    368   G4double ds;
369   G4double A2;                                    369   G4double A2;
370   G4double A, c;                                  370   G4double A, c;
371                                                   371 
372   G4int trials = 0;                               372   G4int trials = 0;
373   G4int index;                                    373   G4int index;
374                                                   374 
375   // find table index                             375   // find table index
376   index = G4int(std::log10(Ek)*100) - index_mi    376   index = G4int(std::log10(Ek)*100) - index_min;
377   Eel = Ek + electron_mass_c2;                    377   Eel = Ek + electron_mass_c2;
378                                                   378 
379   c = ctab[index];                                379   c = ctab[index];
380   A = Atab[index];                                380   A = Atab[index];
381   if(index < index_max){                          381   if(index < index_max){
382     A2 = Atab[index+1];                           382     A2 = Atab[index+1];
383     if(A2 > A) A = A2;                            383     if(A2 > A) A = A2;
384   }                                               384   }
385                                                   385 
386   do{                                             386   do{
387   // generate k accordimg to std::pow(k,-b)       387   // generate k accordimg to std::pow(k,-b)
388   trials++;                                       388   trials++;
389                                                   389 
390   // normalization constant                       390   // normalization constant 
391 //  cte1 = (1-b)/(std::pow(kmax,(1-b))-std::po    391 //  cte1 = (1-b)/(std::pow(kmax,(1-b))-std::pow(kmin2,(1-b)));
392 //  y = G4UniformRand();                          392 //  y = G4UniformRand();
393 //  k = std::pow(((1-b)/cte1*y+std::pow(kmin2,    393 //  k = std::pow(((1-b)/cte1*y+std::pow(kmin2,(1-b))),(1/(1-b)));
394                                                   394 
395   // generate theta accordimg to theta/(1+c*st    395   // generate theta accordimg to theta/(1+c*std::pow(theta,2)
396   // Normalization constant                       396   // Normalization constant
397   cte2 = 2*c/std::log(1+c*pi2);                   397   cte2 = 2*c/std::log(1+c*pi2);
398                                                   398 
399   y = G4UniformRand();                            399   y = G4UniformRand();
400   t = std::sqrt((std::exp(2*c*y/cte2)-1)/c);      400   t = std::sqrt((std::exp(2*c*y/cte2)-1)/c);
401   u = G4UniformRand();                            401   u = G4UniformRand();
402                                                   402 
403   // point acceptance algorithm                   403   // point acceptance algorithm
404   fk = std::pow(k,-b);                            404   fk = std::pow(k,-b);
405   ft = t/(1+c*t*t);                               405   ft = t/(1+c*t*t);
406   ds = Calculatedsdkdt(k,t,Eel);                  406   ds = Calculatedsdkdt(k,t,Eel);
407                                                   407 
408   // violation point                              408   // violation point
409   if(ds > (A*fk*ft)) G4cout << "WARNING IN G4R    409   if(ds > (A*fk*ft)) G4cout << "WARNING IN G4RDGenerator2BN !!!" << Ek << " " <<  (ds-A*fk*ft)/ds << G4endl;
410                                                   410 
411   }while(u*A*fk*ft > ds);                         411   }while(u*A*fk*ft > ds);
412                                                   412 
413   return t;                                       413   return t;
414                                                   414 
415 }                                                 415 }
416                                                   416 
417 void G4RDGenerator2BN::PrintGeneratorInformati    417 void G4RDGenerator2BN::PrintGeneratorInformation() const
418 {                                                 418 {
419   G4cout << "\n" << G4endl;                       419   G4cout << "\n" << G4endl;
420   G4cout << "Bremsstrahlung Angular Generator     420   G4cout << "Bremsstrahlung Angular Generator is 2BN Generator from 2BN Koch & Motz distribution (Rev Mod Phys 31(4), 920 (1959))" << G4endl;
421   G4cout << "\n" << G4endl;                       421   G4cout << "\n" << G4endl;
422 }                                                 422 } 
423                                                   423 
424                                                   424