Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (14 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

Location:
trunk/source/processes/electromagnetic/utils
Files:
23 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/electromagnetic/utils/History

    r1228 r1315  
    1 $Id: History,v 1.400 2009/11/22 19:48:30 vnivanch Exp $
     1$Id: History,v 1.424 2010/06/04 15:36:39 vnivanch Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     204 June 10: V.Ivant (emutils-V09-03-16)
     21Previous tag was set from wrong directory
     22
     234 June 10: V.Ivant (emutils-V09-03-15)
     24G4EmConfigurator - fixed case of more than one model is added per a process
     25
     264 June 10: V.Ivant (emutils-V09-03-14)
     27G4EmCorrections - moved G4AtomicShell header into source
     28
     2926 May 10: V.Ivant (emutils-V09-03-13)
     30G4VEmModel - added method ChargeSquareRatio to access current charge of an ion
     31G4VEnergyLossProcess - use this new method
     32
     3310 May 10: V.Ivant (emutils-V09-03-12)
     34G4VEmProcess - cleanup printout at initialisation for scattering process
     35
     3628 April 10: V.Ivant (emutils-V09-03-11)
     37G4VEmProcess, G4VEnergyLossProcess, G4VEmModel - provided GetCurrentElement method
     38                       (addressed bug report #1115 and HyperNews request)
     39
     4027 April 10: V.Ivant (emutils-V09-03-10)
     41G4LossTableManager   - added class member and a method GetNumberOfBinsPerDecade
     42G4VEmModel           - use GetNumberOfBinsPerDecade and spline flag to initialise
     43                       G4EmElementSelector (addressed bug report #1115)
     44G4EmElementSelector  - use spline flag to construct vectors probabilities
     45G4EmProcessOptions   - removed double implementation of initialisation code,
     46                       which already exist in G4LossTableManager
     47G4VEnergyLossProcess - call CorrectionsAlongStep only for ions (minor CPU saving)
     48
     4923 April 10: V.Ivant (emutils-V09-03-09)
     50G4VEnergyLossProcess - removed unused variable
     51
     5212 April 10: V.Ivant (emutils-V09-03-08)
     53G4EmModelManager - do not use min energy cut defined by models allowing
     54                   decreasing of cuts in limit to zero
     55G4EmCalculator - fixed GetCrossSection method
     56
     5712 April 10: V.Ivant (emutils-V09-03-07)
     58G4LossTableManager - added methods PreparePhsyicsTables, BuildPhysicsTables,
     59                     and changed initialisation of models via G4EmConfigurator
     60G4VEnergyLossProcess, G4VEmProcess, G4VMultipleScattering - added
     61                     calls of new G4LossTableManager methods
     62                     PreparePhsyicsTables, BuildPhysicsTables
     63G4EmConfigarator - upgraded and fixed old problem
     64
     6506 April 10: V.Ivant (emutils-V09-03-06)
     66G4VEnergyLossProcess - use the same method to build cross section table
     67                       as DEDX table (use copy constructors to reduce
     68                       number of calls to std::exp)
     69G4EmModelManager - cleanup comments
     70
     7122 March 10: V.Ivant (emutils-V09-03-05)
     72G4EmCorrections - added protection against large Barkas and Bloch
     73                  corrections in the case of large negatively charged
     74                  particle (100*e-) - fixed problem reported by ATLAS
     75G4EmCalculator - cleanup
     76
     7710 March 10: V.Ivant (emutils-V09-03-04)
     78G4VEmModel, G4VEmProcess, G4VEnergyLossProcess, G4VMultipleScattering,
     79G4LossTableManager - reorder inline methods and add comments
     80
     8105 March 10: V.Ivant (emutils-V09-03-03)
     82G4VMscModel, G4VMultipleScattering - set skin=1.0 as a default
     83
     8424 February 10: V.Ivant (emutils-V09-03-02)
     85G4VEmProcess - move SetBuildTableFlag method from protected to public
     86
     8717 February 10: V.Ivant (emutils-V09-03-01)
     88G4VEmProcess - fixed problem for ion processes by adding pointer to
     89               currentParticle which may be different from GenericIon
     90
     9122 January 10: V.Ivant (emutils-V09-03-00)
     92G4VEmProcess - added protection against negative cross section
     93G4VEnergyLossProcess - added protection against negative cross section;
     94                     - improved logic in RetrieveTable method
    1995
    209623 November 09: V.Ivant (emutils-V09-02-24)

  • trunk/source/processes/electromagnetic/utils/include/G4EmCalculator.hh

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmCalculator.hh,v 1.19 2009/11/11 23:59:48 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4EmCalculator.hh,v 1.20 2010/04/13 10:58:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929//
    … …  
    249249  const G4Material*            currentMaterial;
    250250  const G4ParticleDefinition*  currentParticle;
     251  const G4ParticleDefinition*  lambdaParticle;
    251252  const G4ParticleDefinition*  baseParticle;
    252253  const G4PhysicsTable*        currentLambda;
    … …  
    260261
    261262  G4String                     currentName;
     263  G4String                     lambdaName;
    262264  G4double                     currentCut;
    263265  G4double                     chargeSquare;

  • trunk/source/processes/electromagnetic/utils/include/G4EmConfigurator.hh

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmConfigurator.hh,v 1.2 2008/11/21 12:30:29 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4EmConfigurator.hh,v 1.3 2010/04/12 11:44:40 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5858//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    5959
     60class G4VEnergyLossProcess;
     61class G4VEmProcess;
     62class G4VMultipleScattering;
     63
    6064class G4EmConfigurator
    6165{
    6266public:
    6367 
    64   G4EmConfigurator();
     68  G4EmConfigurator(G4int verboseLevel = 1);
    6569 
    6670  ~G4EmConfigurator();
    6771
    68   // Add EM model to the list of extra models potentially to be
    69   // declared for the G4Region and energy interval
    70   //
    71   void AddExtraEmModel(const G4String& particleName,
    72                        G4VEmModel*, G4VEmFluctuationModel* fm = 0);
    73 
    74   // Declare EM model for particle type and process to
    75   // be active for the G4Region and energy interval
    76   // The model should be previously added to the configurator
    77   // or be "dummy"
    78   //
    79   void AddModelForRegion(const G4String& particleName,
    80                          const G4String& processName,
    81                          const G4String& modelName,
    82                          const G4String& regionName = "",
    83                          G4double emin = 0.0,
    84                          G4double emax = DBL_MAX,
    85                          const G4String& flucModelName = "");
    86 
    8772  // Set EM model for particle type and process to
    8873  // be active for the G4Region and energy interval
     74  // The model will be added to the list
    8975  //
    9076  void SetExtraEmModel(const G4String& particleName,
    … …  
    9783
    9884  // Add all previously declared models to corresponding processes
     85  // Can be called in ConstructPhysics
     86  //
    9987  void AddModels();
     88
     89  // These methods called by G4LossTableManager
     90  //
     91  void PrepareModels(const G4ParticleDefinition* aParticle,
     92                     G4VEnergyLossProcess* p);
     93
     94  void PrepareModels(const G4ParticleDefinition* aParticle,
     95                     G4VEmProcess* p);
     96
     97  void PrepareModels(const G4ParticleDefinition* aParticle,
     98                     G4VMultipleScattering* p);
     99
     100  void Clear();
     101
     102  inline void SetVerbose(G4int value);
    100103
    101104private:
    102105
    103   void SetModelForRegion(const G4String& particleName,
     106  G4Region* FindRegion(const G4String&);
     107
     108  void SetModelForRegion(G4VEmModel* model,
     109                         G4VEmFluctuationModel* fm,
     110                         G4Region* reg,
     111                         const G4String& particleName,
    104112                         const G4String& processName,
    105                          const G4String& modelName,
    106                          const G4String& regionName,
    107                          const G4String& flucModelName,
    108113                         G4double emin,
    109114                         G4double emax);
     115
     116  G4bool UpdateModelEnergyRange(G4VEmModel* mod,
     117                                G4double emin, G4double emax);
    110118
    111119  // hide assignment operator
    … …  
    113121  G4EmConfigurator(const G4EmConfigurator&);
    114122
     123  std::vector<G4VEmModel*> models; 
     124  std::vector<G4VEmFluctuationModel*> flucModels; 
    115125  std::vector<G4String> particles; 
    116126  std::vector<G4String> processes; 
    117   std::vector<G4String> models; 
    118127  std::vector<G4String> regions; 
    119   std::vector<G4String> flucModels; 
    120128  std::vector<G4double> lowEnergy;
    121129  std::vector<G4double> highEnergy;
    122130 
    123   std::vector<G4String> particleList; 
    124   std::vector<G4VEmModel*> modelList;
    125   std::vector<G4VEmFluctuationModel*> flucModelList;
    126 
    127131  G4int index;
     132  G4int verbose;
    128133};
    129134
    130135//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     136
     137inline void G4EmConfigurator::SetVerbose(G4int value)
     138{
     139  verbose = value;
     140}
    131141
    132142#endif

  • trunk/source/processes/electromagnetic/utils/include/G4EmCorrections.hh

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmCorrections.hh,v 1.24 2008/09/12 14:44:48 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4EmCorrections.hh,v 1.25 2010/06/04 09:28:46 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5656
    5757#include "globals.hh"
    58 #include "G4AtomicShells.hh"
    5958#include "G4ionEffectiveCharge.hh"
    6059#include "G4Material.hh"
    … …  
    263262  G4int     nbinCorr;
    264263
    265   G4AtomicShells        shells;
    266264  G4ionEffectiveCharge  effCharge;
    267265
    … …  
    356354    tmax  = 2.0*electron_mass_c2*bg2 /(1. + 2.0*gamma*ratio + ratio*ratio);
    357355    charge  = p->GetPDGCharge()/eplus;
    358     if(charge < 1.5)  {q2 = charge*charge;}
    359     else {
    360       q2 = effCharge.EffectiveChargeSquareRatio(p,mat,kinEnergy);
    361       charge = std::sqrt(q2);
    362     }
     356    //if(charge < 1.5)  {q2 = charge*charge;}
     357    //else {
     358    //  q2 = effCharge.EffectiveChargeSquareRatio(p,mat,kinEnergy);
     359    //  charge = std::sqrt(q2);
     360    //}
     361    if(charge > 1.5) { charge = effCharge.EffectiveCharge(p,mat,kinEnergy); }
     362    q2 = charge*charge;
    363363  }
    364364  if(mat != material) {

  • trunk/source/processes/electromagnetic/utils/include/G4LossTableManager.hh

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LossTableManager.hh,v 1.55 2009/10/29 19:25:28 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4LossTableManager.hh,v 1.58 2010/04/27 16:59:52 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929//
    … …  
    6161// 12-02-07 Add SetSkin, SetLinearLossLimit (V.Ivanchenko)
    6262// 18-06-07 Move definition of msc parameters to G4EmProcessOptions (V.Ivanchenko)
     63// 12-04-10 Added PreparePhsyicsTables and BuildPhysicsTables entries (V.Ivanchenko)
    6364//
    6465// Class Description:
    … …  
    9394class G4EmConfigurator;
    9495class G4LossTableBuilder;
     96class G4Region;
    9597
    9698class G4LossTableManager
    … …  
    103105  ~G4LossTableManager();
    104106
     107  //-------------------------------------------------
     108  // called from destructor
     109  //-------------------------------------------------
     110
    105111  void Clear();
    106112
    107   // get the DEDX or the range for a given particle/energy/material
     113  //-------------------------------------------------
     114  // initialisation before a new run
     115  //-------------------------------------------------
     116
     117  void PreparePhysicsTable(const G4ParticleDefinition* aParticle,
     118                           G4VEnergyLossProcess* p);
     119  void PreparePhysicsTable(const G4ParticleDefinition* aParticle,
     120                           G4VEmProcess* p);
     121  void PreparePhysicsTable(const G4ParticleDefinition* aParticle,
     122                           G4VMultipleScattering* p);
     123  void BuildPhysicsTable(const G4ParticleDefinition* aParticle);
     124  void BuildPhysicsTable(const G4ParticleDefinition* aParticle,
     125                         G4VEnergyLossProcess* p);
     126
     127  //-------------------------------------------------
     128  // Run time access to DEDX, range, energy for a given particle,
     129  // energy, and G4MaterialCutsCouple
     130  //-------------------------------------------------
     131
    108132  inline G4double GetDEDX(
    109133    const G4ParticleDefinition *aParticle,
    … …  
    141165          G4double& length);
    142166
    143   // to be called only by energy loss processes
     167  //-------------------------------------------------
     168  // Methods to be called only at initialisation
     169  //-------------------------------------------------
     170
    144171  void Register(G4VEnergyLossProcess* p);
    145172
    … …  
    162189  void DeRegister(G4VEmFluctuationModel* p);
    163190
    164   void EnergyLossProcessIsInitialised(const G4ParticleDefinition* aParticle,
    165                                       G4VEnergyLossProcess* p);
    166  
    167191  void RegisterIon(const G4ParticleDefinition* aParticle,
    168192                   G4VEnergyLossProcess* p);
    … …  
    171195                             G4VEnergyLossProcess* p);
    172196
    173   void BuildPhysicsTable(const G4ParticleDefinition* aParticle,
    174                          G4VEnergyLossProcess* p);
    175 
    176197  void SetLossFluctuations(G4bool val);
    177198
    178   void SetSubCutoff(G4bool val);
     199  void SetSubCutoff(G4bool val, const G4Region* r=0);
    179200
    180201  void SetIntegral(G4bool val);
    … …  
    198219  void SetLambdaBinning(G4int val);
    199220
     221  G4int GetNumberOfBinsPerDecade() const;
     222
    200223  void SetStepFunction(G4double v1, G4double v2);
    201224
    … …  
    214237  void SetVerbose(G4int val);
    215238
     239  //-------------------------------------------------
     240  // Access methods
     241  //-------------------------------------------------
     242
    216243  G4EnergyLossMessenger* GetMessenger();
    217244
    … …  
    241268
    242269private:
     270
     271  //-------------------------------------------------
     272  // Private methods and members
     273  //-------------------------------------------------
    243274
    244275  G4LossTableManager();
    … …  
    251282  void ParticleHaveNoLoss(const G4ParticleDefinition* aParticle);
    252283
    253   void SetParameters(G4VEnergyLossProcess*);
     284  void SetParameters(const G4ParticleDefinition* aParticle,
     285                     G4VEnergyLossProcess*);
    254286
    255287  void CopyDEDXTables();
    256288
    257 private:
    258 
    259289  static G4LossTableManager* theInstance;
    260290
    261291  typedef const G4ParticleDefinition* PD;
     292
    262293  std::map<PD,G4VEnergyLossProcess*,std::less<PD> > loss_map;
    263294
    … …  
    279310  PD                    currentParticle;
    280311  PD                    theElectron;
     312  PD                    firstParticle;
    281313
    282314  G4int n_loss;
    … …  
    284316
    285317  G4bool all_tables_are_built;
    286   //  G4bool first_entry;
     318  G4bool startInitialisation;
     319
    287320  G4bool lossFluctuationFlag;
    288321  G4bool subCutoffFlag;
    … …  
    290323  G4bool integral;
    291324  G4bool integralActive;
    292   G4bool all_tables_are_stored;
    293325  G4bool buildCSDARange;
    294326  G4bool minEnergyActive;
    … …  
    314346  G4EmConfigurator*           emConfigurator;
    315347
    316   const G4ParticleDefinition* firstParticle;
     348  G4int nbinsLambda;
     349  G4int nbinsPerDecade;
    317350  G4int verbose;
    318351
    … …  
    322355//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    323356
     357inline G4VEnergyLossProcess* G4LossTableManager::GetEnergyLossProcess(
     358          const G4ParticleDefinition *aParticle)
     359{
     360  if(aParticle != currentParticle) {
     361    currentParticle = aParticle;
     362    std::map<PD,G4VEnergyLossProcess*,std::less<PD> >::const_iterator pos;
     363    if ((pos = loss_map.find(aParticle)) != loss_map.end()) {
     364      currentLoss = (*pos).second;
     365    } else {
     366      currentLoss = 0;
     367     // ParticleHaveNoLoss(aParticle);
     368    }
     369  }
     370  return currentLoss;
     371}
     372
     373//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     374
    324375inline G4double G4LossTableManager::GetDEDX(
    325376          const G4ParticleDefinition *aParticle,
    … …  
    327378          const G4MaterialCutsCouple *couple)
    328379{
    329   if(aParticle != currentParticle) GetEnergyLossProcess(aParticle);
     380  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
    330381  G4double x;
    331   if(currentLoss) x = currentLoss->GetDEDX(kineticEnergy, couple);
    332   else            x = G4EnergyLossTables::GetDEDX(
    333                       currentParticle,kineticEnergy,couple,false);
     382  if(currentLoss) { x = currentLoss->GetDEDX(kineticEnergy, couple); }
     383  else            { x = G4EnergyLossTables::GetDEDX(currentParticle,
     384                                                    kineticEnergy,couple,false); }
    334385  return x;
    335386}
    … …  
    342393          const G4MaterialCutsCouple *couple)
    343394{
    344   if(aParticle != currentParticle) GetEnergyLossProcess(aParticle);
     395  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
    345396  G4double x = 0.0;
    346   if(currentLoss) x = currentLoss->GetDEDXForSubsec(kineticEnergy, couple);
     397  if(currentLoss) { x = currentLoss->GetDEDXForSubsec(kineticEnergy, couple); }
    347398  return x;
    348399}
    … …  
    355406          const G4MaterialCutsCouple *couple)
    356407{
    357   if(aParticle != currentParticle) GetEnergyLossProcess(aParticle);
     408  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
    358409  G4double x = DBL_MAX;
    359   if(currentLoss) x = currentLoss->GetCSDARange(kineticEnergy, couple);
     410  if(currentLoss) { x = currentLoss->GetCSDARange(kineticEnergy, couple); }
    360411  return x;
    361412}
    … …  
    368419          const G4MaterialCutsCouple *couple)
    369420{
    370   if(aParticle != currentParticle) GetEnergyLossProcess(aParticle);
     421  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
    371422  G4double x;
    372   if(currentLoss) x = currentLoss->GetRangeForLoss(kineticEnergy, couple);
    373   else    
    374     x = G4EnergyLossTables::GetRange(currentParticle,kineticEnergy,couple,false);
     423  if(currentLoss) { x = currentLoss->GetRangeForLoss(kineticEnergy, couple); }
     424  else { x = G4EnergyLossTables::GetRange(currentParticle,kineticEnergy,
     425                                          couple,false); }
    375426  return x;
    376427}
    … …  
    383434          const G4MaterialCutsCouple *couple)
    384435{
    385   if(aParticle != currentParticle) GetEnergyLossProcess(aParticle);
     436  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
    386437  G4double x;
    387   if(currentLoss) x = currentLoss->GetRange(kineticEnergy, couple);
    388   else    
    389     x = G4EnergyLossTables::GetRange(currentParticle,kineticEnergy,couple,false);
     438  if(currentLoss) { x = currentLoss->GetRange(kineticEnergy, couple); }
     439  else { x = G4EnergyLossTables::GetRange(currentParticle,kineticEnergy,
     440                                          couple,false); }
    390441  return x;
    391442}
    … …  
    398449          const G4MaterialCutsCouple *couple)
    399450{
    400   if(aParticle != currentParticle) GetEnergyLossProcess(aParticle);
     451  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
    401452  G4double x;
    402   if(currentLoss) x = currentLoss->GetKineticEnergy(range, couple);
    403   else            x = G4EnergyLossTables::GetPreciseEnergyFromRange(
    404                       currentParticle,range,couple,false);
     453  if(currentLoss) { x = currentLoss->GetKineticEnergy(range, couple); }
     454  else { x = G4EnergyLossTables::GetPreciseEnergyFromRange(currentParticle,range,
     455                                                           couple,false); }
    405456  return x;
    406457}
    … …  
    428479//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    429480
    430 inline G4VEnergyLossProcess* G4LossTableManager::GetEnergyLossProcess(
    431           const G4ParticleDefinition *aParticle)
    432 {
    433   if(aParticle != currentParticle) {
    434     currentParticle = aParticle;
    435     std::map<PD,G4VEnergyLossProcess*,std::less<PD> >::const_iterator pos;
    436     if ((pos = loss_map.find(aParticle)) != loss_map.end()) {
    437       currentLoss = (*pos).second;
    438     } else {
    439       currentLoss = 0;
    440      // ParticleHaveNoLoss(aParticle);
    441     }
    442   }
    443   return currentLoss;
    444 }
    445 
    446 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    447 
    448481#endif
    449482

  • trunk/source/processes/electromagnetic/utils/include/G4VAtomDeexcitation.hh

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VAtomDeexcitation.hh,v 1.1 2009/07/09 11:42:52 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4VAtomDeexcitation.hh,v 1.3 2010/03/30 09:19:56 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6363  virtual ~G4VAtomDeexcitation();
    6464
    65   //initialization
     65  //========== initialization ==========
     66
    6667  virtual void PreparePhysicsTable(const G4ParticleDefinition&);
    6768  virtual void BuildPhysicsTable(const G4ParticleDefinition&);
    68 
    69   // Get atomic shell by shell index, used by discrete processes
    70   // (for example, photoelectric), when shell vacancy sampled by the model
    71   virtual const G4AtomicChell* GetAtomicShell(G4int Z, G4int ShellIndex);
    72 
    73   // selection of random shell for ionisation process
    74   virtual const G4AtomicShell* SelectRandomShell(const G4DynamicParticle*,
    75                                                  G4int Z);
    76 
    77   // generation of deexcitation for given atom and shell vacancy
    78   virtual void GenerateParticles(std::vector<G4DynamicParticle*>*,
    79                                  const G4AtomicChell*, G4int Z);
    80 
    81   // access or compute PIXE cross section
    82   virtual G4double GetPIXECrossSection (const G4ParticleDefinition*,
    83                                         G4int Z, G4double kinE);
    84 
    85   // calculate PIXE cross section from the models
    86   virtual G4double CalculatePIXECrossSection(const G4ParticleDefinition*,
    87                                              G4int Z, G4double kinE);
    88 
    89   // Sampling of PIXE for ionisation processes
    90   virtual void
    91   AlongStepDeexcitation(std::vector<G4DynamicParticle*>* secVect,
    92                         const G4DynamicParticle* icidentParticle,
    93                         const G4MaterialCutsCouple*,
    94                         G4double trueStepLenght,
    95                         G4double eLoss);
    96 
    97   // Check if deexcitation is active for a given geometry volume
    98   G4bool CheckActiveRegion(G4int coupleIndex);
    99 
    100   // Access flags defined in the CheckActiveVolume method
    101   inline G4bool IsFluorescenceActive() const;
    102   inline G4bool IsPIXECrossSectionActive() const;
    10369
    10470  // PIXE model name
    … …  
    11581  void SetPIXECrossSectionActiveRegion(const G4String& rname = "");
    11682
     83  //========== Run time methods ==========
     84
     85  // Check if deexcitation is active for a given geometry volume
     86  G4bool CheckFluorescenceActiveRegion(G4int coupleIndex);
     87
     88  // Check if deexcitation is active for a given geometry volume
     89  G4bool CheckPIXEActiveRegion(G4int coupleIndex);
     90
     91  // Get atomic shell by shell index, used by discrete processes
     92  // (for example, photoelectric), when shell vacancy sampled by the model
     93  const G4AtomicShell* GetAtomicShell(G4int Z, G4int ShellIndex);
     94
     95  // selection of random shell for ionisation process
     96  virtual const G4AtomicShell* SelectRandomShell(const G4DynamicParticle*,
     97                                                 G4int Z);
     98
     99  // generation of deexcitation for given atom and shell vacancy
     100  virtual void GenerateParticles(std::vector<G4DynamicParticle*>*,
     101                                 const G4AtomicShell*, G4int Z) = 0;
     102
     103  // access or compute PIXE cross section
     104  virtual G4double GetPIXECrossSection (const G4ParticleDefinition*,
     105                                        G4int Z, G4double kinE) = 0;
     106
     107  // calculate PIXE cross section from the models
     108  virtual G4double CalculatePIXECrossSection(const G4ParticleDefinition*,
     109                                             G4int Z, G4double kinE) = 0;
     110
     111  // Sampling of PIXE for ionisation processes
     112  virtual void
     113  AlongStepDeexcitation(std::vector<G4DynamicParticle*>* secVect,
     114                        const G4DynamicParticle* icidentParticle,
     115                        const G4MaterialCutsCouple*,
     116                        G4double trueStepLenght,
     117                        G4double eLoss) = 0;
     118
     119
     120
    117121private:
    118122
    … …  
    122126
    123127  G4String namePIXE;
    124   G4bool isFluoActive;
    125   G4bool isPIXEActive;
    126128
    127129};
    128 
    129 inline G4bool IsFluorescenceActive() const
    130 {
    131   return isFluoActive;
    132 }
    133 
    134 inline G4bool IsPIXECrossSectionActive() const
    135 {
    136   return isPIXEActive;
    137 }
    138130
    139131inline

  • trunk/source/processes/electromagnetic/utils/include/G4VEmModel.hh

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmModel.hh,v 1.72 2009/09/23 14:42:47 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4VEmModel.hh,v 1.75 2010/05/26 10:41:34 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    150150
    151151  // Compute effective ion charge square
     152  virtual G4double ChargeSquareRatio(const G4Track&);
     153
     154  // Compute effective ion charge square
    152155  virtual G4double GetChargeSquareRatio(const G4ParticleDefinition*,
    153156                                        const G4Material*,
    … …  
    228231                                             G4double maxEnergy = DBL_MAX);
    229232
     233  // select isotope in order to have precise mass of the nucleus
     234  inline G4int SelectIsotopeNumber(const G4Element*);
     235
    230236  // atom can be selected effitiantly if element selectors are initialised
    231237  inline const G4Element* SelectRandomAtom(const G4MaterialCutsCouple*,
    … …  
    236242
    237243  // to select atom cross section per volume is recomputed for each element
    238   inline const G4Element* SelectRandomAtom(const G4Material*,
    239                                            const G4ParticleDefinition*,
    240                                            G4double kineticEnergy,
    241                                            G4double cutEnergy = 0.0,
    242                                            G4double maxEnergy = DBL_MAX);
    243 
    244   // select isotope in order to have precise mass of the nucleus
    245   inline G4int SelectIsotopeNumber(const G4Element*);
     244  const G4Element* SelectRandomAtom(const G4Material*,
     245                                    const G4ParticleDefinition*,
     246                                    G4double kineticEnergy,
     247                                    G4double cutEnergy = 0.0,
     248                                    G4double maxEnergy = DBL_MAX);
    246249
    247250  //------------------------------------------------------------------------
    … …  
    291294  inline void SetCurrentCouple(const G4MaterialCutsCouple*);
    292295
     296  inline const G4Element* GetCurrentElement() const;
     297
    293298protected:
    294299
    … …  
    296301
    297302  inline void SetCurrentElement(const G4Element*);
    298 
    299   inline const G4Element* GetCurrentElement() const;
    300303
    301304private:
    … …  
    341344};
    342345
    343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     346// ======== Run time inline methods ================
     347
     348inline void G4VEmModel::SetCurrentCouple(const G4MaterialCutsCouple* p)
     349{
     350  currentCouple = p;
     351}
     352
     353//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     354
     355inline const G4MaterialCutsCouple* G4VEmModel::CurrentCouple() const
     356{
     357  return currentCouple;
     358}
     359
     360//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     361
     362inline void G4VEmModel::SetCurrentElement(const G4Element* elm)
     363{
     364  currentElement = elm;
     365}
     366
     367//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     368
     369inline const G4Element* G4VEmModel::GetCurrentElement() const
     370{
     371  return currentElement;
     372}
     373
     374//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     375
     376inline
     377G4double G4VEmModel::MaxSecondaryKinEnergy(const G4DynamicParticle* dynPart)
     378{
     379  return MaxSecondaryEnergy(dynPart->GetDefinition(),
     380                            dynPart->GetKineticEnergy());
     381}
     382
    344383//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    345384
    … …  
    375414  G4double mfp = DBL_MAX;
    376415  G4double cross = CrossSectionPerVolume(material,p,ekin,emin,emax);
    377   if (cross > DBL_MIN) mfp = 1./cross;
     416  if (cross > DBL_MIN) { mfp = 1./cross; }
    378417  return mfp;
    379418}
    … …  
    415454//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    416455
    417 inline
    418 const G4Element* G4VEmModel::SelectRandomAtom(const G4Material* material,
    419                                               const G4ParticleDefinition* pd,
    420                                               G4double kinEnergy,
    421                                               G4double tcut,
    422                                               G4double tmax)
    423 {
    424   const G4ElementVector* theElementVector = material->GetElementVector();
    425   G4int n = material->GetNumberOfElements() - 1;
    426   currentElement = (*theElementVector)[n];
    427   if (n > 0) {
    428     G4double x = G4UniformRand()*
    429                  G4VEmModel::CrossSectionPerVolume(material,pd,kinEnergy,tcut,tmax);
    430     for(G4int i=0; i<n; i++) {
    431       if (x <= xsec[i]) {
    432         currentElement = (*theElementVector)[i];
    433         break;
    434       }
    435     }
    436   }
    437   return currentElement;
    438 }
    439 
    440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    441 
    442456inline G4int G4VEmModel::SelectIsotopeNumber(const G4Element* elm)
    443457{
    … …  
    450464      G4double* ab = elm->GetRelativeAbundanceVector();
    451465      G4double x = G4UniformRand();
    452       for(; idx<ni; idx++) {
     466      for(; idx<ni; ++idx) {
    453467        x -= ab[idx];
    454         if (x <= 0.0) break;
     468        if (x <= 0.0) { break; }
    455469      }
    456       if(idx >= ni) idx = ni - 1;
     470      if(idx >= ni) { idx = ni - 1; }
    457471    }
    458472    N = elm->GetIsotope(idx)->GetN();
    … …  
    461475}
    462476
    463 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     477// ======== Get/Set inline methods used at initialisation ================
    464478
    465479inline G4VEmFluctuationModel* G4VEmModel::GetModelOfFluctuations()
    … …  
    578592{
    579593  nuclearStopping = val;
    580 }
    581 
    582 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    583 
    584 inline
    585 G4double G4VEmModel::MaxSecondaryKinEnergy(const G4DynamicParticle* dynPart)
    586 {
    587   return MaxSecondaryEnergy(dynPart->GetDefinition(),
    588                             dynPart->GetKineticEnergy());
    589594}
    590595
    … …  
    601606                                          G4VEmFluctuationModel* f = 0)
    602607{
    603   if(p && pParticleChange != p) pParticleChange = p;
     608  if(p && pParticleChange != p) { pParticleChange = p; }
    604609  fluc = f;
    605610}
    … …  
    607612//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    608613
    609 inline void G4VEmModel::SetCurrentCouple(const G4MaterialCutsCouple* p)
    610 {
    611   currentCouple = p;
    612 }
    613 
    614 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    615 
    616 inline const G4MaterialCutsCouple* G4VEmModel::CurrentCouple() const
    617 {
    618   return currentCouple;
    619 }
    620 
    621 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    622 
    623 inline void G4VEmModel::SetCurrentElement(const G4Element* elm)
    624 {
    625   currentElement = elm;
    626 }
    627 
    628 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    629 
    630 inline const G4Element* G4VEmModel::GetCurrentElement() const
    631 {
    632   return currentElement;
    633 }
    634 
    635 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    636 
    637614#endif
    638615

  • trunk/source/processes/electromagnetic/utils/include/G4VEmProcess.hh

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmProcess.hh,v 1.55 2009/09/23 14:42:47 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4VEmProcess.hh,v 1.60 2010/04/28 14:43:13 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5656// 27-10-07 Virtual functions moved to source (V.Ivanchenko)
    5757// 15-07-08 Reorder class members for further multi-thread development (VI)
     58// 17-02-10 Added pointer currentParticle (VI)
    5859//
    5960// Class Description:
    … …  
    158159
    159160  // It returns the cross section of the process per atom
    160   inline G4double ComputeCrossSectionPerAtom(G4double kineticEnergy,
    161                                              G4double Z, G4double A=0.,
    162                                              G4double cut=0.0);
    163 
    164   inline G4double MeanFreePath(const G4Track& track);
     161  G4double ComputeCrossSectionPerAtom(G4double kineticEnergy,
     162                                      G4double Z, G4double A=0.,
     163                                      G4double cut=0.0);
     164
     165  G4double MeanFreePath(const G4Track& track);
    165166
    166167  // It returns cross section per volume
    … …  
    225226  G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false);
    226227
     228  // access atom on which interaction happens
     229  const G4Element* GetCurrentElement() const;
     230
    227231  inline void SetLambdaFactor(G4double val);
    228232
    … …  
    231235
    232236  inline void SetApplyCuts(G4bool val);
     237
     238  inline void SetBuildTableFlag(G4bool val);
    233239
    234240  //------------------------------------------------------------------------
    … …  
    245251
    246252  inline G4double RecalculateLambda(G4double kinEnergy,
    247                                     const G4MaterialCutsCouple* couple);
     253                                    const G4MaterialCutsCouple* couple);
    248254
    249255  inline G4ParticleChangeForGamma* GetParticleChange();
    … …  
    258264
    259265  inline G4double GetElectronEnergyCut();
    260 
    261   inline void SetBuildTableFlag(G4bool val);
    262266
    263267  inline void SetStartFromNullFlag(G4bool val);
    … …  
    340344
    341345  const G4ParticleDefinition*  particle;
     346  const G4ParticleDefinition*  currentParticle;
    342347
    343348  // cash
    … …  
    352357};
    353358
    354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    356 
    357 inline G4double G4VEmProcess::ComputeCrossSectionPerAtom(
    358                  G4double kineticEnergy, G4double Z, G4double A, G4double cut)
    359 {
    360   SelectModel(kineticEnergy, currentCoupleIndex);
    361   G4double x = 0.0;
    362   if(currentModel) {
    363    x = currentModel->ComputeCrossSectionPerAtom(particle,kineticEnergy,
    364                                                  Z,A,cut);
    365   }
    366   return x;
    367 }
    368 
    369 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    370 
    371 inline G4double G4VEmProcess::MeanFreePath(const G4Track& track)
    372 {
    373   DefineMaterial(track.GetMaterialCutsCouple());
    374   preStepLambda = GetCurrentLambda(track.GetKineticEnergy());
    375   G4double x = DBL_MAX;
    376   if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
    377   return x;
    378 }
    379 
    380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    381 
    382 inline G4double G4VEmProcess::GetLambda(G4double& kineticEnergy,
    383                                         const G4MaterialCutsCouple* couple)
    384 {
    385   DefineMaterial(couple);
    386   return GetCurrentLambda(kineticEnergy);
    387 }
    388 
    389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    390 
    391 inline void G4VEmProcess::SetLambdaBinning(G4int nbins)
    392 {
    393   nLambdaBins = nbins;
    394 }
    395 
    396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    397 
    398 inline G4int G4VEmProcess::LambdaBinning() const
    399 {
    400   return nLambdaBins;
    401 }
    402 
    403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    404 
    405 inline void G4VEmProcess::SetMinKinEnergy(G4double e)
    406 {
    407   minKinEnergy = e;
    408 }
    409 
    410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    411 
    412 inline G4double G4VEmProcess::MinKinEnergy() const
    413 {
    414   return minKinEnergy;
    415 }
    416 
    417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    418 
    419 inline void G4VEmProcess::SetMaxKinEnergy(G4double e)
    420 {
    421   maxKinEnergy = e;
    422 }
    423 
    424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    425 
    426 inline G4double G4VEmProcess::MaxKinEnergy() const
    427 {
    428   return maxKinEnergy;
    429 }
    430 
    431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    432 
    433 inline void G4VEmProcess::SetPolarAngleLimit(G4double val)
    434 {
    435   if(val < 0.0)     polarAngleLimit = 0.0;
    436   else if(val > pi) polarAngleLimit = pi;
    437   else              polarAngleLimit = val;
    438 }
    439 
    440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    441 
    442 inline G4double G4VEmProcess::PolarAngleLimit() const
    443 {
    444   return polarAngleLimit;
    445 }
    446 
    447 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    448 
    449 inline const G4PhysicsTable* G4VEmProcess::LambdaTable() const
    450 {
    451   return theLambdaTable;
    452 }
    453 
    454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    455 
    456 inline const G4ParticleDefinition* G4VEmProcess::Particle() const
    457 {
    458   return particle;
    459 }
    460 
    461 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    462 
    463 inline const G4ParticleDefinition* G4VEmProcess::SecondaryParticle() const
    464 {
    465   return secondaryParticle;
    466 }
    467 
    468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    469 
    470 inline
    471 G4VEmModel* G4VEmProcess::SelectModel(G4double& kinEnergy, size_t index)
    472 {
    473   currentModel = modelManager->SelectModel(kinEnergy, index);
    474   currentModel->SetCurrentCouple(currentCouple);
    475   return currentModel;
    476 }
    477 
    478 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    479 
    480 inline
    481 G4VEmModel* G4VEmProcess::SelectModelForMaterial(G4double kinEnergy,
    482                                                  size_t& idxRegion) const
    483 {
    484   return modelManager->SelectModel(kinEnergy, idxRegion);
    485 }
    486 
    487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    488 
    489 inline void G4VEmProcess::SetLambdaFactor(G4double val)
    490 {
    491   if(val > 0.0 && val <= 1.0) lambdaFactor = val;
    492 }
    493 
    494 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    495 
    496 inline void G4VEmProcess::SetIntegral(G4bool val)
    497 {
    498   if(particle && particle != theGamma) integral = val;
    499   if(integral) buildLambdaTable = true;
    500 }
    501 
    502 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    503 
    504 inline G4bool G4VEmProcess::IsIntegral() const
    505 {
    506   return integral;
    507 }
    508 
    509 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    510 
    511 inline void G4VEmProcess::SetApplyCuts(G4bool val)
    512 {
    513   applyCuts = val;
    514 }
    515 
    516 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    517 
    518 inline G4double G4VEmProcess::RecalculateLambda(G4double e,
    519                                                 const G4MaterialCutsCouple* couple)
    520 {
    521   DefineMaterial(couple);
    522   return ComputeCurrentLambda(e);
    523 }
    524 
    525 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    526 
    527 inline G4ParticleChangeForGamma* G4VEmProcess::GetParticleChange()
    528 {
    529   return &fParticleChange;
    530 }
    531 
    532 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    533 
    534 inline void G4VEmProcess::SetParticle(const G4ParticleDefinition* p)
    535 {
    536   particle = p;
    537 }
    538 
    539 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    540 
    541 inline void G4VEmProcess::SetSecondaryParticle(const G4ParticleDefinition* p)
    542 {
    543   secondaryParticle = p;
    544 }
    545 
    546 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     359// ======== Run time inline methods ================
    547360
    548361inline size_t G4VEmProcess::CurrentMaterialCutsCoupleIndex() const
    … …  
    563376{
    564377  return (*theCutsElectron)[currentCoupleIndex];
    565 }
    566 
    567 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    568 
    569 inline void G4VEmProcess::SetBuildTableFlag(G4bool val)
    570 {
    571   buildLambdaTable = val;
    572   if(!val) integral = false;
    573 }
    574 
    575 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    576 
    577 inline void G4VEmProcess::SetStartFromNullFlag(G4bool val)
    578 {
    579   startFromNull = val;
    580 }
    581 
    582 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    583 
    584 inline void G4VEmProcess::InitialiseStep(const G4Track& track)
    585 {
    586   preStepKinEnergy = track.GetKineticEnergy();
    587   DefineMaterial(track.GetMaterialCutsCouple());
    588   SelectModel(preStepKinEnergy, currentCoupleIndex);
    589   if (theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX;
    590378}
    591379
    … …  
    600388    mfpKinEnergy = DBL_MAX;
    601389  }
     390}
     391
     392//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     393
     394inline
     395G4VEmModel* G4VEmProcess::SelectModel(G4double& kinEnergy, size_t index)
     396{
     397  currentModel = modelManager->SelectModel(kinEnergy, index);
     398  currentModel->SetCurrentCouple(currentCouple);
     399  return currentModel;
     400}
     401
     402//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     403
     404inline
     405G4VEmModel* G4VEmProcess::SelectModelForMaterial(G4double kinEnergy,
     406                                                 size_t& idxRegion) const
     407{
     408  return modelManager->SelectModel(kinEnergy, idxRegion);
     409}
     410
     411//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     412
     413inline void G4VEmProcess::InitialiseStep(const G4Track& track)
     414{
     415  currentParticle = track.GetDefinition();
     416  preStepKinEnergy = track.GetKineticEnergy();
     417  DefineMaterial(track.GetMaterialCutsCouple());
     418  SelectModel(preStepKinEnergy, currentCoupleIndex);
     419  if (theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX;
     420}
     421
     422//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     423
     424inline G4double G4VEmProcess::GetLambdaFromTable(G4double e)
     425{
     426  return (((*theLambdaTable)[currentCoupleIndex])->Value(e));
     427}
     428
     429//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     430
     431inline G4double G4VEmProcess::ComputeCurrentLambda(G4double e)
     432{
     433  SelectModel(e, currentCoupleIndex);
     434  return currentModel->CrossSectionPerVolume(currentMaterial,currentParticle,
     435                                             e,(*theCuts)[currentCoupleIndex]);
     436}
     437
     438//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     439
     440inline G4double G4VEmProcess::GetCurrentLambda(G4double e)
     441{
     442  G4double x = 0.0;
     443  if(theLambdaTable) { x = GetLambdaFromTable(e); }
     444  else               { x = ComputeCurrentLambda(e); }
     445  return x;
     446}
     447
     448//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     449
     450inline G4double G4VEmProcess::GetLambda(G4double& kineticEnergy,
     451                                        const G4MaterialCutsCouple* couple)
     452{
     453  DefineMaterial(couple);
     454  return GetCurrentLambda(kineticEnergy);
     455}
     456
     457//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     458
     459inline G4double G4VEmProcess::RecalculateLambda(G4double e,
     460                                                const G4MaterialCutsCouple* couple)
     461{
     462  DefineMaterial(couple);
     463  return ComputeCurrentLambda(e);
    602464}
    603465
    … …  
    625487}
    626488
    627 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    628 
    629 inline G4double G4VEmProcess::GetLambdaFromTable(G4double e)
    630 {
    631   return (((*theLambdaTable)[currentCoupleIndex])->Value(e));
    632 }
    633 
    634 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    635 
    636 inline G4double G4VEmProcess::GetCurrentLambda(G4double e)
    637 {
    638   G4double x = 0.0;
    639   if(theLambdaTable) { x = GetLambdaFromTable(e); }
    640   else               { x = ComputeCurrentLambda(e); }
    641   return x;
    642 }
    643 
    644 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    645 
    646 inline G4double G4VEmProcess::ComputeCurrentLambda(G4double e)
    647 {
    648   SelectModel(e, currentCoupleIndex);
    649   return currentModel->CrossSectionPerVolume(currentMaterial,particle,
    650                                              e,(*theCuts)[currentCoupleIndex]);
     489// ======== Get/Set inline methods used at initialisation ================
     490
     491inline void G4VEmProcess::SetLambdaBinning(G4int nbins)
     492{
     493  nLambdaBins = nbins;
     494}
     495
     496//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     497
     498inline G4int G4VEmProcess::LambdaBinning() const
     499{
     500  return nLambdaBins;
     501}
     502
     503//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     504
     505inline void G4VEmProcess::SetMinKinEnergy(G4double e)
     506{
     507  minKinEnergy = e;
     508}
     509
     510//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     511
     512inline G4double G4VEmProcess::MinKinEnergy() const
     513{
     514  return minKinEnergy;
     515}
     516
     517//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     518
     519inline void G4VEmProcess::SetMaxKinEnergy(G4double e)
     520{
     521  maxKinEnergy = e;
     522}
     523
     524//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     525
     526inline G4double G4VEmProcess::MaxKinEnergy() const
     527{
     528  return maxKinEnergy;
     529}
     530
     531//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     532
     533inline void G4VEmProcess::SetPolarAngleLimit(G4double val)
     534{
     535  if(val < 0.0)     polarAngleLimit = 0.0;
     536  else if(val > pi) polarAngleLimit = pi;
     537  else              polarAngleLimit = val;
     538}
     539
     540//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     541
     542inline G4double G4VEmProcess::PolarAngleLimit() const
     543{
     544  return polarAngleLimit;
     545}
     546
     547//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     548
     549inline const G4PhysicsTable* G4VEmProcess::LambdaTable() const
     550{
     551  return theLambdaTable;
     552}
     553
     554//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     555
     556inline const G4ParticleDefinition* G4VEmProcess::Particle() const
     557{
     558  return particle;
     559}
     560
     561//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     562
     563inline const G4ParticleDefinition* G4VEmProcess::SecondaryParticle() const
     564{
     565  return secondaryParticle;
     566}
     567
     568//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     569
     570inline void G4VEmProcess::SetLambdaFactor(G4double val)
     571{
     572  if(val > 0.0 && val <= 1.0) { lambdaFactor = val; }
     573}
     574
     575//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     576
     577inline void G4VEmProcess::SetIntegral(G4bool val)
     578{
     579  if(particle && particle != theGamma) { integral = val; }
     580  if(integral) { buildLambdaTable = true; }
     581}
     582
     583//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     584
     585inline G4bool G4VEmProcess::IsIntegral() const
     586{
     587  return integral;
     588}
     589
     590//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     591
     592inline void G4VEmProcess::SetApplyCuts(G4bool val)
     593{
     594  applyCuts = val;
     595}
     596
     597//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     598
     599inline void G4VEmProcess::SetBuildTableFlag(G4bool val)
     600{
     601  buildLambdaTable = val;
     602  if(!val) { integral = false; }
     603}
     604
     605//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     606
     607inline G4ParticleChangeForGamma* G4VEmProcess::GetParticleChange()
     608{
     609  return &fParticleChange;
     610}
     611
     612//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     613
     614inline void G4VEmProcess::SetParticle(const G4ParticleDefinition* p)
     615{
     616  particle = p;
     617  currentParticle = p;
     618}
     619
     620//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     621
     622inline void G4VEmProcess::SetSecondaryParticle(const G4ParticleDefinition* p)
     623{
     624  secondaryParticle = p;
     625}
     626
     627//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     628
     629inline void G4VEmProcess::SetStartFromNullFlag(G4bool val)
     630{
     631  startFromNull = val;
    651632}
    652633

  • trunk/source/processes/electromagnetic/utils/include/G4VEnergyLossProcess.hh

    r1196 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEnergyLossProcess.hh,v 1.89 2009/07/03 14:39:17 vnivanch Exp $
     26// $Id: G4VEnergyLossProcess.hh,v 1.92 2010/04/28 14:43:13 vnivanch Exp $
    2727// GEANT4 tag $Name:
    2828//
    … …  
    196196  // Sampling of secondaries in vicinity of geometrical boundary
    197197  void SampleSubCutSecondaries(std::vector<G4Track*>&, const G4Step&,
    198                                G4VEmModel* model, G4int matIdx,
    199                                G4double& extraEdep);
     198                               G4VEmModel* model, G4int matIdx);
    200199
    201200  // PostStep sampling of secondaries
    … …  
    418417  // Run time method for simulation of ionisation
    419418  //------------------------------------------------------------------------
     419
     420  // access atom on which interaction happens
     421  const G4Element* GetCurrentElement() const;
    420422
    421423  // sample range at the end of a step
    … …  
    553555};
    554556
    555 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    556 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     557// ======== Run time inline methods ================
    557558
    558559inline size_t G4VEnergyLossProcess::CurrentMaterialCutsCoupleIndex() const
    … …  
    586587//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    587588
    588 inline void G4VEnergyLossProcess::SetFluctModel(G4VEmFluctuationModel* p)
    589 {
    590   fluctModel = p;
    591 }
    592 
    593 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    594 
    595 inline G4VEmFluctuationModel* G4VEnergyLossProcess::FluctModel()
    596 {
    597   return fluctModel;
    598 }
    599 
    600 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    601 
    602 inline void G4VEnergyLossProcess::SetParticle(const G4ParticleDefinition* p)
    603 {
    604   particle = p;
    605 }
    606 
    607 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    608 
    609 inline void G4VEnergyLossProcess::SetSecondaryParticle(const G4ParticleDefinition* p)
    610 {
    611   secondaryParticle = p;
    612 }
    613 
    614 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    615 
    616 inline void G4VEnergyLossProcess::SetBaseParticle(const G4ParticleDefinition* p)
    617 {
    618   baseParticle = p;
    619 }
    620 
    621 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    622 
    623 inline const G4ParticleDefinition* G4VEnergyLossProcess::Particle() const
    624 {
    625   return particle;
    626 }
    627 
    628 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    629 
    630 inline const G4ParticleDefinition* G4VEnergyLossProcess::BaseParticle() const
    631 {
    632   return baseParticle;
    633 }
    634 
    635 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    636 
    637 inline const G4ParticleDefinition* G4VEnergyLossProcess::SecondaryParticle() const
    638 {
    639   return secondaryParticle;
    640 }
    641 
    642 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    643 
    644 inline void G4VEnergyLossProcess::SetLossFluctuations(G4bool val)
    645 {
    646   lossFluctuationFlag = val;
    647 }
    648 
    649 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    650 
    651 inline void G4VEnergyLossProcess::SetRandomStep(G4bool val)
    652 {
    653   rndmStepFlag = val;
    654 }
    655 
    656 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    657 
    658 inline void G4VEnergyLossProcess::SetIntegral(G4bool val)
    659 {
    660   integral = val;
    661 }
    662 
    663 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    664  
    665 inline G4bool G4VEnergyLossProcess::IsIntegral() const
    666 {
    667   return integral;
    668 }
    669 
    670 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    671 
    672 inline void G4VEnergyLossProcess::SetIonisation(G4bool val)
    673 {
    674   isIonisation = val;
    675   if(val) aGPILSelection = CandidateForSelection;
    676   else    aGPILSelection = NotCandidateForSelection;
    677 }
    678 
    679 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    680 
    681 inline G4bool G4VEnergyLossProcess::IsIonisationProcess() const
    682 {
    683   return isIonisation;
    684 }
    685 
    686 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    687 
    688 inline void G4VEnergyLossProcess::SetLinearLossLimit(G4double val)
    689 {
    690   linLossLimit = val;
    691 }
    692 
    693 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    694 
    695 inline void G4VEnergyLossProcess::SetMinSubRange(G4double val)
    696 {
    697   minSubRange = val;
    698 }
    699 
    700 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    701 
    702 inline void G4VEnergyLossProcess::SetLambdaFactor(G4double val)
    703 {
    704   if(val > 0.0 && val <= 1.0) lambdaFactor = val;
    705 }
    706 
    707 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    708 
    709 void G4VEnergyLossProcess::SetStepFunction(G4double v1, G4double v2)
    710 {
    711   dRoverRange = v1;
    712   finalRange = v2;
    713   if (dRoverRange > 0.999) dRoverRange = 1.0;
    714   currentCouple = 0;
    715   mfpKinEnergy  = DBL_MAX;
    716 }
    717 
    718 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    719 
    720 inline void G4VEnergyLossProcess::SetLowestEnergyLimit(G4double val)
    721 {
    722   lowestKinEnergy = val;
    723 }
    724 
    725 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    726 
    727 inline G4int G4VEnergyLossProcess::NumberOfSubCutoffRegions() const
    728 {
    729   return nSCoffRegions;
    730 }
    731 
    732 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    733 
    734 inline G4int G4VEnergyLossProcess::NumberOfDERegions() const
    735 {
    736   return nDERegions;
    737 }
    738 
    739 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    740 
    741 inline void G4VEnergyLossProcess::SetDEDXBinning(G4int nbins)
    742 {
    743   nBins = nbins;
    744 }
    745 
    746 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    747 
    748 inline void G4VEnergyLossProcess::SetLambdaBinning(G4int nbins)
    749 {
    750   nBins = nbins;
    751 }
    752 
    753 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    754 
    755 inline void G4VEnergyLossProcess::SetDEDXBinningForCSDARange(G4int nbins)
    756 {
    757   nBinsCSDA = nbins;
    758 }
    759 
    760 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    761 
    762 inline void G4VEnergyLossProcess::SetMinKinEnergy(G4double e)
    763 {
    764   minKinEnergy = e;
    765 }
    766 
    767 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    768 
    769 inline G4double G4VEnergyLossProcess::MinKinEnergy() const
    770 {
    771   return minKinEnergy;
    772 }
    773 
    774 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    775 
    776 inline void G4VEnergyLossProcess::SetMaxKinEnergy(G4double e)
    777 {
    778   maxKinEnergy = e;
    779   if(e < maxKinEnergyCSDA) maxKinEnergyCSDA = e;
    780 }
    781 
    782 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    783 
    784 inline G4double G4VEnergyLossProcess::MaxKinEnergy() const
    785 {
    786   return maxKinEnergy;
    787 }
    788 
    789 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    790 
    791 inline void G4VEnergyLossProcess::SetMaxKinEnergyForCSDARange(G4double e)
    792 {
    793   maxKinEnergyCSDA = e;
    794 }
    795 
    796 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    797 
    798 inline G4double G4VEnergyLossProcess::GetDEDX(G4double& kineticEnergy,
    799                                         const G4MaterialCutsCouple* couple)
     589inline void
     590G4VEnergyLossProcess::DefineMaterial(const G4MaterialCutsCouple* couple)
     591{
     592  if(couple != currentCouple) {
     593    currentCouple   = couple;
     594    currentMaterial = couple->GetMaterial();
     595    currentMaterialIndex = couple->GetIndex();
     596    mfpKinEnergy = DBL_MAX;
     597  }
     598}
     599
     600//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     601
     602inline void G4VEnergyLossProcess::SetDynamicMassCharge(G4double massratio,
     603                                                       G4double charge2ratio)
     604{
     605  massRatio     = massratio;
     606  chargeSqRatio = charge2ratio;
     607  reduceFactor  = 1.0/(chargeSqRatio*massRatio);
     608}
     609
     610//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     611
     612inline G4double G4VEnergyLossProcess::GetDEDXForScaledEnergy(G4double e)
     613{
     614  G4double x = ((*theDEDXTable)[currentMaterialIndex]->Value(e))*chargeSqRatio;
     615  if(e < minKinEnergy) { x *= std::sqrt(e/minKinEnergy); }
     616  return x;
     617}
     618
     619//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     620
     621inline G4double G4VEnergyLossProcess::GetSubDEDXForScaledEnergy(G4double e)
     622{
     623  G4double x = ((*theDEDXSubTable)[currentMaterialIndex]->Value(e))*chargeSqRatio;
     624  if(e < minKinEnergy) { x *= std::sqrt(e/minKinEnergy); }
     625  return x;
     626}
     627
     628//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     629
     630inline G4double G4VEnergyLossProcess::GetIonisationForScaledEnergy(G4double e)
     631{
     632  //G4double x = 0.0;
     633  //  if(theIonisationTable) {
     634  G4double x = ((*theIonisationTable)[currentMaterialIndex]->Value(e))*chargeSqRatio;
     635  if(e < minKinEnergy) { x *= std::sqrt(e/minKinEnergy); }
     636  //}
     637  return x;
     638}
     639
     640//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     641
     642inline
     643G4double G4VEnergyLossProcess::GetSubIonisationForScaledEnergy(G4double e)
     644{
     645  //  G4double x = 0.0;
     646  //if(theIonisationSubTable) {
     647  G4double x = ((*theIonisationSubTable)[currentMaterialIndex]->Value(e))*chargeSqRatio;
     648  if(e < minKinEnergy) { x *= std::sqrt(e/minKinEnergy); }
     649  //}
     650  return x;
     651}
     652
     653//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     654
     655inline G4double G4VEnergyLossProcess::GetScaledRangeForScaledEnergy(G4double e)
     656{
     657  G4double x = ((*theRangeTableForLoss)[currentMaterialIndex])->Value(e);
     658  if(e < minKinEnergy) { x *= std::sqrt(e/minKinEnergy); }
     659  return x;
     660}
     661
     662//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     663
     664inline G4double
     665G4VEnergyLossProcess::GetLimitScaledRangeForScaledEnergy(G4double e)
     666{
     667  G4double x;
     668
     669  if (e < maxKinEnergyCSDA) {
     670    x = ((*theCSDARangeTable)[currentMaterialIndex])->Value(e);
     671    if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
     672  } else {
     673    x = theRangeAtMaxEnergy[currentMaterialIndex] +
     674         (e - maxKinEnergyCSDA)/theDEDXAtMaxEnergy[currentMaterialIndex];
     675  }
     676  return x;
     677}
     678
     679//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     680
     681inline G4double G4VEnergyLossProcess::ScaledKinEnergyForLoss(G4double r)
     682{
     683  G4PhysicsVector* v = (*theInverseRangeTable)[currentMaterialIndex];
     684  G4double rmin = v->Energy(0);
     685  G4double e = 0.0;
     686  if(r >= rmin) { e = v->Value(r); }
     687  else if(r > 0.0) {
     688    G4double x = r/rmin;
     689    e = minKinEnergy*x*x;
     690  }
     691  return e;
     692}
     693
     694//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     695
     696inline G4double G4VEnergyLossProcess::GetLambdaForScaledEnergy(G4double e)
     697{
     698  return chargeSqRatio*(((*theLambdaTable)[currentMaterialIndex])->Value(e));
     699}
     700
     701//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     702
     703inline G4double
     704G4VEnergyLossProcess::GetDEDX(G4double& kineticEnergy,
     705                              const G4MaterialCutsCouple* couple)
    800706{
    801707  DefineMaterial(couple);
    … …  
    805711//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    806712
    807 inline G4double G4VEnergyLossProcess::GetDEDXForSubsec(G4double& kineticEnergy,
    808                                         const G4MaterialCutsCouple* couple)
     713inline G4double
     714G4VEnergyLossProcess::GetDEDXForSubsec(G4double& kineticEnergy,
     715                                                const G4MaterialCutsCouple* couple)
    809716{
    810717  DefineMaterial(couple);
    … …  
    814721//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    815722
    816 inline G4double G4VEnergyLossProcess::GetRange(G4double& kineticEnergy,
    817                                          const G4MaterialCutsCouple* couple)
     723inline G4double
     724G4VEnergyLossProcess::GetRange(G4double& kineticEnergy,
     725                               const G4MaterialCutsCouple* couple)
    818726{
    819727  G4double x = fRange;
    … …  
    831739//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    832740
    833 inline G4double G4VEnergyLossProcess::GetCSDARange(
    834        G4double& kineticEnergy, const G4MaterialCutsCouple* couple)
     741inline G4double
     742G4VEnergyLossProcess::GetCSDARange(G4double& kineticEnergy,
     743                                   const G4MaterialCutsCouple* couple)
    835744{
    836745  DefineMaterial(couple);
    … …  
    844753//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    845754
    846 inline G4double G4VEnergyLossProcess::GetRangeForLoss(
    847                 G4double& kineticEnergy,
    848                 const G4MaterialCutsCouple* couple)
     755inline G4double
     756G4VEnergyLossProcess::GetRangeForLoss(G4double& kineticEnergy,
     757                                      const G4MaterialCutsCouple* couple)
    849758{
    850759  DefineMaterial(couple);
    … …  
    859768//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    860769
    861 inline G4double G4VEnergyLossProcess::GetKineticEnergy(
    862                 G4double& range,
    863                 const G4MaterialCutsCouple* couple)
     770inline G4double
     771G4VEnergyLossProcess::GetKineticEnergy(G4double& range,
     772                                       const G4MaterialCutsCouple* couple)
    864773{
    865774  DefineMaterial(couple);
    … …  
    871780//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    872781
    873 inline G4double G4VEnergyLossProcess::GetLambda(G4double& kineticEnergy,
    874                                           const G4MaterialCutsCouple* couple)
     782inline G4double
     783G4VEnergyLossProcess::GetLambda(G4double& kineticEnergy,
     784                                const G4MaterialCutsCouple* couple)
    875785{
    876786  DefineMaterial(couple);
    877787  G4double x = 0.0;
    878   if(theLambdaTable) x = GetLambdaForScaledEnergy(kineticEnergy*massRatio);
     788  if(theLambdaTable) { x = GetLambdaForScaledEnergy(kineticEnergy*massRatio); }
    879789  return x;
    880 }
    881 
    882 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    883 
    884 inline G4bool G4VEnergyLossProcess::TablesAreBuilt() const
    885 {
    886   return  tablesAreBuilt;
    887 }
    888 
    889 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    890 
    891 inline G4PhysicsTable* G4VEnergyLossProcess::DEDXTable() const
    892 {
    893   return theDEDXTable;
    894 }
    895 
    896 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    897 
    898 inline G4PhysicsTable* G4VEnergyLossProcess::DEDXTableForSubsec() const
    899 {
    900   return theDEDXSubTable;
    901 }
    902 
    903 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    904 
    905 inline G4PhysicsTable* G4VEnergyLossProcess::DEDXunRestrictedTable() const
    906 {
    907   return theDEDXunRestrictedTable;
    908 }
    909 
    910 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    911 
    912 inline G4PhysicsTable* G4VEnergyLossProcess::IonisationTable() const
    913 {
    914   G4PhysicsTable* t = theDEDXTable;
    915   if(theIonisationTable) t = theIonisationTable;
    916   return t;
    917 }
    918 
    919 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    920 
    921 inline G4PhysicsTable* G4VEnergyLossProcess::IonisationTableForSubsec() const
    922 {
    923   G4PhysicsTable* t = theDEDXSubTable;
    924   if(theIonisationSubTable) t = theIonisationSubTable;
    925   return t;
    926 }
    927 
    928 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    929 
    930 inline G4PhysicsTable* G4VEnergyLossProcess::CSDARangeTable() const
    931 {
    932   return theCSDARangeTable;
    933 }
    934 
    935 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    936 
    937 inline G4PhysicsTable* G4VEnergyLossProcess::RangeTableForLoss() const
    938 {
    939   return theRangeTableForLoss;
    940 }
    941 
    942 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    943 
    944 inline G4PhysicsTable* G4VEnergyLossProcess::InverseRangeTable() const
    945 {
    946   return theInverseRangeTable;
    947 }
    948 
    949 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    950 
    951 inline G4PhysicsTable* G4VEnergyLossProcess::LambdaTable()
    952 {
    953   return theLambdaTable;
    954 }
    955 
    956 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    957 
    958 inline G4PhysicsTable* G4VEnergyLossProcess::SubLambdaTable()
    959 {
    960   return theSubLambdaTable;
    961 }
    962 
    963 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    964 
    965 inline G4double G4VEnergyLossProcess::SampleRange()
    966 {
    967   G4double e = amu_c2*preStepKinEnergy/particle->GetPDGMass();
    968   G4double s = fRange*std::pow(10.,vstrag->Value(e));
    969   G4double x = fRange + G4RandGauss::shoot(0.0,s);
    970   if(x > 0.0) fRange = x;
    971   return fRange;
    972 }
    973 
    974 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    975 
    976 inline void G4VEnergyLossProcess::SetDynamicMassCharge(G4double massratio,
    977                                                        G4double charge2ratio)
    978 {
    979   massRatio     = massratio;
    980   chargeSqRatio = charge2ratio;
    981   reduceFactor  = 1.0/(chargeSqRatio*massRatio);
    982 }
    983 
    984 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    985 
    986 inline void G4VEnergyLossProcess::DefineMaterial(
    987             const G4MaterialCutsCouple* couple)
    988 {
    989   if(couple != currentCouple) {
    990     currentCouple   = couple;
    991     currentMaterial = couple->GetMaterial();
    992     currentMaterialIndex = couple->GetIndex();
    993     mfpKinEnergy = DBL_MAX;
    994   }
    995 }
    996 
    997 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    998 
    999 inline G4double G4VEnergyLossProcess::GetDEDXForScaledEnergy(G4double e)
    1000 {
    1001   G4double x = ((*theDEDXTable)[currentMaterialIndex]->Value(e))*chargeSqRatio;
    1002   if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
    1003   return x;
    1004 }
    1005 
    1006 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1007 
    1008 inline G4double G4VEnergyLossProcess::GetSubDEDXForScaledEnergy(G4double e)
    1009 {
    1010   G4double x = ((*theDEDXSubTable)[currentMaterialIndex]->Value(e))*chargeSqRatio;
    1011   if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
    1012   return x;
    1013 }
    1014 
    1015 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1016 
    1017 inline G4double G4VEnergyLossProcess::GetIonisationForScaledEnergy(G4double e)
    1018 {
    1019   //G4double x = 0.0;
    1020   //  if(theIonisationTable) {
    1021   G4double x = ((*theIonisationTable)[currentMaterialIndex]->Value(e))*chargeSqRatio;
    1022   if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
    1023   //}
    1024   return x;
    1025 }
    1026 
    1027 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1028 
    1029 inline
    1030 G4double G4VEnergyLossProcess::GetSubIonisationForScaledEnergy(G4double e)
    1031 {
    1032   //  G4double x = 0.0;
    1033   //if(theIonisationSubTable) {
    1034   G4double x = ((*theIonisationSubTable)[currentMaterialIndex]->Value(e))*chargeSqRatio;
    1035   if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
    1036   //}
    1037   return x;
    1038 }
    1039 
    1040 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1041 
    1042 inline G4double G4VEnergyLossProcess::GetScaledRangeForScaledEnergy(G4double e)
    1043 {
    1044   G4double x = ((*theRangeTableForLoss)[currentMaterialIndex])->Value(e);
    1045   if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
    1046   return x;
    1047 }
    1048 
    1049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1050 
    1051 inline G4double G4VEnergyLossProcess::GetLimitScaledRangeForScaledEnergy(
    1052                 G4double e)
    1053 {
    1054   G4double x;
    1055 
    1056   if (e < maxKinEnergyCSDA) {
    1057     x = ((*theCSDARangeTable)[currentMaterialIndex])->Value(e);
    1058     if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
    1059   } else {
    1060     x = theRangeAtMaxEnergy[currentMaterialIndex] +
    1061          (e - maxKinEnergyCSDA)/theDEDXAtMaxEnergy[currentMaterialIndex];
    1062   }
    1063   return x;
    1064 }
    1065 
    1066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1067 
    1068 inline G4double G4VEnergyLossProcess::ScaledKinEnergyForLoss(G4double r)
    1069 {
    1070   G4PhysicsVector* v = (*theInverseRangeTable)[currentMaterialIndex];
    1071   G4double rmin = v->Energy(0);
    1072   G4double e = 0.0;
    1073   if(r >= rmin) { e = v->Value(r); }
    1074   else if(r > 0.0) {
    1075     G4double x = r/rmin;
    1076     e = minKinEnergy*x*x;
    1077   }
    1078   return e;
    1079 }
    1080 
    1081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1082 
    1083 inline G4double G4VEnergyLossProcess::GetLambdaForScaledEnergy(G4double e)
    1084 {
    1085   return chargeSqRatio*(((*theLambdaTable)[currentMaterialIndex])->Value(e));
    1086790}
    1087791
    … …  
    1111815//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1112816
     817inline G4double G4VEnergyLossProcess::SampleRange()
     818{
     819  G4double e = amu_c2*preStepKinEnergy/particle->GetPDGMass();
     820  G4double s = fRange*std::pow(10.,vstrag->Value(e));
     821  G4double x = fRange + G4RandGauss::shoot(0.0,s);
     822  if(x > 0.0) { fRange = x; }
     823  return fRange;
     824}
     825
     826// ======== Get/Set inline methods used at initialisation ================
     827
     828inline void G4VEnergyLossProcess::SetFluctModel(G4VEmFluctuationModel* p)
     829{
     830  fluctModel = p;
     831}
     832
     833//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     834
     835inline G4VEmFluctuationModel* G4VEnergyLossProcess::FluctModel()
     836{
     837  return fluctModel;
     838}
     839
     840//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     841
     842inline void G4VEnergyLossProcess::SetParticle(const G4ParticleDefinition* p)
     843{
     844  particle = p;
     845}
     846
     847//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     848
     849inline void G4VEnergyLossProcess::SetSecondaryParticle(const G4ParticleDefinition* p)
     850{
     851  secondaryParticle = p;
     852}
     853
     854//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     855
     856inline void G4VEnergyLossProcess::SetBaseParticle(const G4ParticleDefinition* p)
     857{
     858  baseParticle = p;
     859}
     860
     861//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     862
     863inline const G4ParticleDefinition* G4VEnergyLossProcess::Particle() const
     864{
     865  return particle;
     866}
     867
     868//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     869
     870inline const G4ParticleDefinition* G4VEnergyLossProcess::BaseParticle() const
     871{
     872  return baseParticle;
     873}
     874
     875//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     876
     877inline const G4ParticleDefinition* G4VEnergyLossProcess::SecondaryParticle() const
     878{
     879  return secondaryParticle;
     880}
     881
     882//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     883
     884inline void G4VEnergyLossProcess::SetLossFluctuations(G4bool val)
     885{
     886  lossFluctuationFlag = val;
     887}
     888
     889//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     890
     891inline void G4VEnergyLossProcess::SetRandomStep(G4bool val)
     892{
     893  rndmStepFlag = val;
     894}
     895
     896//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     897
     898inline void G4VEnergyLossProcess::SetIntegral(G4bool val)
     899{
     900  integral = val;
     901}
     902
     903//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     904 
     905inline G4bool G4VEnergyLossProcess::IsIntegral() const
     906{
     907  return integral;
     908}
     909
     910//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     911
     912inline void G4VEnergyLossProcess::SetIonisation(G4bool val)
     913{
     914  isIonisation = val;
     915  if(val) { aGPILSelection = CandidateForSelection; }
     916  else    { aGPILSelection = NotCandidateForSelection; }
     917}
     918
     919//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     920
     921inline G4bool G4VEnergyLossProcess::IsIonisationProcess() const
     922{
     923  return isIonisation;
     924}
     925
     926//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     927
     928inline void G4VEnergyLossProcess::SetLinearLossLimit(G4double val)
     929{
     930  linLossLimit = val;
     931}
     932
     933//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     934
     935inline void G4VEnergyLossProcess::SetMinSubRange(G4double val)
     936{
     937  minSubRange = val;
     938}
     939
     940//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     941
     942inline void G4VEnergyLossProcess::SetLambdaFactor(G4double val)
     943{
     944  if(val > 0.0 && val <= 1.0) { lambdaFactor = val; }
     945}
     946
     947//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     948
     949void G4VEnergyLossProcess::SetStepFunction(G4double v1, G4double v2)
     950{
     951  dRoverRange = v1;
     952  finalRange = v2;
     953  if (dRoverRange > 0.999) { dRoverRange = 1.0; }
     954  currentCouple = 0;
     955  mfpKinEnergy  = DBL_MAX;
     956}
     957
     958//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     959
     960inline void G4VEnergyLossProcess::SetLowestEnergyLimit(G4double val)
     961{
     962  lowestKinEnergy = val;
     963}
     964
     965//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     966
     967inline G4int G4VEnergyLossProcess::NumberOfSubCutoffRegions() const
     968{
     969  return nSCoffRegions;
     970}
     971
     972//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     973
     974inline G4int G4VEnergyLossProcess::NumberOfDERegions() const
     975{
     976  return nDERegions;
     977}
     978
     979//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     980
     981inline void G4VEnergyLossProcess::SetDEDXBinning(G4int nbins)
     982{
     983  nBins = nbins;
     984}
     985
     986//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     987
     988inline void G4VEnergyLossProcess::SetLambdaBinning(G4int nbins)
     989{
     990  nBins = nbins;
     991}
     992
     993//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     994
     995inline void G4VEnergyLossProcess::SetDEDXBinningForCSDARange(G4int nbins)
     996{
     997  nBinsCSDA = nbins;
     998}
     999
     1000//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1001
     1002inline void G4VEnergyLossProcess::SetMinKinEnergy(G4double e)
     1003{
     1004  minKinEnergy = e;
     1005}
     1006
     1007//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1008
     1009inline G4double G4VEnergyLossProcess::MinKinEnergy() const
     1010{
     1011  return minKinEnergy;
     1012}
     1013
     1014//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1015
     1016inline void G4VEnergyLossProcess::SetMaxKinEnergy(G4double e)
     1017{
     1018  maxKinEnergy = e;
     1019  if(e < maxKinEnergyCSDA) { maxKinEnergyCSDA = e; }
     1020}
     1021
     1022//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1023
     1024inline G4double G4VEnergyLossProcess::MaxKinEnergy() const
     1025{
     1026  return maxKinEnergy;
     1027}
     1028
     1029//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1030
     1031inline void G4VEnergyLossProcess::SetMaxKinEnergyForCSDARange(G4double e)
     1032{
     1033  maxKinEnergyCSDA = e;
     1034}
     1035
     1036
     1037//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1038
     1039inline G4bool G4VEnergyLossProcess::TablesAreBuilt() const
     1040{
     1041  return  tablesAreBuilt;
     1042}
     1043
     1044//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1045
     1046inline G4PhysicsTable* G4VEnergyLossProcess::DEDXTable() const
     1047{
     1048  return theDEDXTable;
     1049}
     1050
     1051//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1052
     1053inline G4PhysicsTable* G4VEnergyLossProcess::DEDXTableForSubsec() const
     1054{
     1055  return theDEDXSubTable;
     1056}
     1057
     1058//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1059
     1060inline G4PhysicsTable* G4VEnergyLossProcess::DEDXunRestrictedTable() const
     1061{
     1062  return theDEDXunRestrictedTable;
     1063}
     1064
     1065//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1066
     1067inline G4PhysicsTable* G4VEnergyLossProcess::IonisationTable() const
     1068{
     1069  G4PhysicsTable* t = theDEDXTable;
     1070  if(theIonisationTable) { t = theIonisationTable; }
     1071  return t;
     1072}
     1073
     1074//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1075
     1076inline G4PhysicsTable* G4VEnergyLossProcess::IonisationTableForSubsec() const
     1077{
     1078  G4PhysicsTable* t = theDEDXSubTable;
     1079  if(theIonisationSubTable) { t = theIonisationSubTable; }
     1080  return t;
     1081}
     1082
     1083//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1084
     1085inline G4PhysicsTable* G4VEnergyLossProcess::CSDARangeTable() const
     1086{
     1087  return theCSDARangeTable;
     1088}
     1089
     1090//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1091
     1092inline G4PhysicsTable* G4VEnergyLossProcess::RangeTableForLoss() const
     1093{
     1094  return theRangeTableForLoss;
     1095}
     1096
     1097//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1098
     1099inline G4PhysicsTable* G4VEnergyLossProcess::InverseRangeTable() const
     1100{
     1101  return theInverseRangeTable;
     1102}
     1103
     1104//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1105
     1106inline G4PhysicsTable* G4VEnergyLossProcess::LambdaTable()
     1107{
     1108  return theLambdaTable;
     1109}
     1110
     1111//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1112
     1113inline G4PhysicsTable* G4VEnergyLossProcess::SubLambdaTable()
     1114{
     1115  return theSubLambdaTable;
     1116}
     1117
     1118//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1119
    11131120#endif

  • trunk/source/processes/electromagnetic/utils/include/G4VMultipleScattering.hh

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMultipleScattering.hh,v 1.62 2009/10/29 17:56:04 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4VMultipleScattering.hh,v 1.63 2010/03/10 18:29:51 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    154154  // The function overloads the corresponding function of the base
    155155  // class.
    156   inline G4double PostStepGetPhysicalInteractionLength(
     156  G4double PostStepGetPhysicalInteractionLength(
    157157                                            const G4Track&,
    158158                                            G4double  previousStepSize,
    … …  
    160160
    161161  // Along step actions
    162   inline G4VParticleChange* AlongStepDoIt(const G4Track&, const G4Step&);
     162  G4VParticleChange* AlongStepDoIt(const G4Track&, const G4Step&);
    163163
    164164  // Post step actions
    165   inline G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&);
     165  G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&);
    166166
    167167  // This method does not used for tracking, it is intended only for tests
    168   inline G4double ContinuousStepLimit(const G4Track& track,
    169                                       G4double previousStepSize,
    170                                       G4double currentMinimalStep,
    171                                       G4double& currentSafety);
     168  G4double ContinuousStepLimit(const G4Track& track,
     169                               G4double previousStepSize,
     170                               G4double currentMinimalStep,
     171                               G4double& currentSafety);
    172172
    173173  //------------------------------------------------------------------------
    … …  
    318318};
    319319
    320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    322 
    323 inline G4double G4VMultipleScattering::ContinuousStepLimit(
    324                                        const G4Track& track,
    325                                        G4double previousStepSize,
    326                                        G4double currentMinimalStep,
    327                                        G4double& currentSafety)
    328 {
    329   return GetContinuousStepLimit(track,previousStepSize,currentMinimalStep,
    330                                 currentSafety);
    331 }
    332 
    333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    334 
    335 inline void G4VMultipleScattering::SetBinning(G4int nbins)
    336 {
    337   nBins = nbins;
    338 }
    339 
    340 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    341 
    342 inline G4int G4VMultipleScattering::Binning() const
    343 {
    344   return nBins;
    345 }
    346 
    347 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    348 
    349 inline void G4VMultipleScattering::SetMinKinEnergy(G4double e)
    350 {
    351   minKinEnergy = e;
    352 }
    353 
    354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    355 
    356 inline G4double G4VMultipleScattering::MinKinEnergy() const
    357 {
    358   return minKinEnergy;
    359 }
    360 
    361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    362 
    363 inline void G4VMultipleScattering::SetMaxKinEnergy(G4double e)
    364 {
    365   maxKinEnergy = e;
    366 }
    367 
    368 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    369 
    370 inline G4double G4VMultipleScattering::MaxKinEnergy() const
    371 {
    372   return maxKinEnergy;
    373 }
    374 
    375 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    376 
    377 inline  void G4VMultipleScattering::SetBuildLambdaTable(G4bool val)
    378 {
    379   buildLambdaTable = val;
    380 }
    381 
    382 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    383 
    384 inline G4PhysicsTable* G4VMultipleScattering::LambdaTable() const
    385 {
    386   return theLambdaTable;
    387 }
    388 
    389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    390 
    391 inline  const G4ParticleDefinition* G4VMultipleScattering::Particle() const
    392 {
    393   return currentParticle;
    394 }
    395 
    396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    397 
    398 inline G4VEmModel* G4VMultipleScattering::SelectModel(G4double kinEnergy)
    399 {
    400   return modelManager->SelectModel(kinEnergy, currentMaterialIndex);
    401 }
    402 
    403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    404 
    405 inline G4VEmModel* G4VMultipleScattering::SelectModelForMaterial(
    406                    G4double kinEnergy, size_t& idxRegion) const
    407 {
    408   return modelManager->SelectModel(kinEnergy, idxRegion);
    409 }
    410 
    411 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    412 
    413 inline  G4bool G4VMultipleScattering::LateralDisplasmentFlag() const
    414 {
    415   return latDisplasment;
    416 }
    417 
    418 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    419 
    420 inline  void G4VMultipleScattering::SetLateralDisplasmentFlag(G4bool val)
    421 {
    422   latDisplasment = val;
    423 }
    424 
    425 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    426 
    427 inline  G4double G4VMultipleScattering::Skin() const
    428 {
    429   return skin;
    430 }
    431 
    432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    433 
    434 inline  void G4VMultipleScattering::SetSkin(G4double val)
    435 {
    436   if(val < 1.0) skin = 0.0;
    437   else          skin = val;
    438 }
    439 
    440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    441 
    442 inline  G4double G4VMultipleScattering::RangeFactor() const
    443 {
    444   return facrange;
    445 }
    446 
    447 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    448 
    449 inline  void G4VMultipleScattering::SetRangeFactor(G4double val)
    450 {
    451   if(val > 0.0) facrange = val;
    452 }
    453 
    454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    455 
    456 inline  G4double G4VMultipleScattering::GeomFactor() const
    457 {
    458   return facgeom;
    459 }
    460 
    461 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    462 
    463 inline  void G4VMultipleScattering::SetGeomFactor(G4double val)
    464 {
    465   if(val > 0.0) facgeom = val;
    466 }
    467 
    468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    469 
    470 inline  G4double G4VMultipleScattering::PolarAngleLimit() const
    471 {
    472   return polarAngleLimit;
    473 }
    474 
    475 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    476 
    477 inline  void G4VMultipleScattering::SetPolarAngleLimit(G4double val)
    478 {
    479   if(val < 0.0)     polarAngleLimit = 0.0;
    480   else if(val > pi) polarAngleLimit = pi;
    481   else              polarAngleLimit = val;
    482 }
    483 
    484 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    485 
    486 inline G4MscStepLimitType G4VMultipleScattering::StepLimitType() const
    487 {
    488   return stepLimit;
    489 }
    490 
    491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    492 
    493 inline void G4VMultipleScattering::SetStepLimitType(G4MscStepLimitType val)
    494 {
    495   stepLimit = val;
    496   if(val == fMinimal) facrange = 0.2;
     320// ======== Run time inline methods ================
     321
     322inline const G4MaterialCutsCouple*
     323G4VMultipleScattering::CurrentMaterialCutsCouple() const
     324{
     325  return currentCouple;
     326}
     327
     328//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     329
     330inline
     331void G4VMultipleScattering::DefineMaterial(const G4MaterialCutsCouple* couple)
     332{
     333  if(couple != currentCouple) {
     334    currentCouple   = couple;
     335    currentMaterialIndex = couple->GetIndex();
     336  }
    497337}
    498338
    … …  
    509349    x = currentModel->CrossSection(currentCouple,p,e);
    510350  }
    511   if(x > DBL_MIN) x = 1./x;
    512   else            x = DBL_MAX;
     351  if(x > DBL_MIN) { x = 1./x; }
     352  else            { x = DBL_MAX; }
    513353  return x;
    514354}
    … …  
    516356//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    517357
    518 inline
    519 void G4VMultipleScattering::DefineMaterial(const G4MaterialCutsCouple* couple)
    520 {
    521   if(couple != currentCouple) {
    522     currentCouple   = couple;
    523     currentMaterialIndex = couple->GetIndex();
    524   }
    525 }
    526 
    527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    528 
    529 inline const G4MaterialCutsCouple*
    530 G4VMultipleScattering::CurrentMaterialCutsCouple() const
    531 {
    532   return currentCouple;
    533 }
    534 
    535 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    536 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    537 
    538 // Follwoing methods are virtual, they are inlined because they applied at
    539 // each simulation step and some compilers may inline these methods
    540 
    541 inline G4double
    542 G4VMultipleScattering::PostStepGetPhysicalInteractionLength(
    543               const G4Track&, G4double, G4ForceCondition* condition)
    544 {
    545   *condition = Forced;
    546   return DBL_MAX;
    547 }
    548 
    549 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    550 
    551 inline G4VParticleChange*
    552 G4VMultipleScattering::AlongStepDoIt(const G4Track& track, const G4Step& step)
    553 {
    554   if(currentModel->IsActive(track.GetKineticEnergy())) {
    555     fParticleChange.ProposeTrueStepLength(currentModel->ComputeTrueStepLength(step.GetStepLength()));
    556   }
    557   return &fParticleChange;
    558 }
    559 
    560 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    561 
    562 inline G4VParticleChange*
    563 G4VMultipleScattering::PostStepDoIt(const G4Track& track, const G4Step& step)
    564 {
    565   fParticleChange.Initialize(track);
    566   if(currentModel->IsActive(track.GetKineticEnergy())) {
    567     currentModel->SampleScattering(track.GetDynamicParticle(),
    568                                    step.GetPostStepPoint()->GetSafety());
    569   }
    570   return &fParticleChange;
     358inline G4VEmModel* G4VMultipleScattering::SelectModel(G4double kinEnergy)
     359{
     360  return modelManager->SelectModel(kinEnergy, currentMaterialIndex);
     361}
     362
     363//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     364
     365inline G4VEmModel* G4VMultipleScattering::SelectModelForMaterial(
     366                   G4double kinEnergy, size_t& idxRegion) const
     367{
     368  return modelManager->SelectModel(kinEnergy, idxRegion);
     369}
     370
     371// ======== Get/Set inline methods used at initialisation ================
     372
     373inline void G4VMultipleScattering::SetBinning(G4int nbins)
     374{
     375  nBins = nbins;
     376}
     377
     378//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     379
     380inline G4int G4VMultipleScattering::Binning() const
     381{
     382  return nBins;
     383}
     384
     385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     386
     387inline void G4VMultipleScattering::SetMinKinEnergy(G4double e)
     388{
     389  minKinEnergy = e;
     390}
     391
     392//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     393
     394inline G4double G4VMultipleScattering::MinKinEnergy() const
     395{
     396  return minKinEnergy;
     397}
     398
     399//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     400
     401inline void G4VMultipleScattering::SetMaxKinEnergy(G4double e)
     402{
     403  maxKinEnergy = e;
     404}
     405
     406//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     407
     408inline G4double G4VMultipleScattering::MaxKinEnergy() const
     409{
     410  return maxKinEnergy;
     411}
     412
     413//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     414
     415inline  void G4VMultipleScattering::SetBuildLambdaTable(G4bool val)
     416{
     417  buildLambdaTable = val;
     418}
     419
     420//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     421
     422inline G4PhysicsTable* G4VMultipleScattering::LambdaTable() const
     423{
     424  return theLambdaTable;
     425}
     426
     427//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     428
     429inline  const G4ParticleDefinition* G4VMultipleScattering::Particle() const
     430{
     431  return currentParticle;
     432}
     433
     434//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     435
     436inline  G4bool G4VMultipleScattering::LateralDisplasmentFlag() const
     437{
     438  return latDisplasment;
     439}
     440
     441//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     442
     443inline  void G4VMultipleScattering::SetLateralDisplasmentFlag(G4bool val)
     444{
     445  latDisplasment = val;
     446}
     447
     448//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     449
     450inline  G4double G4VMultipleScattering::Skin() const
     451{
     452  return skin;
     453}
     454
     455//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     456
     457inline  void G4VMultipleScattering::SetSkin(G4double val)
     458{
     459  if(val < 1.0) { skin = 0.0; }
     460  else          { skin = val; }
     461}
     462
     463//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     464
     465inline  G4double G4VMultipleScattering::RangeFactor() const
     466{
     467  return facrange;
     468}
     469
     470//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     471
     472inline  void G4VMultipleScattering::SetRangeFactor(G4double val)
     473{
     474  if(val > 0.0) facrange = val;
     475}
     476
     477//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     478
     479inline  G4double G4VMultipleScattering::GeomFactor() const
     480{
     481  return facgeom;
     482}
     483
     484//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     485
     486inline  void G4VMultipleScattering::SetGeomFactor(G4double val)
     487{
     488  if(val > 0.0) facgeom = val;
     489}
     490
     491//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     492
     493inline  G4double G4VMultipleScattering::PolarAngleLimit() const
     494{
     495  return polarAngleLimit;
     496}
     497
     498//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     499
     500inline  void G4VMultipleScattering::SetPolarAngleLimit(G4double val)
     501{
     502  if(val < 0.0)            { polarAngleLimit = 0.0; }
     503  else if(val > CLHEP::pi) { polarAngleLimit = CLHEP::pi; }
     504  else                     { polarAngleLimit = val; }
     505}
     506
     507//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     508
     509inline G4MscStepLimitType G4VMultipleScattering::StepLimitType() const
     510{
     511  return stepLimit;
     512}
     513
     514//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     515
     516inline void G4VMultipleScattering::SetStepLimitType(G4MscStepLimitType val)
     517{
     518  stepLimit = val;
     519  if(val == fMinimal) { facrange = 0.2; }
    571520}
    572521

  • trunk/source/processes/electromagnetic/utils/src/G4EmCalculator.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmCalculator.cc,v 1.49 2009/11/22 17:58:39 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4EmCalculator.cc,v 1.53 2010/04/13 10:58:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9494  currentMaterial    = 0;
    9595  currentParticle    = 0;
     96  lambdaParticle     = 0;
    9697  baseParticle       = 0;
    9798  currentLambda      = 0;
    … …  
    102103  currentParticleName= "";
    103104  currentMaterialName= "";
     105  currentName        = "";
     106  lambdaName         = "";
    104107  theGenericIon      = G4GenericIon::GenericIon();
    105108  ionEffCharge       = new G4ionEffectiveCharge();
    … …  
    113116{
    114117  delete ionEffCharge;
    115   for (G4int i=0; i<nLocalMaterials; i++) {
     118  for (G4int i=0; i<nLocalMaterials; ++i) {
    116119    delete localCouples[i];
    117120  }
    … …  
    310313    G4int idx = couple->GetIndex();
    311314    FindLambdaTable(p, processName);
     315
    312316    if(currentLambda) {
    313317      G4double e = kinEnergy*massRatio;
    314318      res = (((*currentLambda)[idx])->Value(e))*chargeSquare;
    315319      if(verbose>0) {
    316         G4cout << "E(MeV)= " << kinEnergy/MeV
     320        G4cout << "G4EmCalculator::GetXSPerVolume: E(MeV)= " << kinEnergy/MeV
    317321               << " cross(cm-1)= " << res*cm
    318322               << "  " <<  p->GetParticleName()
    319323               << " in " <<  mat->GetName();
    320324        if(verbose>1)
    321           G4cout << "  idx= " << idx << "  e(MeV)= " << e
     325          G4cout << "  idx= " << idx << "  Escaled((MeV)= " << e
    322326                 << "  q2= " << chargeSquare;
    323327        G4cout << G4endl;
    … …  
    350354  G4double res = DBL_MAX;
    351355  G4double x = GetCrossSectionPerVolume(kinEnergy,p, processName, mat,region);
    352   if(x > 0.0) res = 1.0/x;
     356  if(x > 0.0) { res = 1.0/x; }
    353357  if(verbose>1) {
    354358    G4cout << "G4EmCalculator::GetMeanFreePath: E(MeV)= " << kinEnergy/MeV
    … …  
    495499             << " in " <<  currentMaterialName
    496500             << " Zi^2= " << chargeSquare
     501             << " isIon=" << isIon
    497502             << G4endl;
    498503    }
    … …  
    516521      lManager->GetEnergyLossProcessVector();
    517522    G4int n = vel.size();
    518     for(G4int i=0; i<n; i++) {
     523    for(G4int i=0; i<n; ++i) {
    519524      const G4ParticleDefinition* p = (vel[i])->Particle();
    520       if((!isIon && p == part) || (isIon && p == theGenericIon))
     525      if((!isIon && p == part) || (isIon && p == theGenericIon)) {
    521526        dedx += ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),mat,cut);
     527      }
    522528    }
    523529  }
    … …  
    619625      }
    620626      if(verbose>0) {
    621         G4cout << "E(MeV)= " << kinEnergy/MeV
     627        G4cout << "G4EmCalculator::ComputeXSPerVolume: E(MeV)= " << kinEnergy/MeV
    622628               << " cross(cm-1)= " << res*cm
     629               << " cut(keV)= " << cut/keV
    623630               << "  " <<  p->GetParticleName()
    624631               << " in " <<  mat->GetName()
    … …  
    686693  return ComputeCrossSectionPerAtom(kinEnergy,FindParticle(particle),
    687694                                    processName,
    688                                     elm->GetZ(),elm->GetA(),cut);
     695                                    elm->GetZ(),elm->GetN(),cut);
    689696}
    690697
    … …  
    886893  currentMaterial = material;
    887894  currentMaterialName = material->GetName();
    888   for (G4int i=0; i<nLocalMaterials; i++) {
     895  for (G4int i=0; i<nLocalMaterials; ++i) {
    889896    if(material == localMaterials[i] && cut == localCuts[i]) {
    890897      currentCouple = localCouples[i];
    … …  
    911918{
    912919  // Search for the process
    913   if (p != currentParticle || processName != currentName) {
    914     currentName     = processName;
    915     currentLambda   = 0;
     920  if (!currentLambda || p != lambdaParticle || processName != lambdaName) {
     921    lambdaName     = processName;
     922    currentLambda  = 0;
     923    lambdaParticle = p;
    916924
    917925    G4String partname =  p->GetParticleName();
    … …  
    924932    lManager->GetEnergyLossProcessVector();
    925933    G4int n = vel.size();
    926     for(G4int i=0; i<n; i++) {
    927       if((vel[i])->GetProcessName() == currentName &&
     934    for(G4int i=0; i<n; ++i) {
     935      if((vel[i])->GetProcessName() == lambdaName &&
    928936         (vel[i])->Particle() == part)
    929937        {
    930938          currentLambda = (vel[i])->LambdaTable();
    931939          isApplicable    = true;
    932           break;
     940          if(verbose>1) {
     941            G4cout << "G4VEnergyLossProcess is found out: "
     942                   << currentName << G4endl;
     943          }
     944          return;
    933945        }
    934946    }
    … …  
    938950      const std::vector<G4VEmProcess*> vem = lManager->GetEmProcessVector();
    939951      G4int n = vem.size();
    940       for(G4int i=0; i<n; i++) {
    941         if((vem[i])->GetProcessName() == currentName &&
     952      for(G4int i=0; i<n; ++i) {
     953        if((vem[i])->GetProcessName() == lambdaName &&
    942954           (vem[i])->Particle() == part)
    943955        {
    944956          currentLambda = (vem[i])->LambdaTable();
    945957          isApplicable    = true;
    946           break;
     958          if(verbose>1) {
     959            G4cout << "G4VEmProcess is found out: "
     960                   << currentName << G4endl;
     961          }
     962          return;
    947963        }
    948964      }
    … …  
    954970        lManager->GetMultipleScatteringVector();
    955971      G4int n = vmsc.size();
    956       for(G4int i=0; i<n; i++) {
    957         if((vmsc[i])->GetProcessName() == currentName &&
     972      for(G4int i=0; i<n; ++i) {
     973        if((vmsc[i])->GetProcessName() == lambdaName &&
    958974           (vmsc[i])->Particle() == part)
    959975        {
    960976          currentLambda = (vmsc[i])->LambdaTable();
    961977          isApplicable    = true;
    962           break;
     978          if(verbose>1) {
     979            G4cout << "G4VMultipleScattering is found out: "
     980                   << currentName << G4endl;
     981          }
     982          return;
    963983        }
    964984      }
    … …  
    10021022  G4int n = vel.size();
    10031023  G4VEnergyLossProcess* elproc = 0;
    1004   for(G4int i=0; i<n; i++) {
     1024  for(G4int i=0; i<n; ++i) {
    10051025    //    G4cout << "i= " << i << " part= "
    10061026    //  << (vel[i])->Particle()->GetParticleName()
    … …  
    10231043    currentModel = elproc->SelectModelForMaterial(scaledEnergy, idx);
    10241044    G4double eth = currentModel->LowEnergyLimit();
    1025     loweModel = elproc->SelectModelForMaterial(eth - CLHEP::eV, idx);
     1045    if(eth > 0.0) {
     1046      loweModel = elproc->SelectModelForMaterial(eth - CLHEP::eV, idx);
     1047    }
    10261048  }
    10271049
    … …  
    10301052    const std::vector<G4VEmProcess*> vem = lManager->GetEmProcessVector();
    10311053    G4int n = vem.size();
    1032     for(G4int i=0; i<n; i++) {
     1054    for(G4int i=0; i<n; ++i) {
    10331055      if((vem[i])->GetProcessName() == currentName &&
    10341056         (vem[i])->Particle() == part)
    … …  
    10361058        currentModel = (vem[i])->SelectModelForMaterial(kinEnergy, idx);
    10371059        G4double eth = currentModel->LowEnergyLimit();
    1038         loweModel = (vem[i])->SelectModelForMaterial(eth - CLHEP::eV, idx);
     1060        if(eth > 0.0) {
     1061          loweModel = (vem[i])->SelectModelForMaterial(eth - CLHEP::eV, idx);
     1062        }
    10391063        break;
    10401064      }
    … …  
    10471071      lManager->GetMultipleScatteringVector();
    10481072    G4int n = vmsc.size();
    1049     for(G4int i=0; i<n; i++) {
     1073    for(G4int i=0; i<n; ++i) {
    10501074      if((vmsc[i])->GetProcessName() == currentName &&
    10511075         (vmsc[i])->Particle() == part)
    … …  
    10531077        currentModel = (vmsc[i])->SelectModelForMaterial(kinEnergy, idx);
    10541078        G4double eth = currentModel->LowEnergyLimit();
    1055         loweModel = (vmsc[i])->SelectModelForMaterial(eth - CLHEP::eV, idx); 
     1079        if(eth > 0.0) {
     1080          loweModel = (vmsc[i])->SelectModelForMaterial(eth - CLHEP::eV, idx);
     1081        }
    10561082        break;
    10571083      }
    … …  
    10821108  const G4ParticleDefinition* part = p;
    10831109 
    1084   if(p->GetParticleType() == "nucleus" &&
    1085      partname != "deuteron" &&
    1086      partname != "triton") { part = theGenericIon; }
     1110  if(p->GetParticleType() == "nucleus"
     1111     && currentParticleName != "deuteron" 
     1112     && currentParticleName != "triton"
     1113     && currentParticleName != "alpha+"
     1114     && currentParticleName != "helium"
     1115     && currentParticleName != "hydrogen"
     1116     ) { part = theGenericIon; }
    10871117 
    10881118  G4LossTableManager* lManager = G4LossTableManager::Instance();
    … …  
    10901120    lManager->GetEnergyLossProcessVector();
    10911121  G4int n = vel.size();
    1092   for(G4int i=0; i<n; i++) {
     1122  for(G4int i=0; i<n; ++i) {
    10931123    if( (vel[i])->Particle() == part ) {
    10941124      elp = vel[i];

  • trunk/source/processes/electromagnetic/utils/src/G4EmConfigurator.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmConfigurator.cc,v 1.6 2009/11/22 19:48:30 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4EmConfigurator.cc,v 1.8 2010/06/04 15:33:56 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6161#include "G4VMultipleScattering.hh"
    6262
    63 enum PType {unknown=0, eloss, discrete, msc};
    64 
    65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    66 
    67 G4EmConfigurator::G4EmConfigurator()
     63//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     64
     65G4EmConfigurator::G4EmConfigurator(G4int val):verbose(val)
    6866{
    6967  index = -10;
    … …  
    7472G4EmConfigurator::~G4EmConfigurator()
    7573{}
    76 
    77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    78 
    79 void G4EmConfigurator::AddExtraEmModel(const G4String& particleName,
    80                                        G4VEmModel* em,
    81                                        G4VEmFluctuationModel* fm)
    82 {
    83   particleList.push_back(particleName);
    84   modelList.push_back(em);
    85   flucModelList.push_back(fm);
    86 }
    87 
    88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    89 
    90 void G4EmConfigurator::AddModelForRegion(const G4String& particleName,
    91                                          const G4String& processName,
    92                                          const G4String& modelName,
    93                                          const G4String& regionName,
    94                                          G4double emin, G4double emax,
    95                                          const G4String& flucModelName)
    96 {
    97   particles.push_back(particleName);
    98   processes.push_back(processName);
    99   models.push_back(modelName);
    100   regions.push_back(regionName);
    101   flucModels.push_back(flucModelName);
    102   lowEnergy.push_back(emin);
    103   highEnergy.push_back(emax);
    104 }
    10574
    10675//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    11483                                       G4VEmFluctuationModel* fm)
    11584{
    116   AddExtraEmModel(particleName, mod, fm);
    117   G4String fname = "";
    118   if(fm) fname = fm->GetName();
    119   G4String mname = "";
    120   if(mod) mname = mod->GetName();
    121   AddModelForRegion(particleName, processName, mname, regionName,
    122                     emin, emax, fname);
     85  if(1 < verbose) {
     86    G4cout << " G4EmConfigurator::SetExtraEmModel " << mod->GetName()
     87           << " for " << particleName
     88           << " and " << processName
     89           << " in the region <" << regionName
     90           << "> Emin(MeV)= " << emin/MeV
     91           << " Emax(MeV)= " << emax/MeV
     92           << G4endl;
     93  }
     94  if(mod || fm) {
     95    models.push_back(mod);
     96    flucModels.push_back(fm);
     97  } else {
     98    models.push_back(new G4DummyModel());
     99    flucModels.push_back(0);
     100  }
     101
     102  particles.push_back(particleName);
     103  processes.push_back(processName);
     104  regions.push_back(regionName);
     105  lowEnergy.push_back(emin);
     106  highEnergy.push_back(emax);
    123107}
    124108
    … …  
    127111void G4EmConfigurator::AddModels()
    128112{
     113  size_t n = models.size();
     114  if(0 < verbose) {
     115    G4cout << "### G4EmConfigurator::AddModels n= " << n << G4endl;
     116  }
     117  if(n > 0) {
     118    for(size_t i=0; i<n; ++i) {
     119      if(models[i]) {
     120        G4Region* reg = FindRegion(regions[i]);
     121        if(reg) {
     122          --index;
     123          SetModelForRegion(models[i],flucModels[i],reg,
     124                            particles[i],processes[i],
     125                            lowEnergy[i],highEnergy[i]);
     126        }
     127      }
     128    }
     129  }
     130  Clear();
     131}
     132
     133//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     134
     135void G4EmConfigurator::SetModelForRegion(G4VEmModel* mod,
     136                                         G4VEmFluctuationModel* fm,
     137                                         G4Region* reg,
     138                                         const G4String& particleName,
     139                                         const G4String& processName,
     140                                         G4double emin, G4double emax)
     141{
     142  if(1 < verbose) {
     143    G4cout << " G4EmConfigurator::SetModelForRegion: " << mod->GetName()
     144           << G4endl;
     145    G4cout << " For " << particleName
     146           << " and " << processName
     147           << " in the region <" << reg->GetName()
     148           << " Emin(MeV)= " << emin/MeV
     149           << " Emax(MeV)= " << emax/MeV;
     150    if(fm) { G4cout << " FLmodel " << fm->GetName(); }
     151    G4cout << G4endl;
     152  }
     153  G4ParticleTable::G4PTblDicIterator* theParticleIterator =
     154    G4ParticleTable::GetParticleTable()->GetIterator();
     155
     156  theParticleIterator->reset();
     157  while( (*theParticleIterator)() ) {
     158    const G4ParticleDefinition* part = theParticleIterator->value();
     159
     160    //G4cout << particleName << " " << part->GetParticleName() << G4endl;
     161
     162    if((part->GetParticleName() == particleName) ||
     163       (particleName == "all") ||
     164       (particleName == "charged" && part->GetPDGCharge() != 0.0)) {
     165
     166      // search for process
     167      G4ProcessManager* pmanager = part->GetProcessManager();
     168      G4ProcessVector* plist = pmanager->GetProcessList();
     169      G4int np = pmanager->GetProcessListLength();
     170 
     171      //G4cout << processName << " in list of " << np << G4endl;
     172
     173      G4VProcess* proc = 0;
     174      for(G4int i=0; i<np; i++) {
     175        if(processName == (*plist)[i]->GetProcessName()) {
     176          proc = (*plist)[i];
     177          break;
     178        }
     179      }
     180      if(!proc) {
     181        G4cout << "### G4EmConfigurator WARNING: fails to find a process <"
     182               << processName << "> for " << particleName << G4endl;
     183        return;
     184
     185      }
     186
     187      if(mod) {
     188        if(!UpdateModelEnergyRange(mod, emin,emax)) { return; }
     189      }
     190      // classify process
     191      G4int ii = proc->GetProcessSubType();
     192      if(10 == ii && mod) {
     193        G4VMultipleScattering* p = static_cast<G4VMultipleScattering*>(proc);
     194        p->AddEmModel(index,mod,reg);
     195        if(1 < verbose) {
     196          G4cout << "### Added msc model order= " << index << " for "
     197                 << particleName << " and " << processName << G4endl;
     198        }
     199        return;
     200      } else if(2 <= ii && 4 >= ii) {
     201        G4VEnergyLossProcess* p = static_cast<G4VEnergyLossProcess*>(proc);
     202        if(!mod && fm) {
     203          p->SetFluctModel(fm);
     204        } else {
     205          p->AddEmModel(index,mod,fm,reg);
     206          if(1 < verbose) {
     207            G4cout << "### Added eloss model order= " << index << " for "
     208                   << particleName << " and " << processName << G4endl;
     209          }
     210        }
     211        return;
     212      } else if(mod) {
     213        G4VEmProcess* p = static_cast<G4VEmProcess*>(proc);
     214        p->AddEmModel(index,mod,reg);
     215        if(1 < verbose) {
     216          G4cout << "### Added em model order= " << index << " for "
     217                 << particleName << " and " << processName << G4endl;
     218        }
     219        return;
     220      } else {
     221        return;
     222      }
     223    }
     224  }
     225}
     226
     227//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     228
     229void
     230G4EmConfigurator::PrepareModels(const G4ParticleDefinition* aParticle,
     231                                G4VEnergyLossProcess* p)
     232{
    129233  size_t n = particles.size();
    130   //G4cout << " G4EmConfigurator::AddModels n= " << n << G4endl;
     234  if(1 < verbose) {
     235    G4cout << " G4EmConfigurator::PrepareModels for EnergyLoss n= "
     236           << n << G4endl;
     237  }
    131238  if(n > 0) {
    132     for(size_t i=0; i<n; i++) {
    133       SetModelForRegion(particles[i],processes[i],models[i],regions[i],
    134                         flucModels[i],lowEnergy[i],highEnergy[i]);
     239    G4String particleName = aParticle->GetParticleName();
     240    G4String processName  = p->GetProcessName();
     241    //G4cout <<  particleName << "  " <<  processName << G4endl;
     242    for(size_t i=0; i<n; ++i) {
     243      //G4cout <<  particles[i] << "  " <<  processes[i] << G4endl;
     244      if(processName == processes[i]) {
     245        if((particleName == particles[i]) ||
     246           (particles[i] == "all") ||
     247           (particles[i] == "charged" && aParticle->GetPDGCharge() != 0.0)) {
     248          G4Region* reg = FindRegion(regions[i]);
     249          //G4cout << "Region " << reg << G4endl;
     250          if(reg) {
     251            --index;
     252            G4VEmModel* mod = models[i];
     253            G4VEmFluctuationModel* fm = flucModels[i];
     254            if(mod) {
     255              if(UpdateModelEnergyRange(mod, lowEnergy[i], highEnergy[i])) {
     256                p->AddEmModel(index,mod,fm,reg);
     257                if(1 < verbose) {
     258                  G4cout << "### Added eloss model order= " << index << " for "
     259                         << particleName << " and " << processName << G4endl;
     260                }
     261              }
     262            } else if(fm) {
     263              p->SetFluctModel(fm);
     264            }
     265          }
     266        }
     267      }
    135268    }
    136269  }
     270}
     271
     272//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     273
     274void
     275G4EmConfigurator::PrepareModels(const G4ParticleDefinition* aParticle,
     276                                G4VEmProcess* p)
     277{
     278  size_t n = particles.size();
     279  if(1 < verbose) {
     280    G4cout << " G4EmConfigurator::PrepareModels for EM process n= "
     281           << n << G4endl;
     282  }
     283  if(n > 0) {
     284    G4String particleName = aParticle->GetParticleName();
     285    G4String processName  = p->GetProcessName();
     286    //G4cout <<  particleName << "  " <<  particleName << G4endl;
     287    for(size_t i=0; i<n; ++i) {
     288      if(processName == processes[i]) {
     289        if((particleName == particles[i]) ||
     290           (particles[i] == "all") ||
     291           (particles[i] == "charged" && aParticle->GetPDGCharge() != 0.0)) {
     292          G4Region* reg = FindRegion(regions[i]);
     293          //G4cout << "Region " << reg << G4endl;
     294          if(reg) {
     295            --index;
     296            G4VEmModel* mod = models[i];
     297            if(mod) {
     298              if(UpdateModelEnergyRange(mod, lowEnergy[i], highEnergy[i])) {
     299                p->AddEmModel(index,mod,reg);
     300                if(1 < verbose) {
     301                  G4cout << "### Added em model order= " << index << " for "
     302                         << particleName << " and " << processName << G4endl;
     303                }
     304              }
     305            }
     306          }
     307        }
     308      }
     309    }
     310  }
     311}
     312
     313//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     314
     315void
     316G4EmConfigurator::PrepareModels(const G4ParticleDefinition* aParticle,
     317                                G4VMultipleScattering* p)
     318{
     319  size_t n = particles.size();
     320  if(1 < verbose) {
     321    G4cout << " G4EmConfigurator::PrepareModels for MSC process n= "
     322           << n << G4endl;
     323  }
     324
     325  if(n > 0) {
     326    G4String particleName = aParticle->GetParticleName();
     327    G4String processName  = p->GetProcessName();
     328    for(size_t i=0; i<n; ++i) {
     329      if(processName == processes[i]) {
     330        if((particleName == particles[i]) ||
     331           (particles[i] == "all") ||
     332           (particles[i] == "charged" && aParticle->GetPDGCharge() != 0.0)) {
     333          G4Region* reg = FindRegion(regions[i]);
     334          if(reg) {
     335            --index;
     336            G4VEmModel* mod = models[i];
     337            if(mod) {
     338              if(UpdateModelEnergyRange(mod, lowEnergy[i], highEnergy[i])) {
     339                p->AddEmModel(index,mod,reg);
     340                G4cout << "### Added msc model order= " << index << " for "
     341                       << particleName << " and " << processName << G4endl;
     342              }
     343            }
     344          }
     345        }
     346      }
     347    }
     348  }
     349}
     350
     351//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     352
     353void G4EmConfigurator::Clear()
     354{
    137355  particles.clear();
    138356  processes.clear();
    … …  
    146364//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    147365
    148 void G4EmConfigurator::SetModelForRegion(const G4String& particleName,
    149                                          const G4String& processName,
    150                                          const G4String& modelName,
    151                                          const G4String& regionName,
    152                                          const G4String& flucModelName,
    153                                          G4double emin, G4double emax)
    154 {
    155   //G4cout << " G4EmConfigurator::SetModelForRegion" << G4endl;
    156 
    157   // new set
    158   --index;
    159 
    160   G4ParticleTable::G4PTblDicIterator* theParticleIterator =
    161     G4ParticleTable::GetParticleTable()->GetIterator();
    162 
    163   theParticleIterator->reset();
    164   while( (*theParticleIterator)() ) {
    165     const G4ParticleDefinition* part = theParticleIterator->value();
    166 
    167     //G4cout << particleName << " " << part->GetParticleName() << G4endl;
    168 
    169     if(particleName == part->GetParticleName() ||
    170        (particleName == "charged" && part->GetPDGCharge() != 0.0) ) {
    171 
    172 
    173       // search for process
    174       G4ProcessManager* pmanager = part->GetProcessManager();
    175       G4ProcessVector* plist = pmanager->GetProcessList();
    176       G4int np = pmanager->GetProcessListLength();
    177  
    178       //G4cout << processName << " in list of " << np << G4endl;
    179 
    180       G4VProcess* proc = 0;
    181       for(G4int i=0; i<np; i++) {
    182         if(processName == (*plist)[i]->GetProcessName()) {
    183           proc = (*plist)[i];
    184           break;
    185         }
    186       }
    187       if(!proc) {
    188         G4cout << "### G4EmConfigurator WARNING: fails to find a process <"
    189                << processName << "> for " << particleName << G4endl;
    190        
    191       } else {
    192 
    193         // classify process
    194         PType ptype = discrete;
    195         G4int ii = proc->GetProcessSubType();
    196         if(10 == ii) ptype = msc;
    197         else if(2 <= ii && 4 >= ii) ptype = eloss;
    198 
    199         // find out model     
    200         G4VEmModel* mod = 0;
    201         G4VEmFluctuationModel* fluc = 0;
    202 
    203         G4int nm = modelList.size();
    204         //G4cout << "Search model " << modelName << " in " << nm << G4endl;
    205 
    206         for(G4int i=0; i<nm; i++) {
    207           G4String mname = "";
    208           if(modelList[i]) mname = modelList[i]->GetName();
    209           G4String fname = "";
    210           if(flucModelList[i]) fname = flucModelList[i]->GetName();
    211           if(modelName == mname && flucModelName == fname &&
    212              (particleList[i] == "" || particleList[i] == particleName) ) {
    213             mod  = modelList[i];
    214             fluc = flucModelList[i];
    215             break;
    216           }
    217         }
    218 
    219         if("dummy" == modelName) mod = new G4DummyModel();
    220 
    221         if(!mod) {
    222 
    223           // set fluctuation model for ionisation processes
    224           if(fluc && ptype == eloss) {
    225             G4VEnergyLossProcess* p = static_cast<G4VEnergyLossProcess*>(proc);
    226             p->SetFluctModel(fluc);
    227          
    228           } else {
    229             G4cout << "### G4EmConfigurator WARNING: fails to find a model <"
    230                    << modelName << "> for process <"
    231                    << processName << "> and " << particleName
    232                    << G4endl;
    233             if(flucModelName != "") {
    234               G4cout << "                            fluctuation model <"
    235                      << flucModelName << G4endl;
    236             }
    237           }
    238         } else {
    239 
    240           // search for region
    241           G4Region* reg = 0;
    242           G4RegionStore* regStore = G4RegionStore::GetInstance();
    243           G4String r = regionName;
    244           if(r == "" || r == "world" || r == "World") r = "DefaultRegionForTheWorld";
    245           reg = regStore->GetRegion(r, true);
    246           if(!reg) {
    247             G4cout << "### G4EmConfigurator WARNING: fails to find a region <"
    248                    << r << "> for model <" << modelName << "> of the process "
    249                    << processName << " and " << particleName << G4endl;
    250             return;
    251           }
    252 
    253           // energy limits
    254           G4double e1 = std::max(emin,mod->LowEnergyLimit());
    255           G4double e2 = std::min(emax,mod->HighEnergyLimit());
    256           if(e2 < e1) e2 = e1;
    257           mod->SetLowEnergyLimit(e1);
    258           mod->SetHighEnergyLimit(e2);
    259 
    260           //G4cout << "index= " << index << " e1= " << e1 << " e2= " << e2 << G4endl;
    261 
    262           // added model
    263           if(ptype == eloss) {
    264             G4VEnergyLossProcess* p = static_cast<G4VEnergyLossProcess*>(proc);
    265             p->AddEmModel(index,mod,fluc,reg);
    266             //G4cout << "### Added eloss model order= " << index << " for "
    267             //     << particleName << " and " << processName << "  " << mod << G4endl;
    268           } else if(ptype == discrete) {
    269             G4VEmProcess* p = static_cast<G4VEmProcess*>(proc);
    270             p->AddEmModel(index,mod,reg);
    271           } else if(ptype == msc) {
    272             //G4cout << "### Added msc model order= " << index << " for "
    273             //     << particleName << " and " << processName << "  " << mod << G4endl;
    274             G4VMultipleScattering* p = static_cast<G4VMultipleScattering*>(proc);
    275             p->AddEmModel(index,mod,reg);
    276           }
    277         }
    278       }
    279     }
    280   }
    281 }
    282 
    283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    284 
    285 
    286 
    287 
    288 
    289 
    290 
    291 
     366G4Region* G4EmConfigurator::FindRegion(const G4String& regionName)
     367{
     368  // search for region
     369  G4Region* reg = 0;
     370  G4RegionStore* regStore = G4RegionStore::GetInstance();
     371  G4String r = regionName;
     372  if(r == "" || r == "world" || r == "World") {
     373    r = "DefaultRegionForTheWorld";
     374  }
     375  reg = regStore->GetRegion(r, true);
     376  if(!reg) {
     377    G4cout << "### G4EmConfigurator WARNING: fails to find a region <"
     378           << r << G4endl;
     379  } else if(verbose > 1) {
     380    G4cout << "### G4EmConfigurator finds out G4Region <" << r << ">"
     381           << G4endl;
     382  }
     383  return reg; 
     384}
     385
     386//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     387
     388G4bool G4EmConfigurator::UpdateModelEnergyRange(G4VEmModel* mod,
     389                                                G4double emin, G4double emax)
     390{
     391  // energy limits
     392  G4double e1 = std::max(emin,mod->LowEnergyLimit());
     393  G4double e2 = std::min(emax,mod->HighEnergyLimit());
     394  if(e2 <= e1) {
     395    G4cout << "### G4EmConfigurator WARNING: empty energy interval"
     396           << " for <" << mod->GetName()
     397           << ">  Emin(MeV)= " << e1/CLHEP::MeV
     398           << ">  Emax(MeV)= " << e2/CLHEP::MeV
     399           << G4endl;
     400    return false;       
     401  }
     402  mod->SetLowEnergyLimit(e1);
     403  mod->SetHighEnergyLimit(e2);
     404  if(verbose > 1) {
     405    G4cout << "### G4EmConfigurator for " << mod->GetName()
     406           << " Emin(MeV)= " << e1/MeV << " Emax(MeV)= " << e2/MeV
     407           << G4endl;
     408  }
     409  return true;
     410}
     411
     412//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/source/processes/electromagnetic/utils/src/G4EmCorrections.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmCorrections.cc,v 1.54 2009/10/29 17:56:36 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4EmCorrections.cc,v 1.58 2010/06/04 09:28:46 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6868#include "G4ProductionCutsTable.hh"
    6969#include "G4MaterialCutsCouple.hh"
     70#include "G4AtomicShells.hh"
    7071
    7172//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    123124  G4double sum = (2.0*(Barkas + Bloch) + Mott);
    124125
    125   if(verbose > 1)
     126  if(verbose > 1) {
    126127    G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas
    127128           << " Bloch= " << Bloch << " Mott= " << Mott
    128            << " Sum= " << sum << G4endl;
    129 
     129           << " Sum= " << sum << " q2= " << q2 << G4endl;
     130  }
    130131  sum *= material->GetElectronDensity() * q2 *  twopi_mc2_rcl2 /beta2;
    131132  return sum;
    … …  
    165166//   Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
    166167  SetupKinematics(p, mat, e);
    167   if(tau <= 0.0) return 0.0;
     168  if(tau <= 0.0) { return 0.0; }
    168169
    169170  G4double Barkas = BarkasCorrection (p, mat, e);
    … …  
    180181  sum *= material->GetElectronDensity() * q2 *  twopi_mc2_rcl2 /beta2;
    181182
    182   if(verbose > 1) G4cout << " Sum= " << sum << G4endl;
     183  if(verbose > 1) { G4cout << " Sum= " << sum << G4endl; }
    183184  return sum;
    184185}
    … …  
    219220    sum = ComputeIonCorrections(p,couple->GetMaterial(),e) - rest/e;
    220221
    221     if(verbose > 1) G4cout << " Sum= " << sum << " dSum= " << rest/e << G4endl;
     222    if(verbose > 1) { G4cout << " Sum= " << sum << " dSum= " << rest/e << G4endl; }
    222223  }
    223224  return sum;
    … …  
    267268    }
    268269    G4double e0= 13.6*eV*Z2;
    269     term += f*atomDensity[i]*KShell(shells.GetBindingEnergy(iz,0)/e0,ba2/Z2)/Z;
     270    term += f*atomDensity[i]*KShell(G4AtomicShells::GetBindingEnergy(iz,0)/e0,ba2/Z2)/Z;
    270271  }
    271272
    … …  
    293294      G4double e0= 13.6*eV*Z2*0.25;
    294295      G4double f = 0.125;
    295       G4int nmax = std::min(4,shells.GetNumberOfShells(iz));
     296      G4int nmax = std::min(4,G4AtomicShells::GetNumberOfShells(iz));
    296297      for(G4int j=1; j<nmax; j++) {
    297         G4double ne = G4double(shells.GetNumberOfElectrons(iz,j));
    298         G4double e1 = shells.GetBindingEnergy(iz,j);
     298        G4double ne = G4double(G4AtomicShells::GetNumberOfElectrons(iz,j));
     299        G4double e1 = G4AtomicShells::GetBindingEnergy(iz,j);
    299300        //   G4cout << "LShell: j= " << j << " ne= " << ne << " e(eV)= " << e/eV
    300301        //  << " e0(eV)= " << e0/eV << G4endl;
    … …  
    417418  for (G4int i = 0; i<numberOfElements; i++) {
    418419
     420    G4double res = 0.0;
    419421    G4double Z = (*theElementVector)[i]->GetZ();
    420422    G4int   iz = G4int(Z);
    … …  
    426428    }
    427429    G4double e0= 13.6*eV*Z2;
    428     term += f*atomDensity[i]*KShell(shells.GetBindingEnergy(iz,0)/e0,ba2/Z2)/Z;
     430    res += f*KShell(G4AtomicShells::GetBindingEnergy(iz,0)/e0,ba2/Z2);
    429431    if(2 < iz) {
    430432      G4double Zeff = Z - ZD[10];
    … …  
    434436      f = 0.125;
    435437      G4double eta = ba2/Z2;
    436       G4int ntot = shells.GetNumberOfShells(iz);
     438      G4int ntot = G4AtomicShells::GetNumberOfShells(iz);
    437439      G4int nmax = std::min(4, ntot);
    438440      G4double norm   = 0.0;
    439441      G4double eshell = 0.0;
    440       for(G4int j=1; j<nmax; j++) {
    441         G4double x = G4double(shells.GetNumberOfElectrons(iz,j));
    442         G4double e1 = shells.GetBindingEnergy(iz,j);
    443         norm   += x;
    444         eshell += e1*x;
    445         term += f*x*atomDensity[i]*LShell(e1/e0,eta)/Z;
     442      for(G4int j=1; j<nmax; ++j) {
     443        G4double e1 = G4AtomicShells::GetBindingEnergy(iz,j);
     444        G4int ne = G4AtomicShells::GetNumberOfElectrons(iz,j);
     445        norm   += ne;
     446        eshell += e1*ne;
     447        res += f*ne*LShell(e1/e0,eta);
    446448      }
    447       if(10 < iz) {
     449      if(ntot > nmax) {
    448450        eshell /= norm;
    449451        G4double eeff = eshell*eta;
    450         for(G4int k=nmax; k<ntot; k++) {
    451           G4double x = G4double(shells.GetNumberOfElectrons(iz,k));
    452           G4double e1 = shells.GetBindingEnergy(iz,k);
    453           term += f*x*atomDensity[i]*LShell(e1/e0,eeff/e1)/Z;
    454           //          term += f*x*atomDensity[i]*LShell(eshell/e0,eeff/e1)/Z;
     452        for(G4int k=nmax; k<ntot; ++k) {
     453          G4int ne = G4AtomicShells::GetNumberOfElectrons(iz,k);
     454          G4double e1 = G4AtomicShells::GetBindingEnergy(iz,k);
     455          res += f*ne*LShell(e1/e0,eeff/e1);
     456          //res += f*ne*LShell(e1/eshell,eeff/e1);
    455457        }
     458      }
    456459        /*
    457460        if(28 >= iz) {
    458           term += f*(Z - 10.)*atomDensity[i]*LShell(eshell,HM[iz-11]*eta)/Z;
     461          res += f*(Z - 10.)*LShell(eshell,HM[iz-11]*eta);
    459462        } else if(32 >= iz) {
    460           term += f*18.0*atomDensity[i]*LShell(eshell,HM[iz-11]*eta)/Z;
     463          res += f*18.0*LShell(eshell,HM[iz-11]*eta);
    461464        } else if(60 >= iz) {
    462           term += f*18.0*atomDensity[i]*LShell(eshell,HM[iz-11]*eta)/Z;
    463           term += f*(Z - 28.)*atomDensity[i]*LShell(eshell,HN[iz-33]*eta)/Z;
     465          res += f*18.0*LShell(eshell,HM[iz-11]*eta);
     466          res += f*(Z - 28.)*LShell(eshell,HN[iz-33]*eta);
    464467        } else {
    465           term += f*18.0*atomDensity[i]*LShell(eshell,HM[53]*eta)/Z;
    466           term += f*32.0*atomDensity[i]*LShell(eshell,HN[30]*eta)/Z;
    467           term += f*(Z - 60.)*atomDensity[i]*LShell(eshell,150.*eta)/Z;
     468          res += f*18.0*LShell(eshell,HM[50]*eta);
     469          res += f*32.0*LShell(eshell,HN[30]*eta);
     470          res += f*(Z - 60.)*LShell(eshell,150.*eta);
    468471        }
    469472        */
    470       }
    471     }
    472     //term += atomDensity[i]*MSH[iz]/(ba2*ba2);
     473    }
     474    //term += atomDensity[i]*(res/Z + MSH[iz]/(ba2*ba2));
     475    term += res*atomDensity[i]/Z;
    473476  }
    474477
    475478  term /= material->GetTotNbOfAtomsPerVolume();
     479  //if(charge < 0.0) { term = -term; }
    476480  return term;
    477481}
    … …  
    552556
    553557  BarkasTerm *= 1.29*charge/material->GetTotNbOfAtomsPerVolume();
     558
     559  // temporary protection
     560  if(charge < -7.0 ) { BarkasTerm *= (-7.0/charge); }
     561
    554562  return BarkasTerm;
    555563}
    … …  
    574582  } while (del > 0.01*term);
    575583
    576   return -y2*term;
     584  G4double res = -y2*term;
     585  // temporary protection
     586  if(q2 > 49. && res < -0.2) { res = -0.2; }
     587
     588  return res;
    577589}
    578590
    … …  
    808820  ionList[idx]  = ion;
    809821  stopData[idx] = vv;
    810   if(verbose>1) G4cout << "End data set " << G4endl;
     822  if(verbose>1) { G4cout << "End data set " << G4endl; }
    811823}
    812824
    … …  
    13951407    18.2
    13961408  };
    1397   for(i=0; i<53; i++) {HM[i] = hm[i];}
    1398   for(i=0; i<31; i++) {HN[i] = hn[i];}
     1409  for(i=0; i<53; ++i) {HM[i] = hm[i];}
     1410  for(i=0; i<31; ++i) {HN[i] = hn[i];}
    13991411
    14001412  const G4double mm[93] = {

  • trunk/source/processes/electromagnetic/utils/src/G4EmElementSelector.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmElementSelector.cc,v 1.11 2009/09/29 11:31:37 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4EmElementSelector.cc,v 1.12 2010/04/27 16:59:52 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5858                                         G4double emin,
    5959                                         G4double emax,
    60                                          G4bool /*spline*/):
     60                                         G4bool spline):
    6161  model(mod), material(mat), nbins(bins), cutEnergy(-1.0),
    6262  lowEnergy(emin), highEnergy(emax)
    … …  
    6868  if(nElmMinusOne > 0) {
    6969    xSections.reserve(n);
    70     for(G4int i=0; i<n; ++i) {
    71       G4PhysicsLogVector* v = new G4PhysicsLogVector(lowEnergy,highEnergy,nbins);
    72       //v->SetSpline(spline);
     70    G4PhysicsLogVector* v0 = new G4PhysicsLogVector(lowEnergy,highEnergy,nbins);
     71    xSections.push_back(v0);
     72    v0->SetSpline(spline);
     73    for(G4int i=1; i<n; ++i) {
     74      G4PhysicsLogVector* v = new G4PhysicsLogVector(*v0);
     75      v->SetSpline(spline);
    7376      xSections.push_back(v);
    7477    }
    … …  
    102105  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
    103106
    104   G4int i;
    105 
    106107  // loop over bins
    107108  for(G4int j=0; j<=nbins; ++j) {
    … …  
    110111    cross = 0.0;
    111112    //G4cout << "j= " << j << " e(MeV)= " << e/MeV << G4endl;
    112     for (i=0; i<=nElmMinusOne; ++i) {
     113    for (G4int i=0; i<=nElmMinusOne; ++i) {
    113114      cross += theAtomNumDensityVector[i]*     
    114115        model->ComputeCrossSectionPerAtom(part, (*theElementVector)[i], e,
    … …  
    119120
    120121  // xSections start from null, so use probabilities from the next bin
    121   if(DBL_MIN >= (*xSections[nElmMinusOne])[0]) {
    122     for (i=0; i<=nElmMinusOne; ++i) {
     122  if(0.0 == (*xSections[nElmMinusOne])[0]) {
     123    for (G4int i=0; i<=nElmMinusOne; ++i) {
    123124      xSections[i]->PutValue(0, (*xSections[i])[1]);
    124125    }
    125126  }
    126127  // xSections ends with null, so use probabilities from the previous bin
    127   if(DBL_MIN >= (*xSections[nElmMinusOne])[nbins]) {
    128     for (i=0; i<=nElmMinusOne; ++i) {
     128  if(0.0 == (*xSections[nElmMinusOne])[nbins]) {
     129    for (G4int i=0; i<=nElmMinusOne; ++i) {
    129130      xSections[i]->PutValue(nbins, (*xSections[i])[nbins-1]);
    130131    }
    … …  
    134135    cross = (*xSections[nElmMinusOne])[j];
    135136    // only for positive X-section
    136     if(cross > DBL_MIN) {
    137       for (i=0; i<nElmMinusOne; ++i) {
     137    if(cross > 0.0) {
     138      for (G4int i=0; i<nElmMinusOne; ++i) {
    138139        G4double x = (*xSections[i])[j]/cross;
    139140        xSections[i]->PutValue(j, x);

  • trunk/source/processes/electromagnetic/utils/src/G4EmModelManager.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmModelManager.cc,v 1.58 2009/10/29 18:07:08 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4EmModelManager.cc,v 1.60 2010/04/12 18:28:40 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    420420  size_t idx = 1;
    421421  if(secondaryParticle) {
    422     if( secondaryParticle == G4Gamma::Gamma() )           idx = 0;
    423     else if( secondaryParticle == G4Positron::Positron()) idx = 2;
    424   }
    425 
    426   if(numOfCouples > theCuts.size()) {theCuts.resize(numOfCouples);}
     422    if( secondaryParticle == G4Gamma::Gamma() )           { idx = 0; }
     423    else if( secondaryParticle == G4Positron::Positron()) { idx = 2; }
     424  }
     425
     426  if(numOfCouples > theCuts.size()) { theCuts.resize(numOfCouples); }
    427427  if(minSubRange < 1.0 && numOfCouples > theSubCuts.size()) {
    428428    theSubCuts.resize(numOfCouples);
    … …  
    460460        G4double tcutmax =
    461461          theCoupleTable->ConvertRangeToEnergy(secondaryParticle,material,rcut);
    462         if(tcutmax < subcut) subcut = tcutmax;
    463       }
    464     }
    465 
     462        if(tcutmax < subcut) { subcut = tcutmax; }
     463      }
     464    }
     465    /*
    466466    G4int nm = setOfRegionModels[reg]->NumberOfModels();
    467467    for(G4int j=0; j<nm; ++j) {
    468468
    469469      G4VEmModel* model = models[setOfRegionModels[reg]->ModelIndex(j)];
    470 
    471470      G4double tcutmin = model->MinEnergyCut(particle, couple);
    472471
    473       if(cut < tcutmin) cut = tcutmin;
    474       if(subcut < tcutmin) subcut = tcutmin;
     472      if(cut < tcutmin)    { cut = tcutmin; }
     473      if(subcut < tcutmin) { subcut = tcutmin; }
    475474      if(1 < verboseLevel) {
    476475            G4cout << "The model # " << j
    … …  
    482481      }
    483482    }
     483    */
    484484    theCuts[i] = cut;
    485     if(minSubRange < 1.0) theSubCuts[i] = subcut;
     485    if(minSubRange < 1.0) { theSubCuts[i] = subcut; }
    486486  }
    487487
    … …  
    497497    }
    498498  }
    499   if(1 == nn) severalModels = false;
     499  if(1 == nn) { severalModels = false; }
    500500
    501501  if(1 < verboseLevel) {

  • trunk/source/processes/electromagnetic/utils/src/G4EmProcessOptions.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmProcessOptions.cc,v 1.27 2009/10/29 19:25:28 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4EmProcessOptions.cc,v 1.28 2010/04/27 16:59:52 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    7979{
    8080  theManager->SetLossFluctuations(val);
     81}
     82
     83//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     84
     85void G4EmProcessOptions::SetSubCutoff(G4bool val, const G4Region* r)
     86{
     87  theManager->SetSubCutoff(val, r);
     88}
     89
     90//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     91
     92void G4EmProcessOptions::SetIntegral(G4bool val)
     93{
     94  theManager->SetIntegral(val);
     95}
     96
     97//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     98
     99void G4EmProcessOptions::SetMinSubRange(G4double val)
     100{
     101  theManager->SetMinSubRange(val);
     102}
     103
     104//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     105
     106void G4EmProcessOptions::SetMinEnergy(G4double val)
     107{
     108  theManager->SetMinEnergy(val);
     109}
     110
     111//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     112
     113void G4EmProcessOptions::SetMaxEnergy(G4double val)
     114{
     115  theManager->SetMaxEnergy(val);
     116}
     117
     118//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     119
     120void G4EmProcessOptions::SetMaxEnergyForCSDARange(G4double val)
     121{
     122  theManager->SetMaxEnergyForCSDARange(val);
     123}
     124
     125//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     126
     127void G4EmProcessOptions::SetMaxEnergyForMuons(G4double val)
     128{
     129  theManager->SetMaxEnergyForMuons(val);
     130}
     131
     132//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     133
     134void G4EmProcessOptions::SetDEDXBinning(G4int val)
     135{
     136  theManager->SetDEDXBinning(val);
     137}
     138
     139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     140
     141void G4EmProcessOptions::SetDEDXBinningForCSDARange(G4int val)
     142{
     143  theManager->SetDEDXBinningForCSDARange(val);
     144}
     145
     146//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     147
     148void G4EmProcessOptions::SetLambdaBinning(G4int val)
     149{
     150  theManager->SetLambdaBinning(val);
     151}
     152
     153//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     154
     155void G4EmProcessOptions::SetStepFunction(G4double v1, G4double v2)
     156{
     157  theManager->SetStepFunction(v1, v2);
     158}
     159
     160//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     161
     162void G4EmProcessOptions::SetRandomStep(G4bool val)
     163{
     164  theManager->SetRandomStep(val);
     165}
     166
     167//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     168
     169void G4EmProcessOptions::SetApplyCuts(G4bool val)
     170{
     171  const std::vector<G4VEmProcess*>& w =
     172        theManager->GetEmProcessVector();
     173  std::vector<G4VEmProcess*>::const_iterator itp;
     174  for(itp = w.begin(); itp != w.end(); itp++) {
     175    G4VEmProcess* q = *itp;
     176    if(q) { q->SetApplyCuts(val); }
     177  }
     178}
     179
     180//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     181
     182void G4EmProcessOptions::SetBuildCSDARange(G4bool val)
     183{
     184  theManager->SetBuildCSDARange(val);
     185}
     186
     187//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     188
     189void G4EmProcessOptions::SetVerbose(G4int val, const G4String& name)
     190{
     191  G4bool all = false;
     192  if("all" == name) { all = true; }
     193  const std::vector<G4VEnergyLossProcess*>& v =
     194        theManager->GetEnergyLossProcessVector();
     195
     196  if(all) {
     197    theManager->SetVerbose(val);
     198    return;
     199  }
     200
     201  std::vector<G4VEnergyLossProcess*>::const_iterator itr;
     202  for(itr = v.begin(); itr != v.end(); ++itr) {
     203    G4VEnergyLossProcess* p = *itr;
     204    if(p) {
     205      if (p->GetProcessName() == name) { p->SetVerboseLevel(val); }
     206    }
     207  }
     208  const std::vector<G4VEmProcess*>& w =
     209        theManager->GetEmProcessVector();
     210  std::vector<G4VEmProcess*>::const_iterator itp;
     211  for(itp = w.begin(); itp != w.end(); itp++) {
     212    G4VEmProcess* q = *itp;
     213    if(q) {
     214      if (q->GetProcessName() == name) { q->SetVerboseLevel(val); }
     215    }
     216  }
     217  const std::vector<G4VMultipleScattering*>& u =
     218        theManager->GetMultipleScatteringVector();
     219  std::vector<G4VMultipleScattering*>::const_iterator itm;
     220  for(itm = u.begin(); itm != u.end(); itm++) {
     221    G4VMultipleScattering* s = *itm;
     222    if(s) {
     223      if (s->GetProcessName() == name) { s->SetVerboseLevel(val); }
     224    }
     225  }
     226}
     227
     228//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     229
     230void G4EmProcessOptions::SetLambdaFactor(G4double val)
     231{
    81232  const std::vector<G4VEnergyLossProcess*>& v =
    82233        theManager->GetEnergyLossProcessVector();
    … …  
    84235  for(itr = v.begin(); itr != v.end(); itr++) {
    85236    G4VEnergyLossProcess* p = *itr;
    86     if(p) p->SetLossFluctuations(val);
    87   }
    88 }
    89 
    90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    91 
    92 void G4EmProcessOptions::SetSubCutoff(G4bool val, const G4Region* r)
    93 {
    94   theManager->SetSubCutoff(val);
    95   const std::vector<G4VEnergyLossProcess*>& v =
    96         theManager->GetEnergyLossProcessVector();
    97   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    98   for(itr = v.begin(); itr != v.end(); itr++) {
    99     G4VEnergyLossProcess* p = *itr;
    100     if(p) p->ActivateSubCutoff(val, r);
    101   }
    102 }
    103 
    104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    105 
    106 void G4EmProcessOptions::SetIntegral(G4bool val)
    107 {
    108   theManager->SetIntegral(val);
    109   const std::vector<G4VEnergyLossProcess*>& v =
    110         theManager->GetEnergyLossProcessVector();
    111   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    112   for(itr = v.begin(); itr != v.end(); itr++) {
    113     G4VEnergyLossProcess* p = *itr;
    114     if(p) p->SetIntegral(val);
    115   }
    116 }
    117 
    118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    119 
    120 void G4EmProcessOptions::SetMinSubRange(G4double val)
    121 {
    122   theManager->SetMinSubRange(val);
    123   const std::vector<G4VEnergyLossProcess*>& v =
    124         theManager->GetEnergyLossProcessVector();
    125   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    126   for(itr = v.begin(); itr != v.end(); itr++) {
    127     G4VEnergyLossProcess* p = *itr;
    128     if(p) p->SetMinSubRange(val);
    129   }
    130 }
    131 
    132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    133 
    134 void G4EmProcessOptions::SetMinEnergy(G4double val)
    135 {
    136   theManager->SetMinEnergy(val);
    137   const std::vector<G4VEnergyLossProcess*>& v =
    138         theManager->GetEnergyLossProcessVector();
    139   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    140   for(itr = v.begin(); itr != v.end(); itr++) {
    141     G4VEnergyLossProcess* p = *itr;
    142     if(p) p->SetMinKinEnergy(val);
    143   }
    144   const std::vector<G4VEmProcess*>& w =
    145         theManager->GetEmProcessVector();
    146   std::vector<G4VEmProcess*>::const_iterator itp;
    147   for(itp = w.begin(); itp != w.end(); itp++) {
    148     G4VEmProcess* q = *itp;
    149     if(q) q->SetMinKinEnergy(val);
    150   }
    151   const std::vector<G4VMultipleScattering*>& u =
    152         theManager->GetMultipleScatteringVector();
    153   std::vector<G4VMultipleScattering*>::const_iterator itm;
    154   for(itm = u.begin(); itm != u.end(); itm++) {
    155     G4VMultipleScattering* s = *itm;
    156     if(s) s->SetMinKinEnergy(val);
    157   }
    158 }
    159 
    160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    161 
    162 void G4EmProcessOptions::SetMaxEnergy(G4double val)
    163 {
    164   theManager->SetMaxEnergy(val);
    165   const std::vector<G4VEnergyLossProcess*>& v =
    166         theManager->GetEnergyLossProcessVector();
    167   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    168   for(itr = v.begin(); itr != v.end(); itr++) {
    169     G4VEnergyLossProcess* p = *itr;
    170     if(p) p->SetMaxKinEnergy(val);
    171   }
    172   const std::vector<G4VEmProcess*>& w =
    173         theManager->GetEmProcessVector();
    174   std::vector<G4VEmProcess*>::const_iterator itp;
    175   for(itp = w.begin(); itp != w.end(); itp++) {
    176     G4VEmProcess* q = *itp;
    177     if(q) q->SetMaxKinEnergy(val);
    178   }
    179   const std::vector<G4VMultipleScattering*>& u =
    180         theManager->GetMultipleScatteringVector();
    181   std::vector<G4VMultipleScattering*>::const_iterator itm;
    182   for(itm = u.begin(); itm != u.end(); itm++) {
    183     G4VMultipleScattering* s = *itm;
    184     if(s) s->SetMaxKinEnergy(val);
    185   }
    186 }
    187 
    188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    189 
    190 void G4EmProcessOptions::SetMaxEnergyForCSDARange(G4double val)
    191 {
    192   theManager->SetMaxEnergyForCSDARange(val);
    193   const std::vector<G4VEnergyLossProcess*>& v =
    194         theManager->GetEnergyLossProcessVector();
    195   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    196   for(itr = v.begin(); itr != v.end(); itr++) {
    197     G4VEnergyLossProcess* p = *itr;
    198     if(p) p->SetMaxKinEnergyForCSDARange(val);
    199   }
    200 }
    201 
    202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    203 
    204 void G4EmProcessOptions::SetMaxEnergyForMuons(G4double val)
    205 {
    206   theManager->SetMaxEnergyForMuons(val);
     237    if(p) { p->SetLambdaFactor(val); }
     238  }
     239}
     240
     241//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     242
     243void G4EmProcessOptions::ActivateDeexcitation(const G4String& pname,
     244                                              G4bool val,
     245                                              const G4String& reg)
     246{
     247  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     248  const G4Region* r = 0;
     249  if(reg == "" || reg == "World") {
     250    r = regionStore->GetRegion("DefaultRegionForTheWorld", false);
     251  } else {
     252    r = regionStore->GetRegion(reg, false);
     253  }
     254  if(!r) {
     255    G4cout << "G4EmProcessOptions::ActivateDeexcitation ERROR: G4Region <"
     256           << reg << "> not found, the command ignored" << G4endl;
     257    return;
     258  }
     259
    207260  const std::vector<G4VEnergyLossProcess*>& v =
    208261        theManager->GetEnergyLossProcessVector();
    … …  
    211264    G4VEnergyLossProcess* p = *itr;
    212265    if(p) {
    213       if(std::abs(p->Particle()->GetPDGMass() - 105.66*MeV) < MeV)
    214         p->SetMaxKinEnergy(val);
     266      if(pname == p->GetProcessName()) { p->ActivateDeexcitation(val,r); }
    215267    }
    216268  }
    … …  
    221273    G4VEmProcess* q = *itp;
    222274    if(q) {
    223       if(std::abs(q->Particle()->GetPDGMass() - 105.66*MeV) < MeV)
    224         q->SetMaxKinEnergy(val);
    225     }
    226   }
    227   /*
    228   const std::vector<G4VMultipleScattering*>& u =
    229         theManager->GetMultipleScatteringVector();
    230   std::vector<G4VMultipleScattering*>::const_iterator itm;
    231   for(itm = u.begin(); itm != u.end(); itm++) {
    232     G4VMultipleScattering* s = *itm;
    233     if(s) {
    234       if(std::abs(s->Particle()->GetPDGMass() - 105.66*MeV) < MeV)
    235         s->SetMaxKinEnergy(val);
    236     }
    237   }
    238   */
    239 }
    240 
    241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    242 
    243 void G4EmProcessOptions::SetDEDXBinning(G4int val)
    244 {
    245   theManager->SetDEDXBinning(val);
    246   const std::vector<G4VEnergyLossProcess*>& v =
    247         theManager->GetEnergyLossProcessVector();
    248   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    249   for(itr = v.begin(); itr != v.end(); itr++) {
    250     G4VEnergyLossProcess* p = *itr;
    251     if(p) p->SetDEDXBinning(val);
    252   }
    253 }
    254 
    255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    256 
    257 void G4EmProcessOptions::SetDEDXBinningForCSDARange(G4int val)
    258 {
    259   theManager->SetDEDXBinningForCSDARange(val);
    260   const std::vector<G4VEnergyLossProcess*>& v =
    261         theManager->GetEnergyLossProcessVector();
    262   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    263   for(itr = v.begin(); itr != v.end(); itr++) {
    264     G4VEnergyLossProcess* p = *itr;
    265     if(p) p->SetDEDXBinningForCSDARange(val);
    266   }
    267 }
    268 
    269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    270 
    271 void G4EmProcessOptions::SetLambdaBinning(G4int val)
    272 {
    273   theManager->SetLambdaBinning(val);
     275      if(pname == q->GetProcessName()) { q->ActivateDeexcitation(val,r); }
     276    }
     277  }
     278}
     279
     280//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     281
     282void G4EmProcessOptions::SetMscStepLimitation(G4MscStepLimitType val)
     283{
     284  const std::vector<G4VMultipleScattering*>& u =
     285        theManager->GetMultipleScatteringVector();
     286  std::vector<G4VMultipleScattering*>::const_iterator itm;
     287  for(itm = u.begin(); itm != u.end(); itm++) {
     288    if(*itm) (*itm)->SetStepLimitType(val);
     289  }
     290}
     291
     292//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     293
     294void G4EmProcessOptions::SetMscLateralDisplacement(G4bool val)
     295{
     296  const std::vector<G4VMultipleScattering*>& u =
     297        theManager->GetMultipleScatteringVector();
     298  std::vector<G4VMultipleScattering*>::const_iterator itm;
     299  for(itm = u.begin(); itm != u.end(); itm++) {
     300    if(*itm) { (*itm)->SetLateralDisplasmentFlag(val); }
     301  }
     302}
     303
     304//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     305
     306void G4EmProcessOptions::SetSkin(G4double val)
     307{
     308  if(val < 0.0) return;
     309  const std::vector<G4VMultipleScattering*>& u =
     310        theManager->GetMultipleScatteringVector();
     311  std::vector<G4VMultipleScattering*>::const_iterator itm;
     312  for(itm = u.begin(); itm != u.end(); itm++) {
     313    if(*itm) { (*itm)->SetSkin(val); }
     314  }
     315}
     316
     317//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     318
     319void G4EmProcessOptions::SetMscRangeFactor(G4double val)
     320{
     321  if(val < 0.0) return;
     322  const std::vector<G4VMultipleScattering*>& u =
     323        theManager->GetMultipleScatteringVector();
     324  std::vector<G4VMultipleScattering*>::const_iterator itm;
     325  for(itm = u.begin(); itm != u.end(); itm++) {
     326    if(*itm) { (*itm)->SetRangeFactor(val); }
     327  }
     328}
     329
     330//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     331
     332void G4EmProcessOptions::SetMscGeomFactor(G4double val)
     333{
     334  if(val < 0.0) { return; }
     335  const std::vector<G4VMultipleScattering*>& u =
     336        theManager->GetMultipleScatteringVector();
     337  std::vector<G4VMultipleScattering*>::const_iterator itm;
     338  for(itm = u.begin(); itm != u.end(); itm++) {
     339    if(*itm) { (*itm)->SetGeomFactor(val); }
     340  }
     341}
     342
     343//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     344
     345void G4EmProcessOptions::SetPolarAngleLimit(G4double val)
     346{
     347  const std::vector<G4VMultipleScattering*>& u =
     348        theManager->GetMultipleScatteringVector();
     349  std::vector<G4VMultipleScattering*>::const_iterator itm;
     350  for(itm = u.begin(); itm != u.end(); itm++) {
     351    if(*itm) { (*itm)->SetPolarAngleLimit(val); }
     352  }
    274353  const std::vector<G4VEmProcess*>& w =
    275354        theManager->GetEmProcessVector();
    … …  
    277356  for(itp = w.begin(); itp != w.end(); itp++) {
    278357    G4VEmProcess* q = *itp;
    279     if(q) q->SetLambdaBinning(val);
    280   }
    281   const std::vector<G4VMultipleScattering*>& u =
    282         theManager->GetMultipleScatteringVector();
    283   std::vector<G4VMultipleScattering*>::const_iterator itm;
    284   for(itm = u.begin(); itm != u.end(); itm++) {
    285     G4VMultipleScattering* s = *itm;
    286     if(s) s->SetBinning(val);
    287   }
    288 }
    289 
    290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    291 
    292 void G4EmProcessOptions::SetStepFunction(G4double v1, G4double v2)
    293 {
    294   theManager->SetStepFunction(v1, v2);
    295   const std::vector<G4VEnergyLossProcess*>& v =
    296         theManager->GetEnergyLossProcessVector();
    297   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    298   for(itr = v.begin(); itr != v.end(); itr++) {
    299     G4VEnergyLossProcess* p = *itr;
    300     if(p) p->SetStepFunction(v1, v2);
    301   }
    302 }
    303 
    304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    305 
    306 void G4EmProcessOptions::SetRandomStep(G4bool val)
    307 {
    308   theManager->SetRandomStep(val);
    309   const std::vector<G4VEnergyLossProcess*>& v =
    310         theManager->GetEnergyLossProcessVector();
    311   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    312   for(itr = v.begin(); itr != v.end(); itr++) {
    313     G4VEnergyLossProcess* p = *itr;
    314     if(p) p->SetRandomStep(val);
    315   }
    316 }
    317 
    318 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    319 
    320 void G4EmProcessOptions::SetApplyCuts(G4bool val)
    321 {
    322   const std::vector<G4VEmProcess*>& w =
    323         theManager->GetEmProcessVector();
    324   std::vector<G4VEmProcess*>::const_iterator itp;
    325   for(itp = w.begin(); itp != w.end(); itp++) {
    326     G4VEmProcess* q = *itp;
    327     if(q) q->SetApplyCuts(val);
    328   }
    329 }
    330 
    331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    332 
    333 void G4EmProcessOptions::SetBuildCSDARange(G4bool val)
    334 {
    335   theManager->SetBuildCSDARange(val);
    336 }
    337 
    338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    339 
    340 void G4EmProcessOptions::SetVerbose(G4int val, const G4String& name)
    341 {
    342   G4bool all = false;
    343   if("all" == name) all = true;
    344   const std::vector<G4VEnergyLossProcess*>& v =
    345         theManager->GetEnergyLossProcessVector();
    346 
    347   if(all) theManager->SetVerbose(val);
    348 
    349   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    350   for(itr = v.begin(); itr != v.end(); itr++) {
    351     G4VEnergyLossProcess* p = *itr;
    352     if(p) {
    353       if(all) p->SetVerboseLevel(val);
    354       else if (p->GetProcessName() == name) p->SetVerboseLevel(val);
    355     }
    356   }
    357   const std::vector<G4VEmProcess*>& w =
    358         theManager->GetEmProcessVector();
    359   std::vector<G4VEmProcess*>::const_iterator itp;
    360   for(itp = w.begin(); itp != w.end(); itp++) {
    361     G4VEmProcess* q = *itp;
    362     if(q) {
    363       if(all) q->SetVerboseLevel(val);
    364       else if (q->GetProcessName() == name) q->SetVerboseLevel(val);
    365     }
    366   }
    367   const std::vector<G4VMultipleScattering*>& u =
    368         theManager->GetMultipleScatteringVector();
    369   std::vector<G4VMultipleScattering*>::const_iterator itm;
    370   for(itm = u.begin(); itm != u.end(); itm++) {
    371     G4VMultipleScattering* s = *itm;
    372     if(s) {
    373       if(all) s->SetVerboseLevel(val);
    374       else if (s->GetProcessName() == name) s->SetVerboseLevel(val);
    375     }
    376   }
    377 }
    378 
    379 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    380 
    381 void G4EmProcessOptions::SetLambdaFactor(G4double val)
    382 {
    383   const std::vector<G4VEnergyLossProcess*>& v =
    384         theManager->GetEnergyLossProcessVector();
    385   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    386   for(itr = v.begin(); itr != v.end(); itr++) {
    387     G4VEnergyLossProcess* p = *itr;
    388     if(p) p->SetLambdaFactor(val);
    389   }
    390 
    391 }
    392 
    393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    394 
    395 void G4EmProcessOptions::ActivateDeexcitation(const G4String& pname,
    396                                               G4bool val,
    397                                               const G4String& reg)
    398 {
    399   G4RegionStore* regionStore = G4RegionStore::GetInstance();
    400   const G4Region* r = 0;
    401   if(reg == "" || reg == "World") {
    402     r = regionStore->GetRegion("DefaultRegionForTheWorld", false);
    403   } else {
    404     r = regionStore->GetRegion(reg, false);
    405   }
    406   if(!r) {
    407     G4cout << "G4EmProcessOptions::ActivateDeexcitation ERROR: G4Region <"
    408            << reg << "> not found, the command ignored" << G4endl;
    409     return;
    410   }
    411 
    412   const std::vector<G4VEnergyLossProcess*>& v =
    413         theManager->GetEnergyLossProcessVector();
    414   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    415   for(itr = v.begin(); itr != v.end(); itr++) {
    416     G4VEnergyLossProcess* p = *itr;
    417     if(p) {
    418       if(pname == p->GetProcessName()) p->ActivateDeexcitation(val,r);
    419     }
    420   }
    421   const std::vector<G4VEmProcess*>& w =
    422         theManager->GetEmProcessVector();
    423   std::vector<G4VEmProcess*>::const_iterator itp;
    424   for(itp = w.begin(); itp != w.end(); itp++) {
    425     G4VEmProcess* q = *itp;
    426     if(q) {
    427       if(pname == q->GetProcessName()) q->ActivateDeexcitation(val,r);
    428     }
    429   }
    430 }
    431 
    432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    433 
    434 void G4EmProcessOptions::SetMscStepLimitation(G4MscStepLimitType val)
    435 {
    436   const std::vector<G4VMultipleScattering*>& u =
    437         theManager->GetMultipleScatteringVector();
    438   std::vector<G4VMultipleScattering*>::const_iterator itm;
    439   for(itm = u.begin(); itm != u.end(); itm++) {
    440     if(*itm) (*itm)->SetStepLimitType(val);
    441   }
    442 }
    443 
    444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    445 
    446 void G4EmProcessOptions::SetMscLateralDisplacement(G4bool val)
    447 {
    448   const std::vector<G4VMultipleScattering*>& u =
    449         theManager->GetMultipleScatteringVector();
    450   std::vector<G4VMultipleScattering*>::const_iterator itm;
    451   for(itm = u.begin(); itm != u.end(); itm++) {
    452     if(*itm) (*itm)->SetLateralDisplasmentFlag(val);
    453   }
    454 }
    455 
    456 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    457 
    458 void G4EmProcessOptions::SetSkin(G4double val)
    459 {
    460   if(val < 0.0) return;
    461   const std::vector<G4VMultipleScattering*>& u =
    462         theManager->GetMultipleScatteringVector();
    463   std::vector<G4VMultipleScattering*>::const_iterator itm;
    464   for(itm = u.begin(); itm != u.end(); itm++) {
    465     if(*itm) {
    466       (*itm)->SetSkin(val);
    467     }
    468   }
    469 }
    470 
    471 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    472 
    473 void G4EmProcessOptions::SetMscRangeFactor(G4double val)
    474 {
    475   if(val < 0.0) return;
    476   const std::vector<G4VMultipleScattering*>& u =
    477         theManager->GetMultipleScatteringVector();
    478   std::vector<G4VMultipleScattering*>::const_iterator itm;
    479   for(itm = u.begin(); itm != u.end(); itm++) {
    480     if(*itm) (*itm)->SetRangeFactor(val);
    481   }
    482 }
    483 
    484 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    485 
    486 void G4EmProcessOptions::SetMscGeomFactor(G4double val)
    487 {
    488   if(val < 0.0) return;
    489   const std::vector<G4VMultipleScattering*>& u =
    490         theManager->GetMultipleScatteringVector();
    491   std::vector<G4VMultipleScattering*>::const_iterator itm;
    492   for(itm = u.begin(); itm != u.end(); itm++) {
    493     if(*itm) (*itm)->SetGeomFactor(val);
    494   }
    495 }
    496 
    497 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    498 
    499 void G4EmProcessOptions::SetPolarAngleLimit(G4double val)
    500 {
    501   const std::vector<G4VMultipleScattering*>& u =
    502         theManager->GetMultipleScatteringVector();
    503   std::vector<G4VMultipleScattering*>::const_iterator itm;
    504   for(itm = u.begin(); itm != u.end(); itm++) {
    505     if(*itm) (*itm)->SetPolarAngleLimit(val);
    506   }
    507   const std::vector<G4VEmProcess*>& w =
    508         theManager->GetEmProcessVector();
    509   std::vector<G4VEmProcess*>::const_iterator itp;
    510   for(itp = w.begin(); itp != w.end(); itp++) {
    511     G4VEmProcess* q = *itp;
    512     if(q) q->SetPolarAngleLimit(val);
     358    if(q) { q->SetPolarAngleLimit(val); }
    513359  }
    514360}

  • trunk/source/processes/electromagnetic/utils/src/G4LossTableManager.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LossTableManager.cc,v 1.97 2009/10/29 19:25:28 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4LossTableManager.cc,v 1.101 2010/06/04 15:33:56 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    7070// 12-02-07 Add SetSkin, SetLinearLossLimit (V.Ivanchenko)
    7171// 18-06-07 Move definition of msc parameters to G4EmProcessOptions (V.Ivanchenko)
    72 // 21-02-08 Add G4EmSaturation (V.Ivanchenko)
     72// 21-02-08 Added G4EmSaturation (V.Ivanchenko)
     73// 12-04-10 Added PreparePhsyicsTables and BuildPhysicsTables entries (V.Ivanchenko)
    7374//
    7475// Class Description:
    … …  
    9697#include "G4EmTableType.hh"
    9798#include "G4LossTableBuilder.hh"
     99#include "G4Region.hh"
    98100
    99101G4LossTableManager* G4LossTableManager::theInstance = 0;
    … …  
    114116G4LossTableManager::~G4LossTableManager()
    115117{
    116   for (G4int i=0; i<n_loss; i++) {
    117     if( loss_vector[i] ) delete loss_vector[i];
     118  for (G4int i=0; i<n_loss; ++i) {
     119    if( loss_vector[i] ) { delete loss_vector[i]; }
    118120  }
    119121  size_t msc = msc_vector.size();
    120   for (size_t j=0; j<msc; j++) {
    121     if( msc_vector[j] ) delete msc_vector[j];
     122  for (size_t j=0; j<msc; ++j) {
     123    if( msc_vector[j] ) { delete msc_vector[j]; }
    122124  }
    123125  size_t emp = emp_vector.size();
    124   for (size_t k=0; k<emp; k++) {
    125     if( emp_vector[k] ) delete emp_vector[k];
     126  for (size_t k=0; k<emp; ++k) {
     127    if( emp_vector[k] ) { delete emp_vector[k]; }
    126128  }
    127129  size_t mod = mod_vector.size();
    128   for (size_t a=0; a<mod; a++) {
    129     if( mod_vector[a] ) delete mod_vector[a];
     130  for (size_t a=0; a<mod; ++a) {
     131    if( mod_vector[a] ) { delete mod_vector[a]; }
    130132  }
    131133  size_t fmod = fmod_vector.size();
    132   for (size_t b=0; b<fmod; b++) {
    133     if( fmod_vector[b] ) delete fmod_vector[b];
     134  for (size_t b=0; b<fmod; ++b) {
     135    if( fmod_vector[b] ) { delete fmod_vector[b]; }
    134136  }
    135137  Clear();
    … …  
    146148  n_loss = 0;
    147149  run = 0;
    148   //  first_entry = true;
     150  startInitialisation = false;
    149151  all_tables_are_built = false;
    150   all_tables_are_stored = false;
    151152  currentLoss = 0;
    152153  currentParticle = 0;
     154  firstParticle = 0;
    153155  lossFluctuationFlag = true;
    154156  subCutoffFlag = false;
    … …  
    158160  maxFinalStep = 0.0;
    159161  minKinEnergy = 0.1*keV;
    160   maxKinEnergy = 100.0*TeV;
    161   maxKinEnergyForMuons = 100.*TeV;
    162   theMessenger = new G4EnergyLossMessenger();
    163   theElectron  = G4Electron::Electron();
    164   tableBuilder = new G4LossTableBuilder();
    165   emCorrections= new G4EmCorrections();
    166   emSaturation = new G4EmSaturation();
    167   emConfigurator = new G4EmConfigurator();
     162  maxKinEnergy = 10.0*TeV;
     163  nbinsLambda  = 77;
     164  nbinsPerDecade = 7;
     165  maxKinEnergyForMuons = 10.*TeV;
    168166  integral = true;
    169167  integralActive = false;
    … …  
    178176  factorForAngleLimit = 1.0;
    179177  verbose = 1;
     178  theMessenger = new G4EnergyLossMessenger();
     179  theElectron  = G4Electron::Electron();
     180  tableBuilder = new G4LossTableBuilder();
     181  emCorrections= new G4EmCorrections();
     182  emSaturation = new G4EmSaturation();
     183  emConfigurator = new G4EmConfigurator(verbose);
    180184  tableBuilder->SetSplineFlag(splineFlag);
    181185}
    … …  
    207211void G4LossTableManager::Register(G4VEnergyLossProcess* p)
    208212{
    209   n_loss++;
     213  ++n_loss;
    210214  loss_vector.push_back(p);
    211215  part_vector.push_back(0);
    … …  
    217221  isActive.push_back(true);
    218222  all_tables_are_built = false;
    219   if(!lossFluctuationFlag) p->SetLossFluctuations(false);
    220   if(subCutoffFlag)        p->ActivateSubCutoff(true);
    221   if(rndmStepFlag)         p->SetRandomStep(true);
    222   if(stepFunctionActive)   p->SetStepFunction(maxRangeVariation, maxFinalStep);
    223   if(integralActive)       p->SetIntegral(integral);
    224   if(minEnergyActive)      p->SetMinKinEnergy(minKinEnergy);
    225   if(maxEnergyActive)      p->SetMaxKinEnergy(maxKinEnergy);
     223  if(!lossFluctuationFlag) { p->SetLossFluctuations(false); }
     224  if(subCutoffFlag)        { p->ActivateSubCutoff(true); }
     225  if(rndmStepFlag)         { p->SetRandomStep(true); }
     226  if(stepFunctionActive)   { p->SetStepFunction(maxRangeVariation, maxFinalStep); }
     227  if(integralActive)       { p->SetIntegral(integral); }
     228  if(minEnergyActive)      { p->SetMinKinEnergy(minKinEnergy); }
     229  if(maxEnergyActive)      { p->SetMaxKinEnergy(maxKinEnergy); }
    226230  if(verbose > 1)
    227231    G4cout << "G4LossTableManager::Register G4VEnergyLossProcess : "
    … …  
    233237void G4LossTableManager::DeRegister(G4VEnergyLossProcess* p)
    234238{
    235   for (G4int i=0; i<n_loss; i++) {
    236     if(loss_vector[i] == p) loss_vector[i] = 0;
     239  for (G4int i=0; i<n_loss; ++i) {
     240    if(loss_vector[i] == p) { loss_vector[i] = 0; }
    237241  }
    238242}
    … …  
    254258{
    255259  size_t msc = msc_vector.size();
    256   for (size_t i=0; i<msc; i++) {
    257     if(msc_vector[i] == p) msc_vector[i] = 0;
     260  for (size_t i=0; i<msc; ++i) {
     261    if(msc_vector[i] == p) { msc_vector[i] = 0; }
    258262  }
    259263}
    … …  
    275279{
    276280  size_t emp = emp_vector.size();
    277   for (size_t i=0; i<emp; i++) {
    278     if(emp_vector[i] == p) emp_vector[i] = 0;
     281  for (size_t i=0; i<emp; ++i) {
     282    if(emp_vector[i] == p) { emp_vector[i] = 0; }
    279283  }
    280284}
    … …  
    296300{
    297301  size_t n = mod_vector.size();
    298   for (size_t i=0; i<n; i++) {
    299     if(mod_vector[i] == p) mod_vector[i] = 0;
     302  for (size_t i=0; i<n; ++i) {
     303    if(mod_vector[i] == p) { mod_vector[i] = 0; }
    300304  }
    301305}
    … …  
    317321{
    318322  size_t n = fmod_vector.size();
    319   for (size_t i=0; i<n; i++) {
    320     if(fmod_vector[i] == p) fmod_vector[i] = 0;
     323  for (size_t i=0; i<n; ++i) {
     324    if(fmod_vector[i] == p) { fmod_vector[i] = 0; }
    321325  }
    322326}
    … …  
    336340     G4VEnergyLossProcess* p)
    337341{
    338   n_loss++;
     342  ++n_loss;
    339343  loss_vector.push_back(p);
    340344  part_vector.push_back(part);
    … …  
    349353//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    350354
    351 void G4LossTableManager::EnergyLossProcessIsInitialised(
    352                    const G4ParticleDefinition* particle,
    353                    G4VEnergyLossProcess* p)
    354 {
    355   if (run == 0 || (particle == firstParticle && all_tables_are_built) ) {
     355void
     356G4LossTableManager::PreparePhysicsTable(const G4ParticleDefinition* particle,
     357                                        G4VEnergyLossProcess* p)
     358{
     359  if (1 < verbose) {
     360    G4cout << "G4LossTableManager::PreparePhysicsTable for "
     361           << particle->GetParticleName()
     362           << " and " << p->GetProcessName() << " run= " << run << G4endl;
     363  }
     364  // start initialisation for the first run
     365  startInitialisation = true;
     366
     367  if( 0 == run ) {
     368    emConfigurator->PrepareModels(particle, p);
     369
     370    // initialise particles for given process
     371    for (G4int j=0; j<n_loss; ++j) {
     372      if (p == loss_vector[j]) {
     373        if (!part_vector[j]) { part_vector[j] = particle; }
     374      }
     375    }
     376  }
     377}
     378
     379//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     380
     381void
     382G4LossTableManager::PreparePhysicsTable(const G4ParticleDefinition* particle,
     383                                        G4VEmProcess* p)
     384{
     385  if (1 < verbose) {
     386    G4cout << "G4LossTableManager::PreparePhysicsTable for "
     387           << particle->GetParticleName()
     388           << " and " << p->GetProcessName() << G4endl;
     389  }
     390  // start initialisation for the first run
     391  if( 0 == run ) {
     392    emConfigurator->PrepareModels(particle, p);
     393  }
     394  startInitialisation = true;
     395}
     396
     397//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     398
     399void
     400G4LossTableManager::PreparePhysicsTable(const G4ParticleDefinition* particle,
     401                                        G4VMultipleScattering* p)
     402{
     403  if (1 < verbose) {
     404    G4cout << "G4LossTableManager::PreparePhysicsTable for "
     405           << particle->GetParticleName()
     406           << " and " << p->GetProcessName() << G4endl;
     407  }
     408  // start initialisation for the first run
     409  if( 0 == run ) {
     410    emConfigurator->PrepareModels(particle, p);
     411  }
     412}
     413
     414//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     415
     416void
     417G4LossTableManager::BuildPhysicsTable(const G4ParticleDefinition*)
     418{
     419  if(0 == run && startInitialisation) {
     420    emConfigurator->Clear();
     421  }
     422}
     423
     424//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     425
     426void G4LossTableManager::BuildPhysicsTable(
     427     const G4ParticleDefinition* aParticle,
     428     G4VEnergyLossProcess* p)
     429{
     430  if(1 < verbose) {
     431    G4cout << "### G4LossTableManager::BuildDEDXTable() is requested for "
     432           << aParticle->GetParticleName()
     433           << " and process " << p->GetProcessName()
     434           << G4endl;
     435  }
     436  // clear configurator
     437  if(0 == run && startInitialisation) {
     438    emConfigurator->Clear();
     439    firstParticle = aParticle;
     440  }
     441  startInitialisation = false;
     442
     443  // initialisation before any table is built
     444  if ( aParticle == firstParticle ) {
    356445    all_tables_are_built = true;
    357446
    358     if(1 < verbose)
    359       G4cout << "### G4LossTableManager start initilisation of tables"
     447    if(1 < verbose) {
     448      G4cout << "### G4LossTableManager start initilisation for first particle "
     449             << firstParticle->GetParticleName()
    360450             << G4endl;
    361     for (G4int i=0; i<n_loss; i++) {
     451    }
     452    for (G4int i=0; i<n_loss; ++i) {
    362453      G4VEnergyLossProcess* el = loss_vector[i];
    363454
    … …  
    365456        const G4ProcessManager* pm = el->GetProcessManager();
    366457        isActive[i] = pm->GetProcessActivation(el);
     458        if(0 == run) { base_part_vector[i] = el->BaseParticle(); }
    367459        tables_are_built[i] = false;
    368         all_tables_are_built = false;
    369         if(!isActive[i]) el->SetIonisation(false);
     460        all_tables_are_built= false;
     461        if(!isActive[i]) { el->SetIonisation(false); }
    370462 
    371463        if(1 < verbose) {
    … …  
    374466                 << "  active= " << pm->GetProcessActivation(el)
    375467                 << "  table= " << tables_are_built[i]
    376                  << "  isIonisation= " << el->IsIonisationProcess()
    377                  << G4endl;
     468                 << "  isIonisation= " << el->IsIonisationProcess();
     469          if(base_part_vector[i]) {
     470            G4cout << "  base particle " << base_part_vector[i]->GetParticleName();
     471          }
     472          G4cout << G4endl;
    378473        }
    379474      } else {
    … …  
    382477      }
    383478    }
    384     if (0 == run) firstParticle = particle;
    385     run++;
    386   }
    387 
    388   currentParticle = 0;
    389 
    390   SetParameters(p);
    391   for (G4int j=0; j<n_loss; j++) {
    392     if (p == loss_vector[j]) {
    393 
    394       if (!part_vector[j]) {
    395         part_vector[j] = particle;
    396         base_part_vector[j] = p->BaseParticle();
     479    ++run;
     480    currentParticle = 0;
     481  }
     482
     483  // Set run time parameters
     484  SetParameters(aParticle, p);
     485
     486  if (all_tables_are_built) { return; }
     487
     488  // Build tables for given particle
     489  all_tables_are_built = true;
     490
     491  for(G4int i=0; i<n_loss; ++i) {
     492    if(p == loss_vector[i] && !tables_are_built[i] && !base_part_vector[i]) {
     493      const G4ParticleDefinition* curr_part = part_vector[i];
     494      if(1 < verbose) {
     495        G4cout << "### BuildPhysicsTable for " << p->GetProcessName()
     496               << " and " << curr_part->GetParticleName()
     497               << " start BuildTable " << G4endl;
    397498      }
    398       if(maxEnergyForMuonsActive) {
    399         G4double dm = std::abs(particle->GetPDGMass() - 105.7*MeV);
    400         if(dm < 5.*MeV) p->SetMaxKinEnergy(maxKinEnergyForMuons);
    401       }
    402 
    403       if(1 < verbose) {
    404         G4cout << "For " << p->GetProcessName()
    405                << " for " << part_vector[j]->GetParticleName()
    406                << " tables_are_built= " << tables_are_built[j]
    407                << " procFlag= " << loss_vector[j]->TablesAreBuilt()
    408                << " all_tables_are_built= "     << all_tables_are_built
    409                << G4endl;
    410       }
    411     }
    412   }
    413 }
    414 
    415 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    416 
    417 G4EnergyLossMessenger* G4LossTableManager::GetMessenger()
    418 {
    419   return theMessenger;
    420 }
    421 
    422 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    423 
    424 void G4LossTableManager::ParticleHaveNoLoss(
    425      const G4ParticleDefinition* aParticle)
    426 {
    427   G4String s = " dE/dx table not found for "
    428              + aParticle->GetParticleName() + " !";
    429   G4Exception("G4LossTableManager::ParticleHaveNoLoss", "EM01",
    430               FatalException, s);
    431 
    432 }
    433 
    434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    435 
    436 G4bool G4LossTableManager::BuildCSDARange() const
    437 {
    438   return buildCSDARange;
    439 }
    440 
    441 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    442 
    443 void G4LossTableManager::BuildPhysicsTable(
    444      const G4ParticleDefinition* aParticle,
    445      G4VEnergyLossProcess* p)
    446 {
    447   if(1 < verbose) {
    448     G4cout << "### G4LossTableManager::BuildDEDXTable() is requested for "
    449            << aParticle->GetParticleName()
    450            << " and process " << p->GetProcessName()
    451            << G4endl;
    452   }
    453   if (all_tables_are_built) return;
    454   all_tables_are_built = true;
    455 
    456   for(G4int i=0; i<n_loss; i++) {
    457     if(!tables_are_built[i] && !base_part_vector[i]) {
    458       const G4ParticleDefinition* curr_part = part_vector[i];
    459499      G4VEnergyLossProcess* curr_proc = BuildTables(curr_part);
    460       if(curr_proc) CopyTables(curr_part, curr_proc);
    461     }
    462   }
    463 
    464   for (G4int ii=0; ii<n_loss; ii++) {
    465     if ( !tables_are_built[ii] ) {
    466       all_tables_are_built = false;
    467       break;
    468     }
     500      if(curr_proc) { CopyTables(curr_part, curr_proc); }
     501    }
     502    if ( !tables_are_built[i] ) { all_tables_are_built = false; }
    469503  }
    470504
    … …  
    473507           << "all_tables_are_built= " << all_tables_are_built
    474508           << G4endl;
    475 
    476     if(all_tables_are_built)
     509   
     510    if(all_tables_are_built) {
    477511      G4cout << "### All dEdx and Range tables are built #####" << G4endl;
     512    }
    478513  }
    479514}
    … …  
    484519                                          G4VEnergyLossProcess* base_proc)
    485520{
    486   for (G4int j=0; j<n_loss; j++) {
     521  for (G4int j=0; j<n_loss; ++j) {
    487522
    488523    G4VEnergyLossProcess* proc = loss_vector[j];
    489     if(proc == base_proc || proc->Particle() == part)
    490       tables_are_built[j] = true;
     524    //if(proc == base_proc || proc->Particle() == part)
     525    //  tables_are_built[j] = true;
    491526
    492527    if (!tables_are_built[j] && part == base_part_vector[j]) {
    … …  
    511546    }
    512547
    513     if (theElectron == part && theElectron == proc->SecondaryParticle() )
     548    if (theElectron == part && theElectron == proc->SecondaryParticle() ) {
    514549      proc->SetSecondaryRangeTable(base_proc->RangeTableForLoss());
     550    }
    515551  }
    516552}
    … …  
    535571  G4int i;
    536572
    537   for (i=0; i<n_loss; i++) {
     573  for (i=0; i<n_loss; ++i) {
    538574    p = loss_vector[i];
    539575    if (p && aParticle == part_vector[i] && !tables_are_built[i]) {
    … …  
    602638  std::vector<G4PhysicsTable*> listCSDA;
    603639
    604   for (i=0; i<n_dedx; i++) {
     640  for (i=0; i<n_dedx; ++i) {
    605641    p = loss_list[i];
    606642    p->SetIonisation(false);
    … …  
    658694//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    659695
     696G4EnergyLossMessenger* G4LossTableManager::GetMessenger()
     697{
     698  return theMessenger;
     699}
     700
     701//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     702
     703void G4LossTableManager::ParticleHaveNoLoss(
     704     const G4ParticleDefinition* aParticle)
     705{
     706  G4String s = " dE/dx table not found for "
     707             + aParticle->GetParticleName() + " !";
     708  G4Exception("G4LossTableManager::ParticleHaveNoLoss", "EM01",
     709              FatalException, s);
     710
     711}
     712
     713//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     714
     715G4bool G4LossTableManager::BuildCSDARange() const
     716{
     717  return buildCSDARange;
     718}
     719
     720//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     721
    660722void G4LossTableManager::SetLossFluctuations(G4bool val)
    661723{
    662724  lossFluctuationFlag = val;
    663   for(G4int i=0; i<n_loss; i++) {
    664     if(loss_vector[i]) loss_vector[i]->SetLossFluctuations(val);
    665   }
    666 }
    667 
    668 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    669 
    670 void G4LossTableManager::SetSubCutoff(G4bool val)
     725  for(G4int i=0; i<n_loss; ++i) {
     726    if(loss_vector[i]) { loss_vector[i]->SetLossFluctuations(val); }
     727  }
     728}
     729
     730//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     731
     732void G4LossTableManager::SetSubCutoff(G4bool val, const G4Region* r)
    671733{
    672734  subCutoffFlag = val;
    673   for(G4int i=0; i<n_loss; i++) {
    674     if(loss_vector[i]) loss_vector[i]->ActivateSubCutoff(val);
     735  for(G4int i=0; i<n_loss; ++i) {
     736    if(loss_vector[i]) { loss_vector[i]->ActivateSubCutoff(val, r); }
    675737  }
    676738}
    … …  
    682744  integral = val;
    683745  integralActive = true;
    684   for(G4int i=0; i<n_loss; i++) {
    685     if(loss_vector[i]) loss_vector[i]->SetIntegral(val);
     746  for(G4int i=0; i<n_loss; ++i) {
     747    if(loss_vector[i]) { loss_vector[i]->SetIntegral(val); }
    686748  }
    687749  size_t emp = emp_vector.size();
    688   for (size_t k=0; k<emp; k++) {
    689     if(emp_vector[k]) emp_vector[k]->SetIntegral(val);
     750  for (size_t k=0; k<emp; ++k) {
     751    if(emp_vector[k]) { emp_vector[k]->SetIntegral(val); }
    690752  }
    691753}
    … …  
    696758{
    697759  minSubRange = val;
    698   for(G4int i=0; i<n_loss; i++) {
    699     if(loss_vector[i]) loss_vector[i]->SetMinSubRange(val);
     760  for(G4int i=0; i<n_loss; ++i) {
     761    if(loss_vector[i]) { loss_vector[i]->SetMinSubRange(val); }
    700762  }
    701763}
    … …  
    706768{
    707769  rndmStepFlag = val;
    708   for(G4int i=0; i<n_loss; i++) {
    709     if(loss_vector[i]) loss_vector[i]->SetRandomStep(val);
     770  for(G4int i=0; i<n_loss; ++i) {
     771    if(loss_vector[i]) { loss_vector[i]->SetRandomStep(val); }
    710772  }
    711773}
    … …  
    717779  minEnergyActive = true;
    718780  minKinEnergy = val;
    719   for(G4int i=0; i<n_loss; i++) {
    720     if(loss_vector[i]) loss_vector[i]->SetMinKinEnergy(val);
     781  for(G4int i=0; i<n_loss; ++i) {
     782    if(loss_vector[i]) { loss_vector[i]->SetMinKinEnergy(val); }
    721783  }
    722784  size_t msc = msc_vector.size();
    723   for (size_t j=0; j<msc; j++) {
    724     if(msc_vector[j]) msc_vector[j]->SetMinKinEnergy(val);
     785  for (size_t j=0; j<msc; ++j) {
     786    if(msc_vector[j]) { msc_vector[j]->SetMinKinEnergy(val); }
    725787  }
    726788  size_t emp = emp_vector.size();
    727   for (size_t k=0; k<emp; k++) {
    728     if(emp_vector[k]) emp_vector[k]->SetMinKinEnergy(val);
     789  for (size_t k=0; k<emp; ++k) {
     790    if(emp_vector[k]) { emp_vector[k]->SetMinKinEnergy(val); }
    729791  }
    730792}
    … …  
    736798  maxEnergyActive = true;
    737799  maxKinEnergy = val;
    738   for(G4int i=0; i<n_loss; i++) {
    739     if(loss_vector[i]) loss_vector[i]->SetMaxKinEnergy(val);
     800  for(G4int i=0; i<n_loss; ++i) {
     801    if(loss_vector[i]) { loss_vector[i]->SetMaxKinEnergy(val); }
    740802  }
    741803  size_t msc = msc_vector.size();
    742   for (size_t j=0; j<msc; j++) {
    743     if(msc_vector[j]) msc_vector[j]->SetMaxKinEnergy(val);
     804  for (size_t j=0; j<msc; ++j) {
     805    if(msc_vector[j]) { msc_vector[j]->SetMaxKinEnergy(val); }
    744806  }
    745807  size_t emp = emp_vector.size();
    746   for (size_t k=0; k<emp; k++) {
    747     if(emp_vector[k]) emp_vector[k]->SetMaxKinEnergy(val);
     808  for (size_t k=0; k<emp; ++k) {
     809    if(emp_vector[k]) { emp_vector[k]->SetMaxKinEnergy(val); }
    748810  }
    749811}
    … …  
    753815void G4LossTableManager::SetMaxEnergyForCSDARange(G4double val)
    754816{
    755   for(G4int i=0; i<n_loss; i++) {
    756     if(loss_vector[i]) loss_vector[i]->SetMaxKinEnergyForCSDARange(val);
     817  for(G4int i=0; i<n_loss; ++i) {
     818    if(loss_vector[i]) { loss_vector[i]->SetMaxKinEnergyForCSDARange(val); }
    757819  }
    758820}
    … …  
    770832void G4LossTableManager::SetDEDXBinning(G4int val)
    771833{
    772   for(G4int i=0; i<n_loss; i++) {
    773     if(loss_vector[i]) loss_vector[i]->SetDEDXBinning(val);
     834  for(G4int i=0; i<n_loss; ++i) {
     835    if(loss_vector[i]) { loss_vector[i]->SetDEDXBinning(val); }
    774836  }
    775837}
    … …  
    779841void G4LossTableManager::SetDEDXBinningForCSDARange(G4int val)
    780842{
    781   for(G4int i=0; i<n_loss; i++) {
    782     if(loss_vector[i]) loss_vector[i]->SetDEDXBinningForCSDARange(val);
     843  for(G4int i=0; i<n_loss; ++i) {
     844    if(loss_vector[i]) { loss_vector[i]->SetDEDXBinningForCSDARange(val); }
    783845  }
    784846}
    … …  
    788850void G4LossTableManager::SetLambdaBinning(G4int val)
    789851{
     852  G4int n = val/G4int(std::log10(maxKinEnergy/minKinEnergy));
     853  if(n < 5) {
     854    G4cout << "G4LossTableManager::SetLambdaBinning WARNING "
     855           << "too small number of bins " << val << "  ignored"
     856           << G4endl;
     857    return;
     858  }
     859  nbinsLambda = val;
     860  nbinsPerDecade = n;
    790861  size_t msc = msc_vector.size();
    791   for (size_t j=0; j<msc; j++) {
    792     if(msc_vector[j]) msc_vector[j]->SetBinning(val);
     862  for (size_t j=0; j<msc; ++j) {
     863    if(msc_vector[j]) { msc_vector[j]->SetBinning(val); }
    793864  }
    794865  size_t emp = emp_vector.size();
    795   for (size_t k=0; k<emp; k++) {
    796     if(emp_vector[k]) emp_vector[k]->SetLambdaBinning(val);
    797   }
     866  for (size_t k=0; k<emp; ++k) {
     867    if(emp_vector[k]) { emp_vector[k]->SetLambdaBinning(val); }
     868  }
     869}
     870
     871//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     872
     873G4int G4LossTableManager::GetNumberOfBinsPerDecade() const
     874{
     875  return nbinsPerDecade;
    798876}
    799877
    … …  
    803881{
    804882  verbose = val;
    805   for(G4int i=0; i<n_loss; i++) {
    806     if(loss_vector[i]) loss_vector[i]->SetVerboseLevel(val);
     883  for(G4int i=0; i<n_loss; ++i) {
     884    if(loss_vector[i]) { loss_vector[i]->SetVerboseLevel(val); }
    807885  }
    808886  size_t msc = msc_vector.size();
    809   for (size_t j=0; j<msc; j++) {
    810     if(msc_vector[j]) msc_vector[j]->SetVerboseLevel(val);
     887  for (size_t j=0; j<msc; ++j) {
     888    if(msc_vector[j]) { msc_vector[j]->SetVerboseLevel(val); }
    811889  }
    812890  size_t emp = emp_vector.size();
    813   for (size_t k=0; k<emp; k++) {
    814     if(emp_vector[k]) emp_vector[k]->SetVerboseLevel(val);
    815   }
     891  for (size_t k=0; k<emp; ++k) {
     892    if(emp_vector[k]) { emp_vector[k]->SetVerboseLevel(val); }
     893  }
     894  emConfigurator->SetVerbose(val);
     895  //tableBuilder->SetVerbose(val);
     896  //emCorrections->SetVerbose(val);
     897  emSaturation->SetVerbose(val);
    816898}
    817899
    … …  
    823905  maxRangeVariation = v1;
    824906  maxFinalStep = v2;
    825   for(G4int i=0; i<n_loss; i++) {
    826     if(loss_vector[i]) loss_vector[i]->SetStepFunction(v1, v2);
     907  for(G4int i=0; i<n_loss; ++i) {
     908    if(loss_vector[i]) { loss_vector[i]->SetStepFunction(v1, v2); }
    827909  }
    828910}
    … …  
    832914void G4LossTableManager::SetLinearLossLimit(G4double val)
    833915{
    834   for(G4int i=0; i<n_loss; i++) {
    835     if(loss_vector[i]) loss_vector[i]->SetLinearLossLimit(val);
     916  for(G4int i=0; i<n_loss; ++i) {
     917    if(loss_vector[i]) { loss_vector[i]->SetLinearLossLimit(val); }
    836918  }
    837919}
    … …  
    846928//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    847929
    848 void G4LossTableManager::SetParameters(G4VEnergyLossProcess* p)
    849 {
    850   if(stepFunctionActive) p->SetStepFunction(maxRangeVariation, maxFinalStep);
    851   if(integralActive)     p->SetIntegral(integral);
    852   if(minEnergyActive)    p->SetMinKinEnergy(minKinEnergy);
    853   if(maxEnergyActive)    p->SetMaxKinEnergy(maxKinEnergy);
     930void
     931G4LossTableManager::SetParameters(const G4ParticleDefinition* aParticle,
     932                                  G4VEnergyLossProcess* p)
     933{
     934  if(stepFunctionActive) { p->SetStepFunction(maxRangeVariation, maxFinalStep); }
     935  if(integralActive)     { p->SetIntegral(integral); }
     936  if(minEnergyActive)    { p->SetMinKinEnergy(minKinEnergy); }
     937  if(maxEnergyActive)    { p->SetMaxKinEnergy(maxKinEnergy); }
    854938  p->SetVerboseLevel(verbose);
     939  if(maxEnergyForMuonsActive) {
     940    G4double dm = std::abs(aParticle->GetPDGMass() - 105.7*MeV);
     941    if(dm < 5.*MeV) { p->SetMaxKinEnergy(maxKinEnergyForMuons); }
     942  }
    855943}
    856944
    … …  
    858946
    859947const std::vector<G4VEnergyLossProcess*>&
    860       G4LossTableManager::GetEnergyLossProcessVector()
     948G4LossTableManager::GetEnergyLossProcessVector()
    861949{
    862950  return loss_vector;
    … …  
    873961
    874962const std::vector<G4VMultipleScattering*>&
    875       G4LossTableManager::GetMultipleScatteringVector()
     963G4LossTableManager::GetMultipleScatteringVector()
    876964{
    877965  return msc_vector;
    … …  
    9251013void G4LossTableManager::SetFactorForAngleLimit(G4double val)
    9261014{
    927   if(val > 0.0) factorForAngleLimit = val;
     1015  if(val > 0.0) { factorForAngleLimit = val; }
    9281016}
    9291017

  • trunk/source/processes/electromagnetic/utils/src/G4VEmModel.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmModel.cc,v 1.30 2009/09/23 14:42:47 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4VEmModel.cc,v 1.35 2010/05/26 18:06:25 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8080  G4int n = elmSelectors.size();
    8181  if(n > 0) {
    82     for(G4int i=0; i<n; i++) {
     82    for(G4int i=0; i<n; ++i) {
    8383      delete elmSelectors[i];
    8484    }
    … …  
    121121  // initialise before run
    122122  flagDeexcitation = false;
    123 
    124   G4int nbins = G4int(std::log10(highLimit/lowLimit) + 0.5);
    125   if(nbins < 3) nbins = 3;
    126   G4bool spline = G4LossTableManager::Instance()->SplineFlag();
     123  G4LossTableManager* man = G4LossTableManager::Instance();
     124  G4bool spline = man->SplineFlag();
     125
     126  // two times less bins because probability functon is normalized
     127  // so correspondingly is more smooth
     128  G4int nbins = (man->GetNumberOfBinsPerDecade()/3)*
     129    G4int(std::log10(highLimit/lowLimit) + 0.5);
     130  if(nbins < 5) { nbins = 5; }
    127131
    128132  G4ProductionCutsTable* theCoupleTable=
    … …  
    131135
    132136  // prepare vector
    133   if(numOfCouples > nSelectors) elmSelectors.reserve(numOfCouples);
     137  if(numOfCouples > nSelectors) {
     138    elmSelectors.reserve(numOfCouples);
     139    for(G4int i=nSelectors; i<numOfCouples; ++i) { elmSelectors.push_back(0); }
     140    nSelectors = numOfCouples;
     141  }
    134142
    135143  // initialise vector
    136   for(G4int i=0; i<numOfCouples; i++) {
     144  for(G4int i=0; i<numOfCouples; ++i) {
    137145    currentCouple = theCoupleTable->GetMaterialCutsCouple(i);
    138146    const G4Material* material = currentCouple->GetMaterial();
    … …  
    141149    // selector already exist check if should be deleted
    142150    G4bool create = true;
    143     if(i < nSelectors) {
    144       if(elmSelectors[i]) {
    145         if(material == elmSelectors[i]->GetMaterial()) create = false;
    146         else delete elmSelectors[i];
    147       }
    148     } else {
    149       nSelectors++;
    150       elmSelectors.push_back(0);
     151    if(elmSelectors[i]) {
     152      if(material == elmSelectors[i]->GetMaterial()) { create = false; }
     153      else { delete elmSelectors[i]; }
    151154    }
    152155    if(create) {
    … …  
    195198//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    196199
     200const G4Element* G4VEmModel::SelectRandomAtom(const G4Material* material,
     201                                              const G4ParticleDefinition* pd,
     202                                              G4double kinEnergy,
     203                                              G4double tcut,
     204                                              G4double tmax)
     205{
     206  const G4ElementVector* theElementVector = material->GetElementVector();
     207  G4int n = material->GetNumberOfElements() - 1;
     208  currentElement = (*theElementVector)[n];
     209  if (n > 0) {
     210    G4double x = G4UniformRand()*
     211                 G4VEmModel::CrossSectionPerVolume(material,pd,kinEnergy,tcut,tmax);
     212    for(G4int i=0; i<n; ++i) {
     213      if (x <= xsec[i]) {
     214        currentElement = (*theElementVector)[i];
     215        break;
     216      }
     217    }
     218  }
     219  return currentElement;
     220}
     221
     222//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     223
    197224G4double G4VEmModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
    198225                                                G4double, G4double, G4double,
    … …  
    217244//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    218245
     246G4double G4VEmModel::ChargeSquareRatio(const G4Track& track)
     247{
     248  return GetChargeSquareRatio(track.GetDefinition(),
     249                              track.GetMaterial(), track.GetKineticEnergy());
     250}
     251
     252//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     253
    219254G4double G4VEmModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
    220255                                          const G4Material*, G4double)

  • trunk/source/processes/electromagnetic/utils/src/G4VEmProcess.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmProcess.cc,v 1.79 2009/11/10 20:30:55 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4VEmProcess.cc,v 1.87 2010/05/10 13:35:32 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5353// 27-10-07 Virtual functions moved to source (V.Ivanchenko)
    5454// 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
     55// 17-02-10 Added pointer currentParticle (VI)
    5556//
    5657// Class Description:
    … …  
    9899  currentModel(0),
    99100  particle(0),
     101  currentParticle(0),
    100102  currentCouple(0)
    101103{
    … …  
    165167  G4VEmFluctuationModel* fm = 0;
    166168  modelManager->AddEmModel(order, p, fm, region);
    167   if(p) p->SetParticleChange(pParticleChange);
     169  if(p) { p->SetParticleChange(pParticleChange); }
    168170}
    169171
    … …  
    205207void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
    206208{
    207   if(!particle) particle = &part;
     209  if(!particle) { SetParticle(&part); }
    208210  if(1 < verboseLevel) {
    209211    G4cout << "G4VEmProcess::PreparePhysicsTable() for "
    … …  
    214216  }
    215217
    216   (G4LossTableManager::Instance())->EmConfigurator()->AddModels();
     218  (G4LossTableManager::Instance())->PreparePhysicsTable(&part, this);
    217219
    218220  if(particle == &part) {
    … …  
    230232    }
    231233
    232     theCuts = modelManager->Initialise(particle,secondaryParticle,2.,verboseLevel);
     234    theCuts = modelManager->Initialise(particle,secondaryParticle,
     235                                       2.,verboseLevel);
    233236    const G4ProductionCutsTable* theCoupleTable=
    234237          G4ProductionCutsTable::GetProductionCutsTable();
    … …  
    242245    }
    243246  }
    244   // Sub Cutoff and Deexcitation
     247  // Deexcitation
    245248  if (nDERegions>0) {
    246249
    … …  
    287290  }
    288291
     292  (G4LossTableManager::Instance())->BuildPhysicsTable(particle);
    289293  if(buildLambdaTable) {
    290294    BuildLambdaTable();
    … …  
    294298  // reduce printout for nuclear stopping
    295299  G4bool gproc = true;
    296   if(GetProcessName() == "nuclearStopping" &&
     300  G4int st = GetProcessSubType();
     301  if(st == fCoulombScattering && part.GetParticleType() == "nucleus" &&
    297302     partname != "GenericIon" && partname != "alpha") { gproc = false; }
    298303
    … …  
    354359           << particle->GetParticleName()
    355360           << G4endl;
    356     if(2 < verboseLevel) {
    357       G4cout << *theLambdaTable << G4endl;
    358     }
    359361  }
    360362}
    … …  
    385387  if(verboseLevel > 2 && buildLambdaTable) {
    386388    G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
    387     if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
     389    if(theLambdaTable) { G4cout << (*theLambdaTable) << G4endl; }
    388390  }
    389391}
    … …  
    399401  *condition = NotForced;
    400402  G4double x = DBL_MAX;
    401   if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
     403  if(previousStepSize <= DBL_MIN) { theNumberOfInteractionLengthLeft = -1.0; }
    402404  InitialiseStep(track);
    403   if(!currentModel->IsActive(preStepKinEnergy)) return x;
     405  if(!currentModel->IsActive(preStepKinEnergy)) { return x; }
    404406
    405407  if(preStepKinEnergy < mfpKinEnergy) {
    … …  
    428430      G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
    429431      G4cout << "[ " << GetProcessName() << "]" << G4endl;
    430       G4cout << " for " << particle->GetParticleName()
     432      G4cout << " for " << currentParticle->GetParticleName()
    431433             << " in Material  " <<  currentMaterial->GetName()
    432434             << " Ekin(MeV)= " << preStepKinEnergy/MeV
    … …  
    461463  // Do not make anything if particle is stopped, the annihilation then
    462464  // should be performed by the AtRestDoIt!
    463   if (track.GetTrackStatus() == fStopButAlive) return &fParticleChange;
     465  if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; }
    464466
    465467  G4double finalT = track.GetKineticEnergy();
    … …  
    469471    G4double lx = GetLambda(finalT, currentCouple);
    470472    if(preStepLambda<lx && 1 < verboseLevel) {
    471       G4cout << "WARING: for " << particle->GetParticleName()
     473      G4cout << "WARING: for " << currentParticle->GetParticleName()
    472474             << " and " << GetProcessName()
    473475             << " E(MeV)= " << finalT/MeV
    … …  
    483485
    484486  SelectModel(finalT, currentCoupleIndex);
    485   if(!currentModel->IsActive(finalT)) return &fParticleChange;
     487  if(!currentModel->IsActive(finalT)) { return &fParticleChange; }
    486488  if(useDeexcitation) {
    487489    currentModel->SetDeexcitationFlag(idxDERegions[currentCoupleIndex]);
    … …  
    662664    SelectModel(kineticEnergy, currentCoupleIndex);
    663665    cross = currentModel->CrossSectionPerVolume(currentMaterial,
    664                                                 particle,kineticEnergy);
    665   }
    666 
     666                                                currentParticle,kineticEnergy);
     667  }
     668
     669  if(cross < 0.0) { cross = 0.0; }
    667670  return cross;
    668671}
    … …  
    676679  *condition = NotForced;
    677680  return G4VEmProcess::MeanFreePath(track);
     681}
     682
     683//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     684
     685G4double G4VEmProcess::MeanFreePath(const G4Track& track)
     686{
     687  DefineMaterial(track.GetMaterialCutsCouple());
     688  preStepLambda = GetCurrentLambda(track.GetKineticEnergy());
     689  G4double x = DBL_MAX;
     690  if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
     691  return x;
     692}
     693
     694//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     695
     696G4double
     697G4VEmProcess::ComputeCrossSectionPerAtom(G4double kineticEnergy,
     698                                         G4double Z, G4double A, G4double cut)
     699{
     700  SelectModel(kineticEnergy, currentCoupleIndex);
     701  G4double x = 0.0;
     702  if(currentModel) {
     703   x = currentModel->ComputeCrossSectionPerAtom(currentParticle,kineticEnergy,
     704                                                 Z,A,cut);
     705  }
     706  return x;
    678707}
    679708
    … …  
    714743    theEnergyOfCrossSectionMax[i] = emax;
    715744    theCrossSectionMax[i] = smax;
    716     if(2 < verboseLevel) {
     745    if(1 < verboseLevel) {
    717746      G4cout << "For " << particle->GetParticleName()
    718747             << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
    … …  
    733762
    734763//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     764
     765const G4Element* G4VEmProcess::GetCurrentElement() const
     766{
     767  const G4Element* elm = 0;
     768  if(currentModel) {elm = currentModel->GetCurrentElement(); }
     769  return elm;
     770}
     771
     772//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

  • trunk/source/processes/electromagnetic/utils/src/G4VEnergyLossProcess.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEnergyLossProcess.cc,v 1.158 2009/10/29 18:07:08 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4VEnergyLossProcess.cc,v 1.167 2010/05/26 10:41:34 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    336336{
    337337  modelManager->AddEmModel(order, p, fluc, region);
    338   if(p) p->SetParticleChange(pParticleChange, fluc);
     338  if(p) { p->SetParticleChange(pParticleChange, fluc); }
    339339}
    340340
    … …  
    406406    particle = &part;
    407407    if(part.GetParticleType() == "nucleus") {
    408       if(!theGenericIon) theGenericIon = G4GenericIon::GenericIon();
     408      if(!theGenericIon) { theGenericIon = G4GenericIon::GenericIon(); }
    409409      if(particle == theGenericIon) { isIon = true; }
    410410      else if(part.GetPDGCharge() > eplus) {
    … …  
    426426      lManager->RegisterExtraParticle(&part, this);
    427427    }
     428    if(1 < verboseLevel) {
     429      G4cout << "### G4VEnergyLossProcess::PreparePhysicsTable() end for "
     430             << part.GetParticleName() << G4endl;
     431    }
    428432    return;
    429433  }
    430434
    431435  Clean();
    432   lManager->EmConfigurator()->AddModels();
     436  lManager->PreparePhysicsTable(&part, this);
    433437
    434438  // Base particle and set of models can be defined here
    … …  
    502506        G4bool reg = false;
    503507        for(G4int i=0; i<nSCoffRegions; ++i) {
    504           if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = true;
     508          if( pcuts == scoffRegions[i]->GetProductionCuts()) { reg = true; }
    505509        }
    506510        idxSCoffRegions[j] = reg;
    … …  
    509513        G4bool reg = false;
    510514        for(G4int i=0; i<nDERegions; ++i) {
    511           if( pcuts == deRegions[i]->GetProductionCuts()) reg = true;
     515          if( pcuts == deRegions[i]->GetProductionCuts()) { reg = true; }
    512516        }
    513517        idxDERegions[j] = reg;
    … …  
    515519    }
    516520  }
    517 
    518   lManager->EnergyLossProcessIsInitialised(particle, this);
    519521
    520522  if (1 < verboseLevel) {
    … …  
    568570
    569571  // Added tracking cut to avoid tracking artifacts
    570   if(isIonisation) fParticleChange.SetLowEnergyLimit(lowestKinEnergy);
     572  if(isIonisation) { fParticleChange.SetLowEnergyLimit(lowestKinEnergy); }
    571573
    572574  if(1 < verboseLevel) {
    … …  
    574576           << GetProcessName()
    575577           << " and particle " << part.GetParticleName();
    576     if(isIonisation) G4cout << "  isIonisation  flag = 1";
     578    if(isIonisation) { G4cout << "  isIonisation  flag = 1"; }
    577579    G4cout << G4endl;
    578580  }
    … …  
    590592  }
    591593  G4PhysicsTable* table = 0;
    592   G4double emin = minKinEnergy;
    593594  G4double emax = maxKinEnergy;
    594595  G4int bin = nBins;
    … …  
    619620    G4cout << numOfCouples << " materials"
    620621           << " minKinEnergy= " << minKinEnergy
    621            << " maxKinEnergy= " << maxKinEnergy
     622           << " maxKinEnergy= " << emax
     623           << " nbin= " << bin
    622624           << " EmTableType= " << tType
    623625           << " table= " << table
    … …  
    642644        theCoupleTable->GetMaterialCutsCouple(i);
    643645      if(!bVector) {
    644         aVector = new G4PhysicsLogVector(emin, emax, bin);
     646        aVector = new G4PhysicsLogVector(minKinEnergy, emax, bin);
    645647        bVector = aVector;
    646648      } else {
    647649        aVector = new G4PhysicsLogVector(*bVector);
    648650      }
    649       //      G4PhysicsVector* aVector = new G4PhysicsLogVector(emin, emax, bin);
    650651      aVector->SetSpline(splineFlag);
    651652
    652653      modelManager->FillDEDXVector(aVector, couple, tType);
    653       if(splineFlag) aVector->FillSecondDerivatives();
     654      if(splineFlag) { aVector->FillSecondDerivatives(); }
    654655
    655656      // Insert vector for this material into the table
    … …  
    701702
    702703  G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
     704  G4PhysicsLogVector* aVector = 0;
     705  G4PhysicsLogVector* bVector = 0;
    703706
    704707  for(size_t i=0; i<numOfCouples; ++i) {
    … …  
    709712      const G4MaterialCutsCouple* couple =
    710713        theCoupleTable->GetMaterialCutsCouple(i);
    711       G4double cut = (*theCuts)[i];
    712       if(fSubRestricted == tType) cut = (*theSubCuts)[i];
    713       G4PhysicsVector* aVector = LambdaPhysicsVector(couple, cut);
     714      if(!bVector) {
     715        aVector = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nBins);
     716        bVector = aVector;
     717      } else {
     718        aVector = new G4PhysicsLogVector(*bVector);
     719      }
    714720      aVector->SetSpline(splineFlag);
    715721
    716722      modelManager->FillLambdaVector(aVector, couple, true, tType);
    717       if(splineFlag) aVector->FillSecondDerivatives();
     723      if(splineFlag) { aVector->FillSecondDerivatives(); }
    718724
    719725      // Insert vector for this material into the table
    … …  
    881887    if(x > finalRange && y < currentMinStep) {
    882888      x = y + finalRange*(1.0 - dRoverRange)*(2.0 - finalRange/fRange);
    883     } else if (rndmStepFlag) {x = SampleRange();}
     889    } else if (rndmStepFlag) { x = SampleRange(); }
    884890    //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
    885891    //  <<" range= "<<fRange <<" cMinSt="<<currentMinStep
    … …  
    912918  preStepScaledEnergy = preStepKinEnergy*massRatio;
    913919  SelectModel(preStepScaledEnergy);
    914   if(!currentModel->IsActive(preStepScaledEnergy)) return x;
    915 
    916   if(isIon) {
    917     chargeSqRatio =
    918       currentModel->GetChargeSquareRatio(currPart,currentMaterial,preStepKinEnergy);
     920  if(!currentModel->IsActive(preStepScaledEnergy)) { return x; }
     921
     922  if(isIon) {
     923    chargeSqRatio = currentModel->ChargeSquareRatio(track);
    919924    reduceFactor  = 1.0/(chargeSqRatio*massRatio);
    920925  }
    921926  //G4cout << "q2= " << chargeSqRatio << " massRatio= " << massRatio << G4endl;
    922927  // initialisation for sampling of the interaction length
    923   if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
    924   if(theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX;
     928  if(previousStepSize <= DBL_MIN) { theNumberOfInteractionLengthLeft = -1.0; }
     929  if(theNumberOfInteractionLengthLeft < 0.0) { mfpKinEnergy = DBL_MAX; }
    925930
    926931  // compute mean free path
    927932  if(preStepScaledEnergy < mfpKinEnergy) {
    928     if (integral) ComputeLambdaForScaledEnergy(preStepScaledEnergy);
    929     else  preStepLambda = GetLambdaForScaledEnergy(preStepScaledEnergy);
    930     if(preStepLambda <= DBL_MIN) mfpKinEnergy = 0.0;
     933    if (integral) { ComputeLambdaForScaledEnergy(preStepScaledEnergy); }
     934    else  { preStepLambda = GetLambdaForScaledEnergy(preStepScaledEnergy); }
     935    if(preStepLambda <= DBL_MIN) { mfpKinEnergy = 0.0; }
    931936  }
    932937
    … …  
    940945      // subtract NumberOfInteractionLengthLeft
    941946      SubtractNumberOfInteractionLengthLeft(previousStepSize);
    942       if(theNumberOfInteractionLengthLeft < 0.)
     947      if(theNumberOfInteractionLengthLeft < 0.) {
    943948        theNumberOfInteractionLengthLeft = perMillion;
     949      }
    944950    }
    945951
    … …  
    966972      // subtract NumberOfInteractionLengthLeft
    967973      SubtractNumberOfInteractionLengthLeft(previousStepSize);
    968       if(theNumberOfInteractionLengthLeft < 0.)
     974      if(theNumberOfInteractionLengthLeft < 0.) {
    969975        theNumberOfInteractionLengthLeft = perMillion;
     976      }
    970977    }
    971978    currentInteractionLength = DBL_MAX;
    … …  
    987994  // Get the actual (true) Step length
    988995  G4double length = step.GetStepLength();
    989   if(length <= DBL_MIN) return &fParticleChange;
     996  if(length <= DBL_MIN) { return &fParticleChange; }
    990997  G4double eloss  = 0.0;
    991   G4double esecdep = 0.0;
    992998 
    993999  /* 
    … …  
    10711077
    10721078      // Check boundary
    1073       if(prePoint->GetStepStatus() == fGeomBoundary) yes = true;
     1079      if(prePoint->GetStepStatus() == fGeomBoundary) { yes = true; }
    10741080
    10751081      // Check PrePoint
    … …  
    10831089        }
    10841090
    1085         if(preSafety < rcut) yes = true;
     1091        if(preSafety < rcut) { yes = true; }
    10861092
    10871093        // Check PostPoint
    … …  
    10911097            postSafety =
    10921098              safetyHelper->ComputeSafety(step.GetPostStepPoint()->GetPosition());
    1093             if(postSafety < rcut) yes = true;
     1099            if(postSafety < rcut) { yes = true; }
    10941100          }
    10951101        }
    … …  
    11031109        scTracks.clear();
    11041110        SampleSubCutSecondaries(scTracks, step,
    1105                                 currentModel,currentMaterialIndex,
    1106                                 esecdep);
     1111                                currentModel,currentMaterialIndex);
    11071112        // add bremsstrahlung sampling
    11081113        /*
    … …  
    11121117                scTracks, step, (scProcesses[i])->
    11131118                SelectModelForMaterial(preStepKinEnergy, currentMaterialIndex),
    1114                 currentMaterialIndex,esecdep);
     1119                currentMaterialIndex);
    11151120          }
    11161121        }
    … …  
    11231128            G4Track* t = scTracks[i];
    11241129            G4double e = t->GetKineticEnergy();
    1125             if (t->GetDefinition() == thePositron) e += 2.0*electron_mass_c2;
     1130            if (t->GetDefinition() == thePositron) { e += 2.0*electron_mass_c2; }
    11261131            esec += e;
    11271132            pParticleChange->AddSecondary(t);
    … …  
    11331138
    11341139  // Corrections, which cannot be tabulated
    1135   currentModel->CorrectionsAlongStep(currentCouple, dynParticle,
    1136                                      eloss, esecdep, length);
     1140  if(isIon) {
     1141    G4double eadd = 0.0;
     1142    currentModel->CorrectionsAlongStep(currentCouple, dynParticle,
     1143                                       eloss, eadd, length);
     1144  }
    11371145
    11381146  // Sample fluctuations
    … …  
    11401148    G4VEmFluctuationModel* fluc = currentModel->GetModelOfFluctuations();
    11411149    if(fluc &&
    1142       (eloss + esec + esecdep + lowestKinEnergy) < preStepKinEnergy) {
     1150      (eloss + esec + lowestKinEnergy) < preStepKinEnergy) {
    11431151
    11441152      G4double tmax =
    … …  
    11581166    }
    11591167  }
    1160   // add low-energy subcutoff particles
    1161   eloss += esecdep;
    1162   if(eloss < 0.0) eloss = 0.0;
    11631168
    11641169  // deexcitation
    1165   else if (useDeexcitation) {
    1166     if(idxDERegions[currentMaterialIndex]) {
     1170  if (useDeexcitation) {
     1171    if(idxDERegions[currentMaterialIndex] && eloss > 0.0) {
    11671172      currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
    1168       if(eloss < 0.0) eloss = 0.0;
     1173      if(eloss < 0.0) { eloss = 0.0; }
    11691174    }
    11701175  }
    … …  
    12031208       const G4Step& step,
    12041209       G4VEmModel* model,
    1205        G4int idx,
    1206        G4double& /*extraEdep*/)
     1210       G4int idx)
    12071211{
    12081212  // Fast check weather subcutoff can work
    12091213  G4double subcut = (*theSubCuts)[idx];
    12101214  G4double cut = (*theCuts)[idx];
    1211   if(cut <= subcut) return;
     1215  if(cut <= subcut) { return; }
    12121216
    12131217  const G4Track* track = step.GetTrack();
    … …  
    12181222
    12191223  // negligible probability to get any interaction
    1220   if(length*cross < perMillion) return;
     1224  if(length*cross < perMillion) { return; }
    12211225  /*     
    12221226  if(-1 < verboseLevel)
    … …  
    12691273      if(addSec) {
    12701274        G4Track* t = new G4Track((*it), pretime + fragment*dt, r);
    1271         //G4Track* t = new G4Track((*it), pretime, r);
    12721275        t->SetTouchableHandle(track->GetTouchableHandle());
    12731276        tracks.push_back(t);
    … …  
    12961299  G4double postStepScaledEnergy = finalT*massRatio;
    12971300
    1298   if(!currentModel->IsActive(postStepScaledEnergy)) return &fParticleChange;
     1301  if(!currentModel->IsActive(postStepScaledEnergy)) { return &fParticleChange; }
    12991302  /*
    13001303  if(-1 < verboseLevel) {
    … …  
    15061509                                    G4bool mandatory)
    15071510{
    1508   G4bool res = true;
     1511  G4bool isRetrieved = false;
    15091512  G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
    1510   G4bool yes = aTable->ExistPhysicsTable(filename);
    1511   if(yes) {
    1512     if(!aTable) aTable = G4PhysicsTableHelper::PreparePhysicsTable(0);
    1513     yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii);
    1514 
    1515     if((G4LossTableManager::Instance())->SplineFlag()) {
    1516       size_t n = aTable->length();
    1517       for(size_t i=0; i<n; ++i) {
    1518         if((*aTable)[i]) {
    1519           (*aTable)[i]->SetSpline(true);
     1513  if(aTable) {
     1514    if(aTable->ExistPhysicsTable(filename)) {
     1515      if(G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii)) {
     1516        isRetrieved = true;
     1517        if((G4LossTableManager::Instance())->SplineFlag()) {
     1518          size_t n = aTable->length();
     1519          for(size_t i=0; i<n; ++i) {
     1520            if((*aTable)[i]) { (*aTable)[i]->SetSpline(true); }
     1521          }
    15201522        }
    1521       }
    1522     }
    1523   }
    1524   if(yes) {
    1525     if (0 < verboseLevel) {
    1526       G4cout << tname << " table for " << part->GetParticleName()
    1527              << " is Retrieved from <" << filename << ">"
    1528              << G4endl;
    1529     }
    1530   } else {
    1531     if(mandatory) res = false;
    1532     if(mandatory || 1 < verboseLevel) {
     1523        if (0 < verboseLevel) {
     1524          G4cout << tname << " table for " << part->GetParticleName()
     1525                 << " is Retrieved from <" << filename << ">"
     1526                 << G4endl;
     1527        }
     1528      }
     1529    }
     1530  }
     1531  if(mandatory && !isRetrieved) {
     1532    if(0 < verboseLevel) {
    15331533      G4cout << tname << " table for " << part->GetParticleName()
    15341534             << " from file <"
    … …  
    15361536             << G4endl;
    15371537    }
    1538   }
    1539   return res;
     1538    return false;
     1539  }
     1540  return true;
    15401541}
    15411542
    … …  
    15741575                                                (*theCuts)[currentMaterialIndex]);
    15751576  }
    1576 
     1577  if(cross < 0.0) { cross = 0.0; }
    15771578  return cross;
    15781579}
    … …  
    15851586  preStepLambda = GetLambdaForScaledEnergy(track.GetKineticEnergy()*massRatio);
    15861587  G4double x = DBL_MAX;
    1587   if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
     1588  if(DBL_MIN < preStepLambda) { x = 1.0/preStepLambda; }
    15881589  return x;
    15891590}
    … …  
    16221623//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    16231624
    1624 G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
    1625                  const G4MaterialCutsCouple* couple, G4double cut)
    1626 {
     1625G4PhysicsVector*
     1626G4VEnergyLossProcess::LambdaPhysicsVector(const G4MaterialCutsCouple* /*couple*/,
     1627                                          G4double /*cut*/)
     1628{
     1629  /*
    16271630  G4double tmin =
    16281631    std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
    16291632             minKinEnergy);
    1630   if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
     1633  if(tmin >= maxKinEnergy) { tmin = 0.5*maxKinEnergy; }
    16311634  G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
     1635  */
     1636  G4PhysicsVector* v = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nBins);
    16321637  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
    16331638  return v;
    … …  
    16401645{
    16411646  G4bool add = true;
    1642   if(p->GetProcessName() != "eBrem") add = false;
     1647  if(p->GetProcessName() != "eBrem") { add = false; }
    16431648  if(add && nProcesses > 0) {
    16441649    for(G4int i=0; i<nProcesses; ++i) {
    … …  
    16991704void G4VEnergyLossProcess::SetCSDARangeTable(G4PhysicsTable* p)
    17001705{
    1701   if(theCSDARangeTable != p) theCSDARangeTable = p;
     1706  if(theCSDARangeTable != p) { theCSDARangeTable = p; }
    17021707
    17031708  if(p) {
    … …  
    17681773           << " and process " << GetProcessName() << G4endl;
    17691774  }
    1770   if(theLambdaTable != p) theLambdaTable = p;
     1775  if(theLambdaTable != p) { theLambdaTable = p; }
    17711776  tablesAreBuilt = true;
    17721777
    … …  
    18221827//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    18231828
     1829const G4Element* G4VEnergyLossProcess::GetCurrentElement() const
     1830{
     1831  const G4Element* elm = 0;
     1832  if(currentModel) { elm = currentModel->GetCurrentElement(); }
     1833  return elm;
     1834}
     1835
     1836//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1837

  • trunk/source/processes/electromagnetic/utils/src/G4VMscModel.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMscModel.cc,v 1.13 2009/07/20 17:32:47 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4VMscModel.cc,v 1.16 2010/04/06 09:24:46 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6161  facgeom(2.5),
    6262  facsafety(0.3),
    63   skin(3.0),
     63  skin(1.0),
    6464  dtrl(0.05),
    6565  lambdalimit(mm),

  • trunk/source/processes/electromagnetic/utils/src/G4VMultipleScattering.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMultipleScattering.cc,v 1.77 2009/10/29 18:07:08 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4VMultipleScattering.cc,v 1.82 2010/04/12 11:45:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9595  firstParticle(0),
    9696  stepLimit(fUseSafety),
    97   skin(3.0),
     97  skin(1.0),
    9898  facrange(0.04),
    9999  facgeom(2.5),
    … …  
    144144  G4VEmFluctuationModel* fm = 0;
    145145  modelManager->AddEmModel(order, p, fm, region);
    146   if(p) p->SetParticleChange(pParticleChange);
     146  if(p) { p->SetParticleChange(pParticleChange); }
    147147}
    148148//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    170170{
    171171  return modelManager->GetModel(idx, ver);
    172 }
    173 
    174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    175 
    176 void G4VMultipleScattering::BuildPhysicsTable(const G4ParticleDefinition& part)
    177 {
    178   G4String num = part.GetParticleName();
    179   if(1 < verboseLevel) {
    180     G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for "
    181            << GetProcessName()
    182            << " and particle " << num
    183            << G4endl;
    184   }
    185 
    186   if (buildLambdaTable && firstParticle == &part) {
    187 
    188     const G4ProductionCutsTable* theCoupleTable=
    189           G4ProductionCutsTable::GetProductionCutsTable();
    190     size_t numOfCouples = theCoupleTable->GetTableSize();
    191 
    192     G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
    193 
    194     G4PhysicsLogVector* aVector = 0;
    195     G4PhysicsLogVector* bVector = 0;
    196 
    197     for (size_t i=0; i<numOfCouples; ++i) {
    198 
    199       if (theLambdaTable->GetFlag(i)) {
    200         // create physics vector and fill it
    201         const G4MaterialCutsCouple* couple =
    202           theCoupleTable->GetMaterialCutsCouple(i);
    203         if(!bVector) {
    204           aVector = static_cast<G4PhysicsLogVector*>(PhysicsVector(couple));
    205           bVector = aVector;
    206         } else {
    207           aVector = new G4PhysicsLogVector(*bVector);
    208         }       
    209         //G4PhysicsVector* aVector = PhysicsVector(couple);
    210         aVector->SetSpline(splineFlag);
    211         modelManager->FillLambdaVector(aVector, couple, false);
    212         if(splineFlag) aVector->FillSecondDerivatives();
    213         G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
    214       }
    215     }
    216 
    217     if(1 < verboseLevel) {
    218       G4cout << "Lambda table is built for "
    219              << num
    220              << G4endl;
    221     }
    222   }
    223   if(verboseLevel>0 && ( num == "e-" || num == "mu+" || 
    224                          num == "proton" || num == "pi-" ||
    225                          num == "GenericIon")) {
    226     PrintInfoDefinition();
    227     if(2 < verboseLevel && theLambdaTable) G4cout << *theLambdaTable << G4endl;
    228   }
    229 
    230   if(1 < verboseLevel) {
    231     G4cout << "### G4VMultipleScattering::BuildPhysicsTable() done for "
    232            << GetProcessName()
    233            << " and particle " << num
    234            << G4endl;
    235   }
    236172}
    237173
    … …  
    260196  }
    261197
     198  (G4LossTableManager::Instance())->PreparePhysicsTable(&part, this);
     199
    262200  if(1 < verboseLevel) {
    263201    G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for "
    … …  
    267205           << G4endl;
    268206  }
    269 
    270   (G4LossTableManager::Instance())->EmConfigurator()->AddModels();
    271207
    272208  if(firstParticle == &part) {
    … …  
    307243//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    308244
     245void G4VMultipleScattering::BuildPhysicsTable(const G4ParticleDefinition& part)
     246{
     247  G4String num = part.GetParticleName();
     248  if(1 < verboseLevel) {
     249    G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for "
     250           << GetProcessName()
     251           << " and particle " << num
     252           << G4endl;
     253  }
     254
     255  (G4LossTableManager::Instance())->BuildPhysicsTable(firstParticle);
     256
     257  if (buildLambdaTable && firstParticle == &part) {
     258
     259    const G4ProductionCutsTable* theCoupleTable=
     260          G4ProductionCutsTable::GetProductionCutsTable();
     261    size_t numOfCouples = theCoupleTable->GetTableSize();
     262
     263    G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
     264
     265    G4PhysicsLogVector* aVector = 0;
     266    G4PhysicsLogVector* bVector = 0;
     267
     268    for (size_t i=0; i<numOfCouples; ++i) {
     269
     270      if (theLambdaTable->GetFlag(i)) {
     271        // create physics vector and fill it
     272        const G4MaterialCutsCouple* couple =
     273          theCoupleTable->GetMaterialCutsCouple(i);
     274        if(!bVector) {
     275          aVector = static_cast<G4PhysicsLogVector*>(PhysicsVector(couple));
     276          bVector = aVector;
     277        } else {
     278          aVector = new G4PhysicsLogVector(*bVector);
     279        }       
     280        //G4PhysicsVector* aVector = PhysicsVector(couple);
     281        aVector->SetSpline(splineFlag);
     282        modelManager->FillLambdaVector(aVector, couple, false);
     283        if(splineFlag) aVector->FillSecondDerivatives();
     284        G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
     285      }
     286    }
     287
     288    if(1 < verboseLevel) {
     289      G4cout << "Lambda table is built for "
     290             << num
     291             << G4endl;
     292    }
     293  }
     294  if(verboseLevel>0 && ( num == "e-" || num == "mu+" || 
     295                         num == "proton" || num == "pi-" ||
     296                         num == "GenericIon")) {
     297    PrintInfoDefinition();
     298    if(2 < verboseLevel && theLambdaTable) G4cout << *theLambdaTable << G4endl;
     299  }
     300
     301  if(1 < verboseLevel) {
     302    G4cout << "### G4VMultipleScattering::BuildPhysicsTable() done for "
     303           << GetProcessName()
     304           << " and particle " << num
     305           << G4endl;
     306  }
     307}
     308
     309//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     310
    309311void G4VMultipleScattering::PrintInfoDefinition()
    310312{
    … …  
    362364//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    363365
     366G4double
     367G4VMultipleScattering::PostStepGetPhysicalInteractionLength(
     368              const G4Track&, G4double, G4ForceCondition* condition)
     369{
     370  *condition = Forced;
     371  return DBL_MAX;
     372}
     373
     374//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     375
     376G4VParticleChange*
     377G4VMultipleScattering::AlongStepDoIt(const G4Track& track, const G4Step& step)
     378{
     379  if(currentModel->IsActive(track.GetKineticEnergy())) {
     380    fParticleChange.ProposeTrueStepLength(currentModel->ComputeTrueStepLength(step.GetStepLength()));
     381  }
     382  return &fParticleChange;
     383}
     384
     385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     386
     387G4VParticleChange*
     388G4VMultipleScattering::PostStepDoIt(const G4Track& track, const G4Step& step)
     389{
     390  fParticleChange.Initialize(track);
     391  if(currentModel->IsActive(track.GetKineticEnergy())) {
     392    currentModel->SampleScattering(track.GetDynamicParticle(),
     393                                   step.GetPostStepPoint()->GetSafety());
     394  }
     395  return &fParticleChange;
     396}
     397
     398//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     399
    364400G4double G4VMultipleScattering::GetContinuousStepLimit(
    365401                                       const G4Track& track,
    … …  
    371407  return AlongStepGetPhysicalInteractionLength(track,previousStepSize,currentMinimalStep,
    372408                                               currentSafety, selection);
     409}
     410
     411//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     412
     413G4double G4VMultipleScattering::ContinuousStepLimit(
     414                                       const G4Track& track,
     415                                       G4double previousStepSize,
     416                                       G4double currentMinimalStep,
     417                                       G4double& currentSafety)
     418{
     419  return GetContinuousStepLimit(track,previousStepSize,currentMinimalStep,
     420                                currentSafety);
    373421}
    374422

  • trunk/source/processes/electromagnetic/utils/test/testG4EnergyLossTables.cc

    r1199 r1315  
    2626//
    2727// $Id: testG4EnergyLossTables.cc,v 1.7 2006/06/29 19:55:29 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-cand-01 $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//-------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.