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 11.1.2)


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