Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4LivermoreBremsstrahlungModel.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 /processes/electromagnetic/lowenergy/src/G4LivermoreBremsstrahlungModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4LivermoreBremsstrahlungModel.cc (Version 10.5)


  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:     G4LivermoreBremsstrahlungMod     32 // File name:     G4LivermoreBremsstrahlungModel
 33 //                                                 33 //
 34 // Author:        Vladimir Ivanchenko use inhe     34 // Author:        Vladimir Ivanchenko use inheritance from Andreas Schaelicke
 35 //                base class implementing ultr     35 //                base class implementing ultra relativistic bremsstrahlung
 36 //                model                            36 //                model
 37 //                                                 37 //
 38 // Creation date: 04.10.2011                       38 // Creation date: 04.10.2011
 39 //                                                 39 //
 40 // Modifications:                                  40 // Modifications:
 41 //                                                 41 //
 42 // -------------------------------------------     42 // -------------------------------------------------------------------
 43 //                                                 43 //
 44 //....oooOO0OOooo........oooOO0OOooo........oo     44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 45 //....oooOO0OOooo........oooOO0OOooo........oo     45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 46                                                    46 
 47 #include "G4LivermoreBremsstrahlungModel.hh"       47 #include "G4LivermoreBremsstrahlungModel.hh"
 48 #include "G4PhysicalConstants.hh"                  48 #include "G4PhysicalConstants.hh"
 49 #include "G4SystemOfUnits.hh"                      49 #include "G4SystemOfUnits.hh"
 50 #include "G4Electron.hh"                           50 #include "G4Electron.hh"
 51 #include "G4Positron.hh"                           51 #include "G4Positron.hh"
 52 #include "G4Gamma.hh"                              52 #include "G4Gamma.hh"
 53 #include "Randomize.hh"                            53 #include "Randomize.hh"
 54 #include "G4AutoLock.hh"                       << 
 55 #include "G4Material.hh"                           54 #include "G4Material.hh"
 56 #include "G4Element.hh"                            55 #include "G4Element.hh"
 57 #include "G4ElementVector.hh"                      56 #include "G4ElementVector.hh"
 58 #include "G4ProductionCutsTable.hh"                57 #include "G4ProductionCutsTable.hh"
 59 #include "G4ParticleChangeForLoss.hh"              58 #include "G4ParticleChangeForLoss.hh"
 60 #include "G4Generator2BS.hh"                       59 #include "G4Generator2BS.hh"
 61                                                    60 
 62 #include "G4Physics2DVector.hh"                    61 #include "G4Physics2DVector.hh"
 63 #include "G4Exp.hh"                                62 #include "G4Exp.hh"
 64 #include "G4Log.hh"                                63 #include "G4Log.hh"
 65                                                    64 
 66 #include "G4ios.hh"                                65 #include "G4ios.hh"
 67 #include <fstream>                                 66 #include <fstream>
 68 #include <iomanip>                                 67 #include <iomanip>
 69                                                    68 
 70 //....oooOO0OOooo........oooOO0OOooo........oo     69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 71                                                    70 
 72 namespace { G4Mutex LivermoreBremsstrahlungMod << 
 73 using namespace std;                               71 using namespace std;
 74                                                    72 
 75                                                <<  73 G4Physics2DVector* G4LivermoreBremsstrahlungModel::dataSB[] = {0};
 76 G4Physics2DVector* G4LivermoreBremsstrahlungMo << 
 77 G4double G4LivermoreBremsstrahlungModel::ylimi     74 G4double G4LivermoreBremsstrahlungModel::ylimit[] = {0.0};
 78 G4double G4LivermoreBremsstrahlungModel::expnu     75 G4double G4LivermoreBremsstrahlungModel::expnumlim = -12.;
 79                                                    76 
 80 static const G4double emaxlog = 4*G4Log(10.);      77 static const G4double emaxlog = 4*G4Log(10.);
 81 static const G4double alpha = CLHEP::twopi*CLH     78 static const G4double alpha = CLHEP::twopi*CLHEP::fine_structure_const;
 82 static const G4double epeaklimit= 300*CLHEP::M     79 static const G4double epeaklimit= 300*CLHEP::MeV;
 83 static const G4double elowlimit = 20*CLHEP::ke     80 static const G4double elowlimit = 20*CLHEP::keV;
 84                                                    81 
 85 G4LivermoreBremsstrahlungModel::G4LivermoreBre     82 G4LivermoreBremsstrahlungModel::G4LivermoreBremsstrahlungModel(
 86   const G4ParticleDefinition* p, const G4Strin     83   const G4ParticleDefinition* p, const G4String& nam)
 87   : G4eBremsstrahlungRelModel(p,nam),useBicubi     84   : G4eBremsstrahlungRelModel(p,nam),useBicubicInterpolation(false)
 88 {                                                  85 {
 89   SetLowEnergyLimit(10.0*eV);                      86   SetLowEnergyLimit(10.0*eV);
                                                   >>  87   SetLPMFlag(false);
                                                   >>  88   nwarn = 0;
                                                   >>  89   idx = idy = 0;
 90   SetAngularDistribution(new G4Generator2BS())     90   SetAngularDistribution(new G4Generator2BS());
 91 }                                                  91 }
 92                                                    92 
 93 //....oooOO0OOooo........oooOO0OOooo........oo     93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 94                                                    94 
 95 G4LivermoreBremsstrahlungModel::~G4LivermoreBr     95 G4LivermoreBremsstrahlungModel::~G4LivermoreBremsstrahlungModel()
 96 {                                                  96 {
 97   if(IsMaster()) {                                 97   if(IsMaster()) {
 98     for(size_t i=0; i<101; ++i) {                  98     for(size_t i=0; i<101; ++i) {
 99       if(dataSB[i]) {                              99       if(dataSB[i]) {
100   delete dataSB[i];                               100   delete dataSB[i];
101   dataSB[i] = nullptr;                         << 101   dataSB[i] = 0;
102       }                                           102       }
103     }                                             103     }
104   }                                               104   }
105 }                                                 105 }
106                                                   106 
107 //....oooOO0OOooo........oooOO0OOooo........oo    107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108                                                   108 
109 void G4LivermoreBremsstrahlungModel::Initialis    109 void G4LivermoreBremsstrahlungModel::Initialise(const G4ParticleDefinition* p,
110             const G4DataVector& cuts)             110             const G4DataVector& cuts)
111 {                                                 111 {
112   // Access to elements                           112   // Access to elements
113   if(IsMaster()) {                                113   if(IsMaster()) {
                                                   >> 114 
114     // check environment variable                 115     // check environment variable
115     // Build the complete string identifying t    116     // Build the complete string identifying the file with the data set
116     const char* path = G4FindDataDir("G4LEDATA << 117     char* path = getenv("G4LEDATA");
117                                                   118 
118     const G4ElementTable* theElmTable = G4Elem    119     const G4ElementTable* theElmTable = G4Element::GetElementTable();
119     size_t numOfElm = G4Element::GetNumberOfEl    120     size_t numOfElm = G4Element::GetNumberOfElements();
120     if(numOfElm > 0) {                            121     if(numOfElm > 0) {
121       for(size_t i=0; i<numOfElm; ++i) {          122       for(size_t i=0; i<numOfElm; ++i) {
122   G4int Z = (*theElmTable)[i]->GetZasInt();    << 123   G4int Z = G4int(((*theElmTable)[i])->GetZ());
123   if(Z < 1)        { Z = 1; }                     124   if(Z < 1)        { Z = 1; }
124   else if(Z > 100) { Z = 100; }                   125   else if(Z > 100) { Z = 100; }
125   //G4cout << "Z= " << Z << G4endl;               126   //G4cout << "Z= " << Z << G4endl;
126   // Initialisation                               127   // Initialisation
127   if(!dataSB[Z]) { ReadData(Z, path); }           128   if(!dataSB[Z]) { ReadData(Z, path); }
128       }                                           129       }
129     }                                             130     }
130   }                                               131   }
                                                   >> 132 
131   G4eBremsstrahlungRelModel::Initialise(p, cut    133   G4eBremsstrahlungRelModel::Initialise(p, cuts);
132 }                                                 134 }
133                                                   135 
134 //....oooOO0OOooo........oooOO0OOooo........oo    136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135                                                   137 
136 G4String G4LivermoreBremsstrahlungModel::Direc    138 G4String G4LivermoreBremsstrahlungModel::DirectoryPath() const
137 {                                                 139 {
138   return "/livermore/brem/br";                    140   return "/livermore/brem/br";
139 }                                                 141 }
140                                                   142 
141 //....oooOO0OOooo........oooOO0OOooo........oo    143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142                                                   144 
143 void G4LivermoreBremsstrahlungModel::ReadData(    145 void G4LivermoreBremsstrahlungModel::ReadData(G4int Z, const char* path)
144 {                                                 146 {
                                                   >> 147   //  G4cout << "ReadData Z= " << Z << G4endl;
                                                   >> 148   // G4cout << "Status for Z= " << dataSB[Z] << G4endl;
                                                   >> 149   //if(path) { G4cout << path << G4endl; }
145   if(dataSB[Z]) { return; }                       150   if(dataSB[Z]) { return; }
146   const char* datadir = path;                     151   const char* datadir = path;
147                                                   152 
148   if(nullptr == datadir) {                     << 153   if(!datadir) {
149     datadir = G4FindDataDir("G4LEDATA");       << 154     datadir = getenv("G4LEDATA");
150     if(!datadir) {                                155     if(!datadir) {
151       G4Exception("G4LivermoreBremsstrahlungMo    156       G4Exception("G4LivermoreBremsstrahlungModel::ReadData()","em0006",
152       FatalException,"Environment variable G4L    157       FatalException,"Environment variable G4LEDATA not defined");
153       return;                                     158       return;
154     }                                             159     }
155   }                                               160   }
156   std::ostringstream ost;                         161   std::ostringstream ost;
157   ost << datadir << DirectoryPath() << Z;         162   ost << datadir << DirectoryPath() << Z;
158   std::ifstream fin(ost.str().c_str());           163   std::ifstream fin(ost.str().c_str());
159   if( !fin.is_open()) {                           164   if( !fin.is_open()) {
160     G4ExceptionDescription ed;                    165     G4ExceptionDescription ed;
161     ed << "Bremsstrahlung data file <" << ost.    166     ed << "Bremsstrahlung data file <" << ost.str().c_str()
162        << "> is not opened!";                     167        << "> is not opened!";
163     G4Exception("G4LivermoreBremsstrahlungMode    168     G4Exception("G4LivermoreBremsstrahlungModel::ReadData()","em0003",
164     FatalException,ed,                            169     FatalException,ed,
165     "G4LEDATA version should be G4EMLOW8.0 or  << 170     "G4LEDATA version should be G4EMLOW6.23 or later.");
166     return;                                       171     return;
167   }                                               172   }
168   //G4cout << "G4LivermoreBremsstrahlungModel     173   //G4cout << "G4LivermoreBremsstrahlungModel read from <" << ost.str().c_str()
169   //   << ">" << G4endl;                          174   //   << ">" << G4endl;
170   G4Physics2DVector* v = new G4Physics2DVector    175   G4Physics2DVector* v = new G4Physics2DVector();
171   if(v->Retrieve(fin)) {                          176   if(v->Retrieve(fin)) {
172     if(useBicubicInterpolation) { v->SetBicubi    177     if(useBicubicInterpolation) { v->SetBicubicInterpolation(true); }
173     dataSB[Z] = v;                                178     dataSB[Z] = v;
174     ylimit[Z] = v->Value(0.97, emaxlog, idx, i    179     ylimit[Z] = v->Value(0.97, emaxlog, idx, idy);
175   } else {                                        180   } else {
176     G4ExceptionDescription ed;                    181     G4ExceptionDescription ed;
177     ed << "Bremsstrahlung data file <" << ost.    182     ed << "Bremsstrahlung data file <" << ost.str().c_str()
178        << "> is not retrieved!";                  183        << "> is not retrieved!";
179     G4Exception("G4LivermoreBremsstrahlungMode    184     G4Exception("G4LivermoreBremsstrahlungModel::ReadData()","em0005",
180                 FatalException,ed,                185                 FatalException,ed,
181     "G4LEDATA version should be G4EMLOW8.0 or  << 186     "G4LEDATA version should be G4EMLOW6.23 or later.");
182     delete v;                                     187     delete v;
183   }                                               188   }
                                                   >> 189   // G4cout << dataSB[Z] << G4endl;
184 }                                                 190 }
185                                                   191 
186 //....oooOO0OOooo........oooOO0OOooo........oo    192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
187                                                   193 
188 G4double                                          194 G4double
189 G4LivermoreBremsstrahlungModel::ComputeDXSecti    195 G4LivermoreBremsstrahlungModel::ComputeDXSectionPerAtom(G4double gammaEnergy)
190 {                                                 196 {
                                                   >> 197 
191   if(gammaEnergy < 0.0 || fPrimaryKinEnergy <=    198   if(gammaEnergy < 0.0 || fPrimaryKinEnergy <= 0.0) { return 0.0; }
192   G4double x = gammaEnergy/fPrimaryKinEnergy;     199   G4double x = gammaEnergy/fPrimaryKinEnergy;
193   G4double y = G4Log(fPrimaryKinEnergy/MeV);      200   G4double y = G4Log(fPrimaryKinEnergy/MeV);
194   G4int Z = fCurrentIZ;                           201   G4int Z = fCurrentIZ;
195                                                   202 
196   //G4cout << "G4LivermoreBremsstrahlungModel:    203   //G4cout << "G4LivermoreBremsstrahlungModel::ComputeDXSectionPerAtom Z= " << Z
197   //   << " x= " << x << " y= " << y << " " <<    204   //   << " x= " << x << " y= " << y << " " << dataSB[Z] << G4endl;
198   if(!dataSB[Z]) { InitialiseForElement(0, Z);    205   if(!dataSB[Z]) { InitialiseForElement(0, Z); }
199                                                << 206   /*
                                                   >> 207     G4ExceptionDescription ed;
                                                   >> 208     ed << "Bremsstrahlung data for Z= " << Z
                                                   >> 209        << " are not initialized!";
                                                   >> 210     G4Exception("G4LivermoreBremsstrahlungModel::ComputeDXSectionPerAtom()",
                                                   >> 211                 "em0005", FatalException, ed,
                                                   >> 212     "G4LEDATA version should be G4EMLOW6.23 or later.");
                                                   >> 213   }
                                                   >> 214   */
200   G4double invb2 = fPrimaryTotalEnergy*fPrimar    215   G4double invb2 = fPrimaryTotalEnergy*fPrimaryTotalEnergy/(fPrimaryKinEnergy
201                    *(fPrimaryKinEnergy + 2.*fP    216                    *(fPrimaryKinEnergy + 2.*fPrimaryParticleMass));
202   G4double cross = dataSB[Z]->Value(x,y,idx,id    217   G4double cross = dataSB[Z]->Value(x,y,idx,idy)*invb2*millibarn/gBremFactor;
203                                                   218 
204   if(!fIsElectron) {                              219   if(!fIsElectron) {
205     G4double invbeta1 = sqrt(invb2);              220     G4double invbeta1 = sqrt(invb2);
206     G4double e2 = fPrimaryKinEnergy - gammaEne    221     G4double e2 = fPrimaryKinEnergy - gammaEnergy;
207     if(e2 > 0.0) {                                222     if(e2 > 0.0) {
208       G4double invbeta2 = (e2 + fPrimaryPartic    223       G4double invbeta2 = (e2 + fPrimaryParticleMass)
209                           /sqrt(e2*(e2 + 2.*fP    224                           /sqrt(e2*(e2 + 2.*fPrimaryParticleMass));
210       G4double xxx = alpha*fCurrentIZ*(invbeta    225       G4double xxx = alpha*fCurrentIZ*(invbeta1 - invbeta2);
211       if(xxx < expnumlim) { cross = 0.0; }        226       if(xxx < expnumlim) { cross = 0.0; }
212       else { cross *= G4Exp(xxx); }               227       else { cross *= G4Exp(xxx); }
213     } else {                                      228     } else {
214       cross = 0.0;                                229       cross = 0.0;
215     }                                             230     }
216   }                                               231   }
217                                                   232 
218   return cross;                                   233   return cross;
219 }                                                 234 }
220                                                   235 
221 //....oooOO0OOooo........oooOO0OOooo........oo    236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222                                                   237 
223 void                                              238 void
224 G4LivermoreBremsstrahlungModel::SampleSecondar    239 G4LivermoreBremsstrahlungModel::SampleSecondaries(
225                                         std::v    240                                         std::vector<G4DynamicParticle*>* vdp,
226           const G4MaterialCutsCouple* couple,     241           const G4MaterialCutsCouple* couple,
227           const G4DynamicParticle* dp,            242           const G4DynamicParticle* dp,
228           G4double cutEnergy,                     243           G4double cutEnergy,
229           G4double maxEnergy)                     244           G4double maxEnergy)
230 {                                                 245 {
231   G4double kineticEnergy = dp->GetKineticEnerg    246   G4double kineticEnergy = dp->GetKineticEnergy();
232   G4double cut  = std::min(cutEnergy, kineticE    247   G4double cut  = std::min(cutEnergy, kineticEnergy);
233   G4double emax = std::min(maxEnergy, kineticE    248   G4double emax = std::min(maxEnergy, kineticEnergy);
234   if(cut >= emax) { return; }                     249   if(cut >= emax) { return; }
235   // sets total energy, kinetic energy and den    250   // sets total energy, kinetic energy and density correction
236   SetupForMaterial(fPrimaryParticle, couple->G    251   SetupForMaterial(fPrimaryParticle, couple->GetMaterial(), kineticEnergy);
237                                                   252 
238   const G4Element* elm =                          253   const G4Element* elm =
239     SelectRandomAtom(couple,fPrimaryParticle,k    254     SelectRandomAtom(couple,fPrimaryParticle,kineticEnergy,cut,emax);
240   fCurrentIZ = elm->GetZasInt();                  255   fCurrentIZ = elm->GetZasInt();
241   G4int Z = fCurrentIZ;                           256   G4int Z = fCurrentIZ;
242                                                   257 
243   G4double totMomentum = sqrt(kineticEnergy*(f    258   G4double totMomentum = sqrt(kineticEnergy*(fPrimaryTotalEnergy+electron_mass_c2));
244   /*                                              259   /*
245   G4cout << "G4LivermoreBremsstrahlungModel::S    260   G4cout << "G4LivermoreBremsstrahlungModel::SampleSecondaries E(MeV)= "
246    << kineticEnergy/MeV                           261    << kineticEnergy/MeV
247    << " Z= " << Z << " cut(MeV)= " << cut/MeV     262    << " Z= " << Z << " cut(MeV)= " << cut/MeV
248    << " emax(MeV)= " << emax/MeV << " corr= "     263    << " emax(MeV)= " << emax/MeV << " corr= " << fDensityCorr << G4endl;
249   */                                              264   */
250   G4double xmin = G4Log(cut*cut + fDensityCorr    265   G4double xmin = G4Log(cut*cut + fDensityCorr);
251   G4double xmax = G4Log(emax*emax  + fDensityC    266   G4double xmax = G4Log(emax*emax  + fDensityCorr);
252   G4double y = G4Log(kineticEnergy/MeV);          267   G4double y = G4Log(kineticEnergy/MeV);
253                                                   268 
254   G4double gammaEnergy, v;                        269   G4double gammaEnergy, v;
255                                                   270 
256   // majoranta                                    271   // majoranta
257   G4double x0 = cut/kineticEnergy;                272   G4double x0 = cut/kineticEnergy;
258   G4double vmax = dataSB[Z]->Value(x0, y, idx,    273   G4double vmax = dataSB[Z]->Value(x0, y, idx, idy)*1.02;
                                                   >> 274   //  G4double invbeta1 = 0;
259                                                   275 
260   // majoranta corrected for e-                   276   // majoranta corrected for e-
261   if(fIsElectron && x0 < 0.97 &&                  277   if(fIsElectron && x0 < 0.97 &&
262      ((kineticEnergy > epeaklimit) || (kinetic    278      ((kineticEnergy > epeaklimit) || (kineticEnergy < elowlimit))) {
263     G4double ylim = std::min(ylimit[Z],1.1*dat    279     G4double ylim = std::min(ylimit[Z],1.1*dataSB[Z]->Value(0.97,y,idx,idy));
264     if(ylim > vmax) { vmax = ylim; }              280     if(ylim > vmax) { vmax = ylim; }
265   }                                               281   }
266   if(x0 < 0.05) { vmax *= 1.2; }                  282   if(x0 < 0.05) { vmax *= 1.2; }
267                                                   283 
                                                   >> 284   //G4cout<<"y= "<<y<<" xmin= "<<xmin<<" xmax= "<<xmax
                                                   >> 285   //<<" vmax= "<<vmax<<G4endl;
                                                   >> 286   //  G4int ncount = 0;
268   do {                                            287   do {
269     //++ncount;                                   288     //++ncount;
270     G4double x = G4Exp(xmin + G4UniformRand()*    289     G4double x = G4Exp(xmin + G4UniformRand()*(xmax - xmin)) - fDensityCorr;
271     if(x < 0.0) { x = 0.0; }                      290     if(x < 0.0) { x = 0.0; }
272     gammaEnergy = sqrt(x);                        291     gammaEnergy = sqrt(x);
273     G4double x1 = gammaEnergy/kineticEnergy;      292     G4double x1 = gammaEnergy/kineticEnergy;
274     v = dataSB[Z]->Value(x1, y, idx, idy);        293     v = dataSB[Z]->Value(x1, y, idx, idy);
275                                                   294 
276     // correction for positrons                   295     // correction for positrons
277     if(!fIsElectron) {                            296     if(!fIsElectron) {
278       G4double e1 = kineticEnergy - cut;          297       G4double e1 = kineticEnergy - cut;
279       G4double invbeta1 = (e1 + fPrimaryPartic    298       G4double invbeta1 = (e1 + fPrimaryParticleMass)
280                           /sqrt(e1*(e1 + 2*fPr    299                           /sqrt(e1*(e1 + 2*fPrimaryParticleMass));
281       G4double e2 = kineticEnergy - gammaEnerg    300       G4double e2 = kineticEnergy - gammaEnergy;
282       G4double invbeta2 = (e2 + fPrimaryPartic    301       G4double invbeta2 = (e2 + fPrimaryParticleMass)
283                           /sqrt(e2*(e2 + 2*fPr    302                           /sqrt(e2*(e2 + 2*fPrimaryParticleMass));
284       G4double xxx = twopi*fine_structure_cons    303       G4double xxx = twopi*fine_structure_const*fCurrentIZ*(invbeta1 - invbeta2);
285                                                   304 
286       if(xxx < expnumlim) { v = 0.0; }            305       if(xxx < expnumlim) { v = 0.0; }
287       else { v *= G4Exp(xxx); }                   306       else { v *= G4Exp(xxx); }
288     }                                             307     }
289                                                   308 
290     if (v > 1.05*vmax && nwarn < 5) {             309     if (v > 1.05*vmax && nwarn < 5) {
291       ++nwarn;                                    310       ++nwarn;
292       G4ExceptionDescription ed;                  311       G4ExceptionDescription ed;
293       ed << "### G4LivermoreBremsstrahlungMode    312       ed << "### G4LivermoreBremsstrahlungModel Warning: Majoranta exceeded! "
294    << v << " > " << vmax << " by " << v/vmax      313    << v << " > " << vmax << " by " << v/vmax
295    << " Egamma(MeV)= " << gammaEnergy             314    << " Egamma(MeV)= " << gammaEnergy
296    << " Ee(MeV)= " << kineticEnergy               315    << " Ee(MeV)= " << kineticEnergy
297    << " Z= " << Z << "  " << fPrimaryParticle-    316    << " Z= " << Z << "  " << fPrimaryParticle->GetParticleName();
298                                                   317 
299       if ( 20 == nwarn ) {                        318       if ( 20 == nwarn ) {
300   ed << "\n ### G4LivermoreBremsstrahlungModel    319   ed << "\n ### G4LivermoreBremsstrahlungModel Warnings stopped";
301       }                                           320       }
302       G4Exception("G4LivermoreBremsstrahlungMo    321       G4Exception("G4LivermoreBremsstrahlungModel::SampleScattering","em0044",
303       JustWarning, ed,"");                        322       JustWarning, ed,"");
304                                                   323 
305     }                                             324     }
306   } while (v < vmax*G4UniformRand());             325   } while (v < vmax*G4UniformRand());
307                                                   326 
308   //                                              327   //
309   // angles of the emitted gamma. ( Z - axis a    328   // angles of the emitted gamma. ( Z - axis along the parent particle)
310   // use general interface                        329   // use general interface
311   //                                              330   //
312                                                   331 
313   G4ThreeVector gammaDirection =                  332   G4ThreeVector gammaDirection =
314     GetAngularDistribution()->SampleDirection(    333     GetAngularDistribution()->SampleDirection(dp,fPrimaryTotalEnergy-gammaEnergy,
315                 Z, couple->GetMaterial());        334                 Z, couple->GetMaterial());
316                                                   335 
317   // create G4DynamicParticle object for the G    336   // create G4DynamicParticle object for the Gamma
318   G4DynamicParticle* gamma =                      337   G4DynamicParticle* gamma =
319     new G4DynamicParticle(fGammaParticle,gamma    338     new G4DynamicParticle(fGammaParticle,gammaDirection,gammaEnergy);
320   vdp->push_back(gamma);                          339   vdp->push_back(gamma);
321                                                   340 
322   G4ThreeVector direction = (totMomentum*dp->G    341   G4ThreeVector direction = (totMomentum*dp->GetMomentumDirection()
323            - gammaEnergy*gammaDirection).unit(    342            - gammaEnergy*gammaDirection).unit();
324                                                   343 
325   /*                                              344   /*
326   G4cout << "### G4SBModel: v= "                  345   G4cout << "### G4SBModel: v= "
327    << " Eg(MeV)= " << gammaEnergy                 346    << " Eg(MeV)= " << gammaEnergy
328    << " Ee(MeV)= " << kineticEnergy               347    << " Ee(MeV)= " << kineticEnergy
329    << " DirE " << direction << " DirG " << gam    348    << " DirE " << direction << " DirG " << gammaDirection
330    << G4endl;                                     349    << G4endl;
331   */                                              350   */
332   // energy of primary                            351   // energy of primary
333   G4double finalE = kineticEnergy - gammaEnerg    352   G4double finalE = kineticEnergy - gammaEnergy;
334                                                   353 
335   // stop tracking and create new secondary in    354   // stop tracking and create new secondary instead of primary
336   if(gammaEnergy > SecondaryThreshold()) {        355   if(gammaEnergy > SecondaryThreshold()) {
337     fParticleChange->ProposeTrackStatus(fStopA    356     fParticleChange->ProposeTrackStatus(fStopAndKill);
338     fParticleChange->SetProposedKineticEnergy(    357     fParticleChange->SetProposedKineticEnergy(0.0);
339     G4DynamicParticle* el =                       358     G4DynamicParticle* el =
340       new G4DynamicParticle(const_cast<G4Parti    359       new G4DynamicParticle(const_cast<G4ParticleDefinition*>(fPrimaryParticle),
341           direction, finalE);                     360           direction, finalE);
342     vdp->push_back(el);                           361     vdp->push_back(el);
343                                                   362 
344     // continue tracking                          363     // continue tracking
345   } else {                                        364   } else {
346     fParticleChange->SetProposedMomentumDirect    365     fParticleChange->SetProposedMomentumDirection(direction);
347     fParticleChange->SetProposedKineticEnergy(    366     fParticleChange->SetProposedKineticEnergy(finalE);
348   }                                               367   }
349 }                                                 368 }
350                                                   369 
351 //....oooOO0OOooo........oooOO0OOooo........oo    370 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
352                                                   371 
                                                   >> 372 #include "G4AutoLock.hh"
                                                   >> 373 namespace { G4Mutex LivermoreBremsstrahlungModelMutex = G4MUTEX_INITIALIZER; }
353 void G4LivermoreBremsstrahlungModel::Initialis    374 void G4LivermoreBremsstrahlungModel::InitialiseForElement(
354                                      const G4P    375                                      const G4ParticleDefinition*,
355              G4int Z)                             376              G4int Z)
356 {                                                 377 {
357   G4AutoLock l(&LivermoreBremsstrahlungModelMu    378   G4AutoLock l(&LivermoreBremsstrahlungModelMutex);
                                                   >> 379   //G4cout << "G4LivermoreBremsstrahlungModel::InitialiseForElement Z= "
                                                   >> 380   //<< Z << G4endl;
358   if(!dataSB[Z]) { ReadData(Z); }                 381   if(!dataSB[Z]) { ReadData(Z); }
359   l.unlock();                                     382   l.unlock();
360 }                                                 383 }
361                                                   384 
362 //....oooOO0OOooo........oooOO0OOooo........oo    385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 386 
                                                   >> 387 
363                                                   388