Ignore:
Timestamp:
May 28, 2009, 4:26:57 PM (15 years ago)
Author:
garnier
Message:

maj sur la beta de geant 4.9.3

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

Legend:

Unmodified
Added
Removed

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

    r1007 r1055  
    1 $Id: History,v 1.364 2008/11/20 20:32:40 vnivanch Exp $
     1$Id: History,v 1.383 2009/05/27 11:55:02 vnivanch Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     2027 May 09: V.Ivant (emutils-V09-02-11)
     21- G4VMultipleScattering - discarded changes from the previous tag
     22
     2326 May 09: V.Ivant (emutils-V09-02-10)
     24- G4VEmModel: remove protection in  from previous tag and include protection
     25              to G4ParticleChangeForLoss;
     26              more save initialisation of G4EmElementSelectors
     27- G4VMultipleScattering - added a protection against zero kinetic energy 
     28- G4EmElementSelector - do not use spline
     29
     3022 May 09: V.Ivant (emutils-V09-02-09)
     31- Added protection into G4VEmModel for zero input energy
     32
     3315 May 09: V.Ivant (emutils-V09-02-08)
     34- Added new testG4EnergyLossTables and fixed GNUmakefile for tests
     35
     3610 May 09: V.Ivant (emutils-V09-02-07)
     37G4EmElementSelector - added protection for zero cross section at first and last
     38                      bins of physics vector
     39G4VMscModel, G4VMultipleScattering - set default Range Factor to 0.04
     40
     4117 April 09: V.Ivant (emutils-V09-02-06)
     42- G4EmModelManager - fixed energy range selection algorithm for the case of
     43                     a small intersection of model energy intervals
     44- G4VEnergyLossProcess, G4VEmProcess, G4VMultipleScattering - set high enegry
     45                     limit 10 TeV and number of bins 77
     46
     4708 April 09: V.Ivant (emutils-V09-02-05)
     48- G4LossTableManager - added G4EmConfigurator providing easier addition of
     49                       models per region
     50- G4VMultipleScattering, G4VEmProcess, G4VEnergyLossProcess: added
     51                       initialisation of the G4EmConfigurator
     52
     5308 April 09: V.Ivant (emutils-V09-02-04)
     54- G4EmModelManager - fixed energy range selection algorithm for the case
     55                     when there is an intersection of energy regions of standard
     56                     and low-energy models,
     57                   - reduce internal vectors if no model per region are initialized.
     58                   - do not initilise unused models.
     59- G4VEmModel - msc methods are moved to G4VMscModel, added protected
     60               methods for initialisation of ParticleChange             
     61- G4VMultipleScattering, G4VEmProcess, G4VEnergyLossProcess:
     62  methods for initialisations are moved from inline to source
     63
     6426 February 09: V.Ivant (emutils-V09-02-03)
     65G4EmConfigurator - fixed for the case if only fluctuation model is set
     66                   and main model is default
     67
     6822 February 09: V.Ivant (emutils-V09-02-02)
     69- G4VEmModel - make methods to access geometry protected, added new
     70               method SetSampleZ, added geommax private member
     71- G4EmCalculator - added possibility to be used by DNA processes:
     72                   take into account special DNA particles
     73
     7418 February 09: V.Ivant (emutils-V09-02-01)
     75G4VEmModel, G4VEmFluctuationModel, G4VEnegryLossProcess, G4VEmProcess,
     76G4VMultipleScattering - move all virtual methods to source, update comments
     77G4VEmModel - added flagDeexcitation and Get/Set methods
     78G4VEnegryLossProcess, G4VEmProcess - added calls to deexcitation PostStep
     79G4EmProcessOptions - added ActivateDeexcitation method
     80G4EnergyLossMessenger - added /process/em/deexcitation UI command
     81G4LossTableBuilder - added protection in BuildRangeTable against zero dedx
     82
     8327 January 09: V.Ivant (emutils-V09-02-00)
     84G4VEmModel - added method SampleDeexcitationAlongStep
     85G4VEnegryLossProcess - added deexcitation AlongStep per region
     86G4VMscModel - added methdos: InitialiseSafetyHelper, ComputeSafety,
     87              ComputeGeomLimit, ComputeDisplacement
     88G4VEmProcess - added possibility to set more than 1 model
    1989
    209020 November 08: V.Ivant (emutils-V09-01-37)

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DummyModel.hh,v 1.3 2007/05/22 17:31:57 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4DummyModel.hh,v 1.4 2009/04/07 18:39:47 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5151
    5252#include "globals.hh"
    53 #include "G4VEmModel.hh"
     53#include "G4VMscModel.hh"
    5454
    55 class G4DummyModel :  public G4VEmModel
     55class G4DummyModel :  public G4VMscModel
    5656{
    5757

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmElementSelector.hh,v 1.2 2008/07/22 15:55:15 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4EmElementSelector.hh,v 1.4 2009/05/10 18:26:43 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8989  G4VEmModel*       model;
    9090  const G4Material* material;
     91  const G4Element*  element;
    9192  const G4ElementVector* theElementVector;
    9293
    … …  
    107108inline const G4Element* G4EmElementSelector::SelectRandomAtom(G4double e)
    108109{
    109   const G4Element* elm = (*theElementVector)[nElmMinusOne];
    110110  if (nElmMinusOne > 0) {
    111111    G4bool b;
    … …  
    113113    for(G4int i=0; i<nElmMinusOne; i++) {
    114114      if (x <= (xSections[i])->GetValue(e,b)) {
    115         elm = (*theElementVector)[i];
     115        element = (*theElementVector)[i];
    116116        break;
    117117      }
    118118    }
    119119  }
    120   return elm;
     120  return element;
    121121}
    122122

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmModelManager.hh,v 1.25 2008/10/13 14:56:56 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4EmModelManager.hh,v 1.28 2009/04/09 15:53:17 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5353// 13-05-06 Add GetModel by index method (VI)
    5454// 15-03-07 Add maxCutInRange (V.Ivanchenko)
     55// 08-04-08 Simplify Select method for only one G4RegionModel (VI)
    5556//
    5657// Class Description:
    … …  
    9192  ~G4RegionModels();
    9293
    93   G4int SelectIndex(G4double e) const {
     94  inline G4int SelectIndex(G4double e) const {
    9495    G4int idx = 0;
    9596    if (nModelsForRegion>1) {
    … …  
    100101  };
    101102
    102   G4int ModelIndex(G4int n) const {
     103  inline G4int ModelIndex(G4int n) const {
    103104    return theListOfModelIndexes[n];
    104105  };
    105106
    106   G4int NumberOfModels() const {
     107  inline G4int NumberOfModels() const {
    107108    return nModelsForRegion;
    108109  };
    109110
    110   G4double LowEdgeEnergy(G4int n) const {
     111  inline G4double LowEdgeEnergy(G4int n) const {
    111112    return lowKineticEnergy[n];
    112113  };
    113114
    114   const G4Region* Region() const {
     115  inline const G4Region* Region() const {
    115116    return theRegion;
    116117  };
    … …  
    150151                                       G4int);
    151152
    152   const G4DataVector* Cuts() const;
    153 
    154   const G4DataVector* SubCutoff() const;
    155 
    156153  void FillDEDXVector(G4PhysicsVector*, const G4MaterialCutsCouple*,
    157154                      G4EmTableType t = fRestricted);
    … …  
    160157                        G4bool startFromNull = true, G4EmTableType t = fRestricted);
    161158
    162   G4VEmModel* SelectModel(G4double& energy, size_t& index);
    163 
    164159  G4VEmModel* GetModel(G4int, G4bool ver = false);
    165 
    166   G4int NumberOfModels() const;
    167160
    168161  void AddEmModel(G4int, G4VEmModel*, G4VEmFluctuationModel*, const G4Region*);
    … …  
    171164
    172165  void DumpModelList(G4int verb);
     166
     167  inline G4VEmModel* SelectModel(G4double& energy, size_t& index);
     168
     169  inline const G4DataVector* Cuts() const;
     170
     171  inline const G4DataVector* SubCutoff() const;
     172
     173  inline G4int NumberOfModels() const;
    173174
    174175private:
    … …  
    190191  std::vector<const G4Region*>            regions;
    191192  std::vector<G4int>                      orderOfModels;
     193  std::vector<G4int>                      isUsed;
    192194
    193195  G4int                       nEmModels;
    … …  
    209211  G4int                       verboseLevel;
    210212
    211   // cash
    212   G4int                       currentIdx;
    213 
     213  // may be changed in run time
     214  G4RegionModels*             currRegionModel;
    214215};
    215216
    … …  
    220221                                                 size_t& index)
    221222{
    222   currentIdx =
    223     (setOfRegionModels[idxOfRegionModels[index]])->SelectIndex(kinEnergy);
    224   return models[currentIdx];
     223  if(nRegions > 1) currRegionModel = setOfRegionModels[idxOfRegionModels[index]];
     224  return models[currRegionModel->SelectIndex(kinEnergy)];
    225225}
    226226

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmProcessOptions.hh,v 1.14 2008/04/17 10:33:26 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4EmProcessOptions.hh,v 1.15 2009/02/18 14:40:10 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929//
    … …  
    107107  void SetLinearLossLimit(G4double val);
    108108
    109   void ActivateDeexcitation(G4bool val, const G4Region* r = 0);
     109  void ActivateDeexcitation(const G4String& proc, G4bool val,
     110                            const G4String& reg = "");
    110111
    111112  void SetMscStepLimitation(G4MscStepLimitType val);

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

    r1007 r1055  
    2525//
    2626//
    27 // $Id: G4EnergyLossMessenger.hh,v 1.22 2008/10/20 13:27:45 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4EnergyLossMessenger.hh,v 1.23 2009/02/18 14:40:10 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    107107  G4UIcmdWithADouble*        MinSubSecCmd;
    108108  G4UIcommand*               StepFuncCmd;
     109  G4UIcommand*               deexCmd;
    109110  G4UIcmdWithAString*        mscCmd;
    110111  G4UIcmdWithADoubleAndUnit* MinEnCmd;

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LossTableManager.hh,v 1.53 2008/07/15 16:56:38 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4LossTableManager.hh,v 1.54 2009/04/09 16:10:57 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929//
    … …  
    9191class G4EmCorrections;
    9292class G4EmSaturation;
     93class G4EmConfigurator;
    9394class G4LossTableBuilder;
    9495
    … …  
    232233
    233234  G4EmSaturation* EmSaturation();
     235
     236  G4EmConfigurator* EmConfigurator();
    234237
    235238private:
    … …  
    304307  G4EmCorrections*            emCorrections;
    305308  G4EmSaturation*             emSaturation;
     309  G4EmConfigurator*           emConfigurator;
    306310
    307311  const G4ParticleDefinition* firstParticle;

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmFluctuationModel.hh,v 1.11 2008/09/12 14:47:38 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VEmFluctuationModel.hh,v 1.12 2009/02/19 11:25:50 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    101101  //------------------------------------------------------------------------
    102102
    103   G4String GetName() const;
     103  inline G4String GetName() const;
    104104
    105105private:
    … …  
    115115//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    116116
    117 inline void G4VEmFluctuationModel::InitialiseMe(const G4ParticleDefinition*)
    118 {}
    119 
    120 inline
    121 void G4VEmFluctuationModel::SetParticleAndCharge(const G4ParticleDefinition*,
    122                                                  G4double)
    123 {}
    124 
    125117inline G4String G4VEmFluctuationModel::GetName() const
    126118{

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmModel.hh,v 1.59 2008/11/13 19:29:41 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VEmModel.hh,v 1.69 2009/05/26 15:00:49 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6565// 12-09-08 Added methods GetParticleCharge, GetChargeSquareRatio,
    6666//          CorrectionsAlongStep, ActivateNuclearStopping (VI)
     67// 16-02-09 Moved implementations of virtual methods to source (VI)
     68// 07-04-09 Moved msc methods from G4VEmModel to G4VMscModel (VI)
    6769//
    6870// Class Description:
    … …  
    9294class G4Region;
    9395class G4VParticleChange;
     96class G4ParticleChangeForLoss;
     97class G4ParticleChangeForGamma;
    9498class G4Track;
    9599
    … …  
    120124  //------------------------------------------------------------------------
    121125
    122   // dEdx per unit length
    123   virtual G4double ComputeDEDX(const G4MaterialCutsCouple*,
    124                                const G4ParticleDefinition*,
    125                                G4double kineticEnergy,
    126                                G4double cutEnergy = DBL_MAX);
    127 
    128126  // main method to compute dEdx
    129127  virtual G4double ComputeDEDXPerVolume(const G4Material*,
    … …  
    131129                                        G4double kineticEnergy,
    132130                                        G4double cutEnergy = DBL_MAX);
    133 
    134   // cross section per volume
    135   virtual G4double CrossSection(const G4MaterialCutsCouple*,
    136                                 const G4ParticleDefinition*,
    137                                 G4double kineticEnergy,
    138                                 G4double cutEnergy = 0.0,
    139                                 G4double maxEnergy = DBL_MAX);
    140131
    141132  // main method to compute cross section per Volume
    … …  
    175166                                    G4double length);
    176167
     168  // sample PIXE deexcitation
     169  virtual void SampleDeexcitationAlongStep(const G4Material*,
     170                                           const G4Track&,
     171                                           G4double& eloss);
     172
     173  // add region for the model
     174  virtual void DefineForRegion(const G4Region*);
     175
     176  // initilisation at run time for given material
     177  virtual void SetupForMaterial(const G4ParticleDefinition*,
     178                                const G4Material*,
     179                                G4double kineticEnergy);
     180
    177181protected:
     182
     183  // initialisation of the ParticleChange for the model
     184  G4ParticleChangeForLoss* GetParticleChangeForLoss();
     185
     186  // initialisation of the ParticleChange for the model
     187  G4ParticleChangeForGamma* GetParticleChangeForGamma();
    178188
    179189  // kinematically allowed max kinetic energy of a secondary
    180190  virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
    181191                                      G4double kineticEnergy); 
    182                                          
    183   //------------------------------------------------------------------------
    184   // Methods for msc simulation which needs to be overwritten
    185   //------------------------------------------------------------------------
    186192
    187193public:
    188 
    189   virtual void SampleScattering(const G4DynamicParticle*,
    190                                 G4double safety);
    191 
    192   virtual G4double ComputeTruePathLengthLimit(const G4Track& track,
    193                                               G4PhysicsTable* theLambdaTable,
    194                                               G4double currentMinimalStep);
    195 
    196   virtual G4double ComputeGeomPathLength(G4double truePathLength);
    197 
    198   virtual G4double ComputeTrueStepLength(G4double geomPathLength);
    199 
    200   virtual void DefineForRegion(const G4Region*);
    201 
    202   virtual void SetupForMaterial(const G4ParticleDefinition*,
    203                                 const G4Material*,
    204                                 G4double kineticEnergy);
    205194
    206195  //------------------------------------------------------------------------
    … …  
    212201                                  const G4DataVector&);
    213202
     203  // dEdx per unit length
     204  inline G4double ComputeDEDX(const G4MaterialCutsCouple*,
     205                              const G4ParticleDefinition*,
     206                              G4double kineticEnergy,
     207                              G4double cutEnergy = DBL_MAX);
     208
     209  // cross section per volume
     210  inline G4double CrossSection(const G4MaterialCutsCouple*,
     211                                const G4ParticleDefinition*,
     212                                G4double kineticEnergy,
     213                                G4double cutEnergy = 0.0,
     214                                G4double maxEnergy = DBL_MAX);
     215
    214216  // compute mean free path via cross section per volume
    215   G4double ComputeMeanFreePath(const G4ParticleDefinition*,
    216                                G4double kineticEnergy,
    217                                const G4Material*,   
    218                                G4double cutEnergy = 0.0,
    219                                G4double maxEnergy = DBL_MAX);
     217  inline G4double ComputeMeanFreePath(const G4ParticleDefinition*,
     218                                      G4double kineticEnergy,
     219                                      const G4Material*,   
     220                                      G4double cutEnergy = 0.0,
     221                                      G4double maxEnergy = DBL_MAX);
    220222
    221223  // generic cross section per element
    … …  
    233235                                           G4double maxEnergy = DBL_MAX);
    234236
    235   // this method can be used only in the case if generic method to compute
    236   // cross section per volume is used and not overwritten in derived class
     237  // to select atom cross section per volume is recomputed for each element
    237238  inline const G4Element* SelectRandomAtom(const G4Material*,
    238239                                           const G4ParticleDefinition*,
    … …  
    260261  inline G4bool LPMFlag() const;
    261262
     263  inline G4bool DeexcitationFlag() const;
     264
    262265  inline void SetHighEnergyLimit(G4double);
    263266
    … …  
    270273  inline void SetLPMFlag(G4bool val);
    271274
     275  inline void SetDeexcitationFlag(G4bool val);
     276
    272277  inline void ActivateNuclearStopping(G4bool);
    273278
    … …  
    278283  inline void SetParticleChange(G4VParticleChange*, G4VEmFluctuationModel*);
    279284
     285  inline void SetCurrentCouple(const G4MaterialCutsCouple*);
     286
    280287protected:
    281288
     289  inline const G4MaterialCutsCouple* CurrentCouple() const;
     290
     291  inline void SetCurrentElement(const G4Element*);
     292
    282293  inline const G4Element* GetCurrentElement() const;
    283 
    284   inline void SetCurrentElement(const G4Element*);
    285294
    286295private:
    … …  
    315324private:
    316325
    317   const G4Element* currentElement;
     326  const G4MaterialCutsCouple* currentCouple;
     327  const G4Element*            currentElement;
     328
    318329  G4int                  nsec;
     330  G4bool                 flagDeexcitation;
    319331  std::vector<G4double>  xsec;
    320332
    … …  
    324336//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    325337
    326 inline G4double G4VEmModel::HighEnergyLimit() const
    327 {
    328   return highLimit;
    329 }
    330 
    331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    332 
    333 inline G4double G4VEmModel::LowEnergyLimit() const
    334 {
    335   return lowLimit;
    336 }
    337 
    338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    339 
    340 inline G4double G4VEmModel::PolarAngleLimit() const
    341 {
    342   return polarAngleLimit;
    343 }
    344 
    345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    346 
    347 inline G4double G4VEmModel::SecondaryThreshold() const
    348 {
    349   return secondaryThreshold;
    350 }
    351 
    352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    353 
    354 inline G4bool G4VEmModel::LPMFlag() const
    355 {
    356   return theLPMflag;
    357 }
    358 
    359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    360 
    361 inline void G4VEmModel::SetHighEnergyLimit(G4double val)
    362 {
    363   highLimit = val;
    364 }
    365 
    366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    367 
    368 inline void G4VEmModel::SetLowEnergyLimit(G4double val)
    369 {
    370   lowLimit = val;
    371 }
    372 
    373 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    374 
    375 inline void G4VEmModel::SetPolarAngleLimit(G4double val)
    376 {
    377   polarAngleLimit = val;
    378 }
    379 
    380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    381 
    382 inline void G4VEmModel::SetSecondaryThreshold(G4double val)
    383 {
    384   secondaryThreshold = val;
    385 }
    386 
    387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    388 
    389 inline void G4VEmModel::SetLPMFlag(G4bool val)
    390 {
    391   theLPMflag = val;
    392 }
    393 
    394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    395 
    396 inline void G4VEmModel::ActivateNuclearStopping(G4bool val)
    397 {
    398   nuclearStopping = val;
     338inline G4double G4VEmModel::ComputeDEDX(const G4MaterialCutsCouple* c,
     339                                        const G4ParticleDefinition* p,
     340                                        G4double kinEnergy,
     341                                        G4double cutEnergy)
     342{
     343  currentCouple = c;
     344  return ComputeDEDXPerVolume(c->GetMaterial(),p,kinEnergy,cutEnergy);
     345}
     346
     347//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     348
     349inline G4double G4VEmModel::CrossSection(const G4MaterialCutsCouple* c,
     350                                         const G4ParticleDefinition* p,
     351                                         G4double kinEnergy,
     352                                         G4double cutEnergy,
     353                                         G4double maxEnergy)
     354{
     355  currentCouple = c;
     356  return CrossSectionPerVolume(c->GetMaterial(),p,kinEnergy,cutEnergy,maxEnergy);
     357}
     358
     359//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     360
     361inline G4double G4VEmModel::ComputeMeanFreePath(const G4ParticleDefinition* p,
     362                                                G4double ekin,
     363                                                const G4Material* material,     
     364                                                G4double emin,
     365                                                G4double emax)
     366{
     367  G4double mfp = DBL_MAX;
     368  G4double cross = CrossSectionPerVolume(material,p,ekin,emin,emax);
     369  if (cross > DBL_MIN) mfp = 1./cross;
     370  return mfp;
    399371}
    400372
    … …  
    415387//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    416388
    417 inline void G4VEmModel::SetParticleChange(G4VParticleChange* p, 
    418                                           G4VEmFluctuationModel* f = 0)
    419 {
    420   if(p && pParticleChange != p) pParticleChange = p;
    421   fluc = f;
    422 }
    423 
    424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    425 
    426 
    427 inline G4VEmFluctuationModel* G4VEmModel::GetModelOfFluctuations()
    428 {
    429   return fluc;
    430 }
    431 
    432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    433 
    434 inline G4double G4VEmModel::MinEnergyCut(const G4ParticleDefinition*,
    435                                          const G4MaterialCutsCouple*)
    436 {
    437   return 0.0;
    438 }
    439 
    440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    441 
    442 inline G4double G4VEmModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
    443                                                  const G4Material*, G4double)
    444 {
    445   G4double q = p->GetPDGCharge()/CLHEP::eplus;
    446   return q*q;
    447 }
    448 
    449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    450 
    451 inline G4double G4VEmModel::GetParticleCharge(const G4ParticleDefinition* p,
    452                                               const G4Material*, G4double)
    453 {
    454   return p->GetPDGCharge();
    455 }
    456 
    457 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    458 
    459 inline void G4VEmModel::CorrectionsAlongStep(const G4MaterialCutsCouple*,
    460                                              const G4DynamicParticle*,
    461                                              G4double&,G4double&,G4double)
    462 {}
    463 
    464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    465 
    466 inline G4double G4VEmModel::ComputeDEDXPerVolume(const G4Material*,
    467                                                  const G4ParticleDefinition*,
    468                                                  G4double,G4double)
    469 {
    470   return 0.0;
    471 }
    472 
    473 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    474 
    475 inline G4double G4VEmModel::ComputeDEDX(const G4MaterialCutsCouple* c,
    476                                         const G4ParticleDefinition* p,
    477                                         G4double kinEnergy,
    478                                         G4double cutEnergy)
    479 {
    480   return ComputeDEDXPerVolume(c->GetMaterial(),p,kinEnergy,cutEnergy);
    481 }
    482 
    483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    484 
    485 inline G4double G4VEmModel::CrossSection(const G4MaterialCutsCouple* c,
    486                                          const G4ParticleDefinition* p,
    487                                          G4double kinEnergy,
    488                                          G4double cutEnergy,
    489                                          G4double maxEnergy)
    490 {
    491   return CrossSectionPerVolume(c->GetMaterial(),p,
    492                                kinEnergy,cutEnergy,maxEnergy);
    493 }
    494 
    495 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    496 
    497 inline G4double G4VEmModel::ComputeCrossSectionPerAtom(
    498                                          const G4ParticleDefinition*,
    499                                          G4double, G4double, G4double,
    500                                          G4double, G4double)
    501 {
    502   return 0.0;
    503 }
    504 
    505 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    506 
    507389inline
    508390const G4Element* G4VEmModel::SelectRandomAtom(const G4MaterialCutsCouple* couple,
    … …  
    512394                                              G4double maxEnergy)
    513395{
     396  currentCouple = couple;
    514397  if(nSelectors > 0) {
    515398    currentElement =
    … …  
    571454//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    572455
    573 inline const G4Element* G4VEmModel::GetCurrentElement() const
    574 {
    575   return currentElement;
    576 }
    577 
    578 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    579 
    580 inline void G4VEmModel::SetCurrentElement(const G4Element* elm)
    581 {
    582   currentElement = elm;
     456inline G4VEmFluctuationModel* G4VEmModel::GetModelOfFluctuations()
     457{
     458  return fluc;
     459}
     460
     461//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     462
     463inline G4double G4VEmModel::HighEnergyLimit() const
     464{
     465  return highLimit;
     466}
     467
     468//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     469
     470inline G4double G4VEmModel::LowEnergyLimit() const
     471{
     472  return lowLimit;
     473}
     474
     475//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     476
     477inline G4double G4VEmModel::PolarAngleLimit() const
     478{
     479  return polarAngleLimit;
     480}
     481
     482//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     483
     484inline G4double G4VEmModel::SecondaryThreshold() const
     485{
     486  return secondaryThreshold;
     487}
     488
     489//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     490
     491inline G4bool G4VEmModel::LPMFlag() const
     492{
     493  return theLPMflag;
     494}
     495
     496//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     497
     498inline G4bool G4VEmModel::DeexcitationFlag() const
     499{
     500  return flagDeexcitation;
     501}
     502
     503//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     504
     505inline void G4VEmModel::SetHighEnergyLimit(G4double val)
     506{
     507  highLimit = val;
     508}
     509
     510//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     511
     512inline void G4VEmModel::SetLowEnergyLimit(G4double val)
     513{
     514  lowLimit = val;
     515}
     516
     517//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     518
     519inline void G4VEmModel::SetPolarAngleLimit(G4double val)
     520{
     521  polarAngleLimit = val;
     522}
     523
     524//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     525
     526inline void G4VEmModel::SetSecondaryThreshold(G4double val)
     527{
     528  secondaryThreshold = val;
     529}
     530
     531//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     532
     533inline void G4VEmModel::SetLPMFlag(G4bool val)
     534{
     535  theLPMflag = val;
     536}
     537
     538//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     539
     540inline void G4VEmModel::SetDeexcitationFlag(G4bool val)
     541{
     542  flagDeexcitation = val;
     543}
     544
     545//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     546
     547inline void G4VEmModel::ActivateNuclearStopping(G4bool val)
     548{
     549  nuclearStopping = val;
    583550}
    584551
    … …  
    594561//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    595562
    596 inline G4double G4VEmModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
    597                                                G4double kineticEnergy)
    598 {
    599   return kineticEnergy;
    600 }
    601 
    602 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    603 
    604563inline const G4String& G4VEmModel::GetName() const
    605564{
    … …  
    608567
    609568//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    610 // Methods for msc simulation
    611 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    612 
    613 inline void G4VEmModel::SampleScattering(const G4DynamicParticle*, G4double)
    614 {}
    615 
    616 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    617 
    618 inline G4double G4VEmModel::ComputeTruePathLengthLimit(
    619                                 const G4Track&,
    620                                 G4PhysicsTable*,
    621                                 G4double)
    622 {
    623   return DBL_MAX;
    624 }
    625 
    626 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    627 
    628 inline G4double G4VEmModel::ComputeGeomPathLength(G4double truePathLength)
    629 {
    630   return truePathLength;
    631 }
    632 
    633 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    634 
    635 inline G4double G4VEmModel::ComputeTrueStepLength(G4double geomPathLength)
    636 {
    637   return geomPathLength;
    638 }
    639 
    640 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    641 
    642 inline void G4VEmModel::DefineForRegion(const G4Region*)
    643 {}
    644 
    645 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    646 
    647 inline void G4VEmModel::SetupForMaterial(const G4ParticleDefinition*,
    648                                          const G4Material*, G4double)
    649 {}
     569
     570inline void G4VEmModel::SetParticleChange(G4VParticleChange* p, 
     571                                          G4VEmFluctuationModel* f = 0)
     572{
     573  if(p && pParticleChange != p) pParticleChange = p;
     574  fluc = f;
     575}
     576
     577//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     578
     579inline void G4VEmModel::SetCurrentCouple(const G4MaterialCutsCouple* p)
     580{
     581  currentCouple = p;
     582}
     583
     584//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     585
     586inline const G4MaterialCutsCouple* G4VEmModel::CurrentCouple() const
     587{
     588  return currentCouple;
     589}
     590
     591//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     592
     593inline void G4VEmModel::SetCurrentElement(const G4Element* elm)
     594{
     595  currentElement = elm;
     596}
     597
     598//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     599
     600inline const G4Element* G4VEmModel::GetCurrentElement() const
     601{
     602  return currentElement;
     603}
    650604
    651605//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmProcess.hh,v 1.47 2008/07/31 13:01:26 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VEmProcess.hh,v 1.51 2009/04/07 18:39:47 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    108108
    109109  //------------------------------------------------------------------------
    110   // Methods with standard implementation; may be overwritten if needed
    111   //------------------------------------------------------------------------
    112 
    113   inline G4double RecalculateLambda(G4double kinEnergy,
    114                                     const G4MaterialCutsCouple* couple);
    115 
    116   //------------------------------------------------------------------------
    117   // Generic methods common to all Discrete processes
     110  // Implementation of virtual methods common to all Discrete processes
    118111  //------------------------------------------------------------------------
    119112
    120113public:
    121 
    122   void PrintInfoDefinition();
    123 
    124   G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&);
    125114
    126115  // Initialise for build of tables
    … …  
    129118  // Build physics table during initialisation
    130119  void BuildPhysicsTable(const G4ParticleDefinition&);
     120
     121  void PrintInfoDefinition();
     122
     123  // implementation of virtual method, specific for G4VEmProcess
     124  G4double PostStepGetPhysicalInteractionLength(
     125                             const G4Track& track,
     126                             G4double   previousStepSize,
     127                             G4ForceCondition* condition
     128                            );
     129
     130  // implementation of virtual method, specific for G4VEmProcess
     131  G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&);
    131132
    132133  // Store PhysicsTable in a file.
    … …  
    145146                              G4bool ascii);
    146147
     148  // deexcitation activated per G4Region
     149  void ActivateDeexcitation(G4bool, const G4Region* r = 0);
     150
    147151  //------------------------------------------------------------------------
    148152  // Specific methods for Discrete EM post step simulation
    … …  
    152156  G4double CrossSectionPerVolume(G4double kineticEnergy,
    153157                                 const G4MaterialCutsCouple* couple);
    154 
    155   // implementation of virtual method
    156   virtual G4double PostStepGetPhysicalInteractionLength(
    157                              const G4Track& track,
    158                              G4double   previousStepSize,
    159                              G4ForceCondition* condition
    160                             );
    161158
    162159  // It returns the cross section of the process per atom
    … …  
    167164  inline G4double MeanFreePath(const G4Track& track);
    168165
    169   inline G4VEmModel* SelectModelForMaterial(G4double kinEnergy,
    170                                             size_t& idxRegion) const;
    171 
    172166  // It returns cross section per volume
    173167  inline G4double GetLambda(G4double& kinEnergy,
    … …  
    203197
    204198  //------------------------------------------------------------------------
    205   // Specific methods to set, access, modify models
    206   //------------------------------------------------------------------------
    207 
    208   // Add EM model coupled for the region
    209   inline void AddEmModel(G4int, G4VEmModel*, const G4Region* region = 0);
     199  // Specific methods to set, access, modify models and basic parameters
     200  //------------------------------------------------------------------------
     201
     202protected:
     203  // Select model in run time
     204  inline void SelectModel(G4double& kinEnergy);
     205
     206public:
     207  // Select model by energy and region index
     208  inline G4VEmModel* SelectModelForMaterial(G4double kinEnergy,
     209                                            size_t& idxRegion) const;
    210210   
     211  // Add model for region, smaller value of order defines which
     212  // model will be selected for a given energy interval 
     213  void AddEmModel(G4int, G4VEmModel*, const G4Region* region = 0);
     214
    211215  // Assign a model to a process
    212   inline void SetModel(G4VEmModel*);
     216  void SetModel(G4VEmModel*, G4int index = 1);
    213217 
    214218  // return the assigned model
    215   inline G4VEmModel* Model();
     219  G4VEmModel* Model(G4int index = 1);
    216220   
    217221  // Define new energy range for the model identified by the name
    218   inline void UpdateEmModel(const G4String&, G4double, G4double);
     222  void UpdateEmModel(const G4String&, G4double, G4double);
    219223
    220224  // Access to models
    221   inline G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false);
    222 
    223   //------------------------------------------------------------------------
    224   // Get/set parameters used for simulation of energy loss
    225   //------------------------------------------------------------------------
    226 
    227   inline void ActivateDeexcitation(G4bool, const G4Region* r = 0);
     225  G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false);
    228226
    229227  inline void SetLambdaFactor(G4double val);
    … …  
    233231
    234232  inline void SetApplyCuts(G4bool val);
     233
     234  //------------------------------------------------------------------------
     235  // Other generic methods
     236  //------------------------------------------------------------------------
    235237 
    236238protected:
    … …  
    242244  G4PhysicsVector* LambdaPhysicsVector(const G4MaterialCutsCouple*);
    243245
     246  inline G4double RecalculateLambda(G4double kinEnergy,
     247                                    const G4MaterialCutsCouple* couple);
     248
    244249  inline G4ParticleChangeForGamma* GetParticleChange();
    245250
    … …  
    248253  inline void SetSecondaryParticle(const G4ParticleDefinition* p);
    249254
    250   inline G4VEmModel* SelectModel(G4double& kinEnergy);
    251 
    252255  inline size_t CurrentMaterialCutsCoupleIndex() const;
    253256
    … …  
    280283  inline G4double ComputeCurrentLambda(G4double kinEnergy);
    281284
    282   // hide  assignment operator
    283 
     285  // copy constructor and hide assignment operator
    284286  G4VEmProcess(G4VEmProcess &);
    285287  G4VEmProcess & operator=(const G4VEmProcess &right);
    … …  
    297299  // ======== Parameters of the class fixed at initialisation =======
    298300
     301  std::vector<G4VEmModel*>     emModels;
     302
    299303  // tables and vectors
    300304  G4PhysicsTable*              theLambdaTable;
    … …  
    317321  G4bool                       applyCuts;
    318322  G4bool                       startFromNull;
    319 
    320   G4int                        nRegions;
    321   std::vector<G4Region*>       regions;
    322   std::vector<G4bool>          flagsDeexcitation;
     323  G4bool                       useDeexcitation;
     324
     325  G4int                        nDERegions;
     326  std::vector<const G4Region*> deRegions;
     327  G4bool*                      idxDERegions;
    323328
    324329  // ======== Cashed values - may be state dependent ================
    … …  
    332337  std::vector<G4DynamicParticle*> secParticles;
    333338
    334   G4VEmModel*                  selectedModel; 
     339  G4VEmModel*                  currentModel; 
    335340
    336341  const G4ParticleDefinition*  particle;
    … …  
    348353
    349354//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     355//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     356
     357inline G4double G4VEmProcess::ComputeCrossSectionPerAtom(
     358                 G4double kineticEnergy, G4double Z, G4double A, G4double cut)
     359{
     360  SelectModel(kineticEnergy);
     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
     371inline 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
     382inline 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
     391inline void G4VEmProcess::SetLambdaBinning(G4int nbins)
     392{
     393  nLambdaBins = nbins;
     394}
     395
     396//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     397
     398inline G4int G4VEmProcess::LambdaBinning() const
     399{
     400  return nLambdaBins;
     401}
     402
     403//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     404
     405inline void G4VEmProcess::SetMinKinEnergy(G4double e)
     406{
     407  minKinEnergy = e;
     408}
     409
     410//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     411
     412inline G4double G4VEmProcess::MinKinEnergy() const
     413{
     414  return minKinEnergy;
     415}
     416
     417//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     418
     419inline void G4VEmProcess::SetMaxKinEnergy(G4double e)
     420{
     421  maxKinEnergy = e;
     422}
     423
     424//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     425
     426inline G4double G4VEmProcess::MaxKinEnergy() const
     427{
     428  return maxKinEnergy;
     429}
     430
     431//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     432
     433inline 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
     442inline G4double G4VEmProcess::PolarAngleLimit() const
     443{
     444  return polarAngleLimit;
     445}
     446
     447//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     448
     449inline const G4PhysicsTable* G4VEmProcess::LambdaTable() const
     450{
     451  return theLambdaTable;
     452}
     453
     454//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     455
     456inline const G4ParticleDefinition* G4VEmProcess::Particle() const
     457{
     458  return particle;
     459}
     460
     461//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     462
     463inline const G4ParticleDefinition* G4VEmProcess::SecondaryParticle() const
     464{
     465  return secondaryParticle;
     466}
     467
     468//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     469
     470inline void G4VEmProcess::SelectModel(G4double& kinEnergy)
     471{
     472  currentModel = modelManager->SelectModel(kinEnergy, currentMaterialIndex);
     473  currentModel->SetCurrentCouple(currentCouple);
     474}
     475
     476//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     477
     478inline G4VEmModel* G4VEmProcess::SelectModelForMaterial(
     479                                   G4double kinEnergy, size_t& idxRegion) const
     480{
     481  return modelManager->SelectModel(kinEnergy, idxRegion);
     482}
     483
     484//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     485
     486inline void G4VEmProcess::SetLambdaFactor(G4double val)
     487{
     488  if(val > 0.0 && val <= 1.0) lambdaFactor = val;
     489}
     490
     491//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     492
     493inline void G4VEmProcess::SetIntegral(G4bool val)
     494{
     495  if(particle && particle != theGamma) integral = val;
     496  if(integral) buildLambdaTable = true;
     497}
     498
     499//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     500
     501inline G4bool G4VEmProcess::IsIntegral() const
     502{
     503  return integral;
     504}
     505
     506//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     507
     508inline void G4VEmProcess::SetApplyCuts(G4bool val)
     509{
     510  applyCuts = val;
     511}
     512
     513//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     514
     515inline G4double G4VEmProcess::RecalculateLambda(G4double e,
     516                                                const G4MaterialCutsCouple* couple)
     517{
     518  DefineMaterial(couple);
     519  return ComputeCurrentLambda(e);
     520}
     521
     522//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     523
     524inline G4ParticleChangeForGamma* G4VEmProcess::GetParticleChange()
     525{
     526  return &fParticleChange;
     527}
     528
     529//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     530
     531inline void G4VEmProcess::SetParticle(const G4ParticleDefinition* p)
     532{
     533  particle = p;
     534}
     535
     536//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     537
     538inline void G4VEmProcess::SetSecondaryParticle(const G4ParticleDefinition* p)
     539{
     540  secondaryParticle = p;
     541}
     542
     543//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     544
     545inline size_t G4VEmProcess::CurrentMaterialCutsCoupleIndex() const
     546{
     547  return currentMaterialIndex;
     548}
     549
     550//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     551
     552inline G4double G4VEmProcess::GetGammaEnergyCut()
     553{
     554  return (*theCutsGamma)[currentMaterialIndex];
     555}
     556
     557//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     558
     559inline G4double G4VEmProcess::GetElectronEnergyCut()
     560{
     561  return (*theCutsElectron)[currentMaterialIndex];
     562}
     563
     564//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     565
     566inline void G4VEmProcess::SetBuildTableFlag(G4bool val)
     567{
     568  buildLambdaTable = val;
     569  if(!val) integral = false;
     570}
     571
     572//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     573
     574inline void G4VEmProcess::SetStartFromNullFlag(G4bool val)
     575{
     576  startFromNull = val;
     577}
     578
     579//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     580
     581inline void G4VEmProcess::InitialiseStep(const G4Track& track)
     582{
     583  preStepKinEnergy = track.GetKineticEnergy();
     584  DefineMaterial(track.GetMaterialCutsCouple());
     585  if (theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX;
     586}
     587
    350588//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    351589
    … …  
    358596    mfpKinEnergy = DBL_MAX;
    359597  }
    360 }
    361 
    362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    363 
    364 inline void G4VEmProcess::InitialiseStep(const G4Track& track)
    365 {
    366   preStepKinEnergy = track.GetKineticEnergy();
    367   DefineMaterial(track.GetMaterialCutsCouple());
    368   if (theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX;
    369 }
    370 
    371 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    372 
    373 inline G4double G4VEmProcess::GetLambda(G4double& kineticEnergy,
    374                                         const G4MaterialCutsCouple* couple)
    375 {
    376   DefineMaterial(couple);
    377   return GetCurrentLambda(kineticEnergy);
    378 }
    379 
    380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    381 
    382 inline G4double G4VEmProcess::GetCurrentLambda(G4double e)
    383 {
    384   G4double x = 0.0;
    385   if(theLambdaTable) x = GetLambdaFromTable(e);
    386   else               x = ComputeCurrentLambda(e);
    387   return x;
    388 }
    389 
    390 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    391 
    392 inline G4double G4VEmProcess::RecalculateLambda(G4double e,
    393                                                 const G4MaterialCutsCouple* couple)
    394 {
    395   DefineMaterial(couple);
    396   return ComputeCurrentLambda(e);
    397 }
    398 
    399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    400 
    401 inline G4double G4VEmProcess::ComputeCurrentLambda(G4double e)
    402 {
    403   G4VEmModel* currentModel = SelectModel(e);
    404   G4double x = 0.0;
    405   if(currentModel)
    406     x = currentModel->CrossSectionPerVolume(currentMaterial,particle,
    407                                             e,(*theCuts)[currentMaterialIndex]);
    408   return x;
    409 }
    410 
    411 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    412 
    413 inline G4double G4VEmProcess::GetLambdaFromTable(G4double e)
    414 {
    415   G4bool b;
    416   return (((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b));
    417598}
    418599
    … …  
    442623//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    443624
    444 inline G4double G4VEmProcess::MeanFreePath(const G4Track& track)
    445 {
    446   DefineMaterial(track.GetMaterialCutsCouple());
    447   preStepLambda = GetCurrentLambda(track.GetKineticEnergy());
    448   G4double x = DBL_MAX;
    449   if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
     625inline G4double G4VEmProcess::GetLambdaFromTable(G4double e)
     626{
     627  G4bool b;
     628  return (((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b));
     629}
     630
     631//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     632
     633inline G4double G4VEmProcess::GetCurrentLambda(G4double e)
     634{
     635  G4double x = 0.0;
     636  if(theLambdaTable) x = GetLambdaFromTable(e);
     637  else               x = ComputeCurrentLambda(e);
    450638  return x;
    451639}
    … …  
    453641//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    454642
    455 inline G4VEmModel* G4VEmProcess::SelectModel(G4double& kinEnergy)
    456 {
    457   return modelManager->SelectModel(kinEnergy, currentMaterialIndex);
    458 }
    459 
    460 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    461 
    462 inline G4VEmModel* G4VEmProcess::SelectModelForMaterial(
    463                                    G4double kinEnergy, size_t& idxRegion) const
    464 {
    465   return modelManager->SelectModel(kinEnergy, idxRegion);
    466 }
    467 
    468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    469 
    470 inline const G4ParticleDefinition* G4VEmProcess::Particle() const
    471 {
    472   return particle;
    473 }
    474 
    475 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    476 
    477 inline const G4ParticleDefinition* G4VEmProcess::SecondaryParticle() const
    478 {
    479   return secondaryParticle;
    480 }
    481 
    482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    483 
    484 inline G4double G4VEmProcess::GetGammaEnergyCut()
    485 {
    486   return (*theCutsGamma)[currentMaterialIndex];
    487 }
    488 
    489 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    490 
    491 inline G4double G4VEmProcess::GetElectronEnergyCut()
    492 {
    493   return (*theCutsElectron)[currentMaterialIndex];
    494 }
    495 
    496 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    497 
    498 inline void G4VEmProcess::SetLambdaFactor(G4double val)
    499 {
    500   if(val > 0.0 && val <= 1.0) lambdaFactor = val;
    501 }
    502 
    503 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    504 
    505 inline G4VEmModel* G4VEmProcess::GetModelByIndex(G4int idx, G4bool ver)
    506 {
    507   return modelManager->GetModel(idx, ver);
    508 }
    509 
    510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    511 
    512 inline G4ParticleChangeForGamma* G4VEmProcess::GetParticleChange()
    513 {
    514   return &fParticleChange;
    515 }
    516 
    517 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    518 
    519 inline void G4VEmProcess::SetParticle(const G4ParticleDefinition* p)
    520 {
    521   particle = p;
    522 }
    523 
    524 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    525 
    526 inline void G4VEmProcess::SetSecondaryParticle(const G4ParticleDefinition* p)
    527 {
    528   secondaryParticle = p;
    529 }
    530 
    531 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    532 
    533 inline void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* p,
    534                                      const G4Region* region)
    535 {
    536   G4VEmFluctuationModel* fm = 0;
    537   modelManager->AddEmModel(order, p, fm, region);
    538   if(p) p->SetParticleChange(pParticleChange);
    539 }
    540 
    541 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    542 
    543 inline void G4VEmProcess::SetModel(G4VEmModel* model)
    544 {
    545   selectedModel = model;
    546 }
    547 
    548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    549 
    550 inline G4VEmModel* G4VEmProcess::Model()
    551 {
    552   return selectedModel;
    553 }
    554 
    555 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    556 
    557 inline void G4VEmProcess::UpdateEmModel(const G4String& nam,
    558                                         G4double emin, G4double emax)
    559 {
    560   modelManager->UpdateEmModel(nam, emin, emax);
    561 }
    562 
    563 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    564 
    565 inline G4double G4VEmProcess::ComputeCrossSectionPerAtom(
    566                  G4double kineticEnergy, G4double Z, G4double A, G4double cut)
    567 {
    568   G4VEmModel* model = SelectModel(kineticEnergy);
    569   return model->ComputeCrossSectionPerAtom(particle,kineticEnergy,Z,A,cut);
    570 }
    571 
    572 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    573 
    574 inline void G4VEmProcess::SetLambdaBinning(G4int nbins)
    575 {
    576   nLambdaBins = nbins;
    577 }
    578 
    579 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    580 
    581 inline G4int G4VEmProcess::LambdaBinning() const
    582 {
    583   return nLambdaBins;
    584 }
    585 
    586 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    587 
    588 inline void G4VEmProcess::SetMinKinEnergy(G4double e)
    589 {
    590   minKinEnergy = e;
    591 }
    592 
    593 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    594 
    595 inline G4double G4VEmProcess::MinKinEnergy() const
    596 {
    597   return minKinEnergy;
    598 }
    599 
    600 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    601 
    602 inline void G4VEmProcess::SetMaxKinEnergy(G4double e)
    603 {
    604   maxKinEnergy = e;
    605 }
    606 
    607 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    608 
    609 inline G4double G4VEmProcess::MaxKinEnergy() const
    610 {
    611   return maxKinEnergy;
    612 }
    613 
    614 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    615 
    616 inline void G4VEmProcess::SetPolarAngleLimit(G4double val)
    617 {
    618   if(val < 0.0)     polarAngleLimit = 0.0;
    619   else if(val > pi) polarAngleLimit = pi;
    620   else              polarAngleLimit = val;
    621 }
    622 
    623 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    624 
    625 inline G4double G4VEmProcess::PolarAngleLimit() const
    626 {
    627   return polarAngleLimit;
    628 }
    629 
    630 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    631 
    632 inline void G4VEmProcess::ActivateDeexcitation(G4bool, const G4Region*)
    633 {}
    634 
    635 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    636 
    637 inline const G4PhysicsTable* G4VEmProcess::LambdaTable() const
    638 {
    639   return theLambdaTable;
    640 }
    641 
    642 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    643 
    644 inline void G4VEmProcess::SetIntegral(G4bool val)
    645 {
    646   if(particle && particle != theGamma) integral = val;
    647   if(integral) buildLambdaTable = true;
    648 }
    649 
    650 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    651 
    652 inline G4bool G4VEmProcess::IsIntegral() const
    653 {
    654   return integral;
    655 }
    656 
    657 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    658 
    659 inline void G4VEmProcess::SetBuildTableFlag(G4bool val)
    660 {
    661   buildLambdaTable = val;
    662   if(!val) integral = false;
    663 }
    664 
    665 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    666 
    667 inline void G4VEmProcess::SetStartFromNullFlag(G4bool val)
    668 {
    669   startFromNull = val;
    670 }
    671 
    672 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    673 
    674 inline void G4VEmProcess::SetApplyCuts(G4bool val)
    675 {
    676   applyCuts = val;
    677 }
    678 
    679 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    680 
    681 inline size_t G4VEmProcess::CurrentMaterialCutsCoupleIndex() const
    682 {
    683   return currentMaterialIndex;
     643inline G4double G4VEmProcess::ComputeCurrentLambda(G4double e)
     644{
     645  SelectModel(e);
     646  G4double x = 0.0;
     647  if(currentModel) {
     648    x = currentModel->CrossSectionPerVolume(currentMaterial,particle,
     649                                            e,(*theCuts)[currentMaterialIndex]);
     650  }
     651  return x;
    684652}
    685653

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEnergyLossProcess.hh,v 1.83 2008/09/12 16:19:01 vnivanch Exp $
     26// $Id: G4VEnergyLossProcess.hh,v 1.87 2009/04/07 18:39:47 vnivanch Exp $
    2727// GEANT4 tag $Name:
    2828//
    … …  
    126126  virtual ~G4VEnergyLossProcess();
    127127
     128private:
     129  // clean vectors and arrays
     130  void Clean();
     131
    128132  //------------------------------------------------------------------------
    129133  // Virtual methods to be implemented in concrete processes
    130134  //------------------------------------------------------------------------
    131135
     136public:
    132137  virtual G4bool IsApplicable(const G4ParticleDefinition& p) = 0;
    133138 
    … …  
    143148  //------------------------------------------------------------------------
    144149
    145 protected:
    146 
    147150  virtual G4double MinPrimaryEnergy(const G4ParticleDefinition*,
    148151                                    const G4Material*, G4double cut);
    149152
    150153  //------------------------------------------------------------------------
    151   // Virtual methods common to all EM ContinuousDiscrete processes
    152   // Further inheritance is not assumed
     154  // Virtual methods implementation common to all EM ContinuousDiscrete
     155  // processes. Further inheritance is not assumed
    153156  //------------------------------------------------------------------------
    154157
    155158public:
    156159
     160  // prepare all tables
     161  void PreparePhysicsTable(const G4ParticleDefinition&);
     162
     163  // build all tables
     164  void BuildPhysicsTable(const G4ParticleDefinition&);
     165
     166  // build a table
     167  G4PhysicsTable* BuildDEDXTable(G4EmTableType tType = fRestricted);
     168
     169  // build a table
     170  G4PhysicsTable* BuildLambdaTable(G4EmTableType tType = fRestricted);
     171
     172  // summary printout after initialisation
    157173  void PrintInfoDefinition();
    158174
    159   void PreparePhysicsTable(const G4ParticleDefinition&);
    160 
    161   void BuildPhysicsTable(const G4ParticleDefinition&);
    162 
     175  // Add subcutoff option for the region
     176  void ActivateSubCutoff(G4bool val, const G4Region* region = 0);
     177
     178  // Activate deexcitation code for region
     179  void ActivateDeexcitation(G4bool, const G4Region* region = 0);
     180
     181  // Step limit from AlongStep
    163182  G4double AlongStepGetPhysicalInteractionLength(const G4Track&,
    164183                                                 G4double  previousStepSize,
    … …  
    167186                                                 G4GPILSelection* selection);
    168187
     188  // Step limit from cross section
    169189  G4double PostStepGetPhysicalInteractionLength(const G4Track& track,
    170190                                                G4double   previousStepSize,
    171191                                                G4ForceCondition* condition);
    172192
     193  // AlongStep computations
    173194  G4VParticleChange* AlongStepDoIt(const G4Track&, const G4Step&);
    174195
     196  // Sampling of secondaries in vicinity of geometrical boundary
     197  void SampleSubCutSecondaries(std::vector<G4Track*>&, const G4Step&,
     198                               G4VEmModel* model, G4int matIdx,
     199                               G4double& extraEdep);
     200
     201  // PostStep sampling of secondaries
    175202  G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&);
    176203
    177   // Store PhysicsTable in a file.
    178   // Return false in case of failure at I/O
     204  // Store all PhysicsTable in files.
     205  // Return false in case of any fatal failure at I/O 
    179206  G4bool StorePhysicsTable(const G4ParticleDefinition*,
    180207                           const G4String& directory,
    181208                           G4bool ascii = false);
    182209
    183   // Retrieve Physics from a file.
    184   // (return true if the Physics Table can be build by using file)
    185   // (return false if the process has no functionality or in case of failure)
    186   // File name should is constructed as processName+particleName and the
    187   // should be placed under the directory specifed by the argument.
     210  // Retrieve all Physics from a files.
     211  // Return true if all the Physics Table are built.
     212  // Return false if any fatal failure.
    188213  G4bool RetrievePhysicsTable(const G4ParticleDefinition*,
    189214                              const G4String& directory,
    190215                              G4bool ascii);
    191216
     217private:
     218  // store a table
     219  G4bool StoreTable(const G4ParticleDefinition* p,
     220                    G4PhysicsTable*, G4bool ascii,
     221                    const G4String& directory,
     222                    const G4String& tname);
     223
     224  // retrieve a table
     225  G4bool RetrieveTable(const G4ParticleDefinition* p,
     226                       G4PhysicsTable*, G4bool ascii,
     227                       const G4String& directory,
     228                       const G4String& tname,
     229                       G4bool mandatory);
     230
     231  //------------------------------------------------------------------------
     232  // Public interface to cross section, mfp and sampling of fluctuations
     233  // These methods are not used in run time
     234  //------------------------------------------------------------------------
     235
     236public:
     237  // access to dispersion of restricted energy loss
     238  G4double GetDEDXDispersion(const G4MaterialCutsCouple *couple,
     239                             const G4DynamicParticle* dp,
     240                             G4double length);
     241
     242  // Access to cross section table
     243  G4double CrossSectionPerVolume(G4double kineticEnergy,
     244                                 const G4MaterialCutsCouple* couple);
     245
     246  // access to cross section
     247  G4double MeanFreePath(const G4Track& track);
     248
     249  // access to step limit
     250  G4double ContinuousStepLimit(const G4Track& track,
     251                               G4double previousStepSize,
     252                               G4double currentMinimumStep,
     253                               G4double& currentSafety);
     254
    192255protected:
    193256
     257  // implementation of the pure virtual method
    194258  G4double GetMeanFreePath(const G4Track& track,
    195259                           G4double previousStepSize,
    196260                           G4ForceCondition* condition);
    197261
     262  // implementation of the pure virtual method
    198263  G4double GetContinuousStepLimit(const G4Track& track,
    199264                                  G4double previousStepSize,
    … …  
    202267
    203268  //------------------------------------------------------------------------
    204   // Specific methods for along/post step EM processes
    205   //------------------------------------------------------------------------
     269  // Run time method which may be also used by derived processes
     270  //------------------------------------------------------------------------
     271
     272  // creeation of an empty vector for cross section
     273  G4PhysicsVector* LambdaPhysicsVector(const G4MaterialCutsCouple*,
     274                                       G4double cut);
     275
     276  inline size_t CurrentMaterialCutsCoupleIndex() const;
     277
     278  inline G4double GetCurrentRange() const;
     279
     280  //------------------------------------------------------------------------
     281  // Specific methods to set, access, modify models
     282  //------------------------------------------------------------------------
     283
     284  // Select model in run time
     285  inline void SelectModel(G4double kinEnergy);
    206286
    207287public:
    208 
     288  // Select model by energy and region index
     289  inline G4VEmModel* SelectModelForMaterial(G4double kinEnergy,
     290                                            size_t& idx) const;
     291
     292  // Add EM model coupled with fluctuation model for region, smaller value
     293  // of order defines which pair of models will be selected for a given
     294  // energy interval 
     295  void AddEmModel(G4int, G4VEmModel*,
     296                  G4VEmFluctuationModel* fluc = 0,
     297                  const G4Region* region = 0);
     298
     299  // Define new energy range for the model identified by the name
     300  void UpdateEmModel(const G4String&, G4double, G4double);
     301
     302  // Assign a model to a process
     303  void SetEmModel(G4VEmModel*, G4int index=1);
     304 
     305  // return the assigned model
     306  G4VEmModel* EmModel(G4int index=1);
     307 
     308  // Access to models
     309  G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false);
     310
     311  G4int NumberOfModels();
     312
     313  // Assign a fluctuation model to a process
     314  void SetFluctModel(G4VEmFluctuationModel*);
     315 
     316  // return the assigned fluctuation model
     317  inline G4VEmFluctuationModel* FluctModel();
     318   
     319  //------------------------------------------------------------------------
     320  // Define and access particle type
     321  //------------------------------------------------------------------------
     322
     323protected:
     324  inline void SetParticle(const G4ParticleDefinition* p);
     325  inline void SetSecondaryParticle(const G4ParticleDefinition* p);
     326
     327public:
     328  inline void SetBaseParticle(const G4ParticleDefinition* p);
     329  inline const G4ParticleDefinition* Particle() const;
     330  inline const G4ParticleDefinition* BaseParticle() const;
     331  inline const G4ParticleDefinition* SecondaryParticle() const;
     332
     333  //------------------------------------------------------------------------
     334  // Get/set parameters to configure the process at initialisation time
     335  //------------------------------------------------------------------------
     336
     337  // Add subcutoff process (bremsstrahlung) to sample secondary
     338  // particle production in vicinity of the geometry boundary
    209339  void AddCollaborativeProcess(G4VEnergyLossProcess*);
    210340
    211   void SampleSubCutSecondaries(std::vector<G4Track*>&, const G4Step&,
    212                                G4VEmModel* model, G4int matIdx,
    213                                G4double& extraEdep);
    214 
    215   G4double GetDEDXDispersion(const G4MaterialCutsCouple *couple,
    216                              const G4DynamicParticle* dp,
    217                              G4double length);
    218 
    219   //------------------------------------------------------------------------
    220   // Specific methods to build and access Physics Tables
    221   //------------------------------------------------------------------------
    222 
    223   G4PhysicsTable* BuildDEDXTable(G4EmTableType tType = fRestricted);
    224 
    225   G4PhysicsTable* BuildLambdaTable(G4EmTableType tType = fRestricted);
     341  inline void SetLossFluctuations(G4bool val);
     342  inline void SetRandomStep(G4bool val);
     343
     344  inline void SetIntegral(G4bool val);
     345  inline G4bool IsIntegral() const;
     346
     347  // Set/Get flag "isIonisation"
     348  inline void SetIonisation(G4bool val);
     349  inline G4bool IsIonisationProcess() const;
     350
     351  // Redefine parameteters for stepping control
     352  //
     353  inline void SetLinearLossLimit(G4double val);
     354  inline void SetMinSubRange(G4double val);
     355  inline void SetLambdaFactor(G4double val);
     356  inline void SetStepFunction(G4double v1, G4double v2);
     357
     358  inline G4int NumberOfSubCutoffRegions() const;
     359  inline G4int NumberOfDERegions() const;
     360
     361  //------------------------------------------------------------------------
     362  // Specific methods to path Physics Tables to the process
     363  //------------------------------------------------------------------------
    226364
    227365  void SetDEDXTable(G4PhysicsTable* p, G4EmTableType tType);
    228366  void SetCSDARangeTable(G4PhysicsTable* pRange);
    229367  void SetRangeTableForLoss(G4PhysicsTable* p);
     368  void SetSecondaryRangeTable(G4PhysicsTable* p);
    230369  void SetInverseRangeTable(G4PhysicsTable* p);
    231   void SetSecondaryRangeTable(G4PhysicsTable* p);
    232370
    233371  void SetLambdaTable(G4PhysicsTable* p);
    … …  
    251389  // Max kinetic energy for tables
    252390  inline void SetMaxKinEnergyForCSDARange(G4double e);
     391
     392  // Return values for given G4MaterialCutsCouple
     393  inline G4double GetDEDX(G4double& kineticEnergy, const G4MaterialCutsCouple*);
     394  inline G4double GetDEDXForSubsec(G4double& kineticEnergy,
     395                                   const G4MaterialCutsCouple*);
     396  inline G4double GetRange(G4double& kineticEnergy, const G4MaterialCutsCouple*);
     397  inline G4double GetCSDARange(G4double& kineticEnergy, const G4MaterialCutsCouple*);
     398  inline G4double GetRangeForLoss(G4double& kineticEnergy, const G4MaterialCutsCouple*);
     399  inline G4double GetKineticEnergy(G4double& range, const G4MaterialCutsCouple*);
     400  inline G4double GetLambda(G4double& kineticEnergy, const G4MaterialCutsCouple*);
     401
     402  inline G4bool TablesAreBuilt() const;
    253403
    254404  // Access to specific tables
    … …  
    264414  inline G4PhysicsTable* SubLambdaTable();
    265415
    266   // Return values for given G4MaterialCutsCouple
    267   inline G4double GetDEDX(G4double& kineticEnergy, const G4MaterialCutsCouple*);
    268   inline G4double GetDEDXForSubsec(G4double& kineticEnergy,
    269                                    const G4MaterialCutsCouple*);
    270   inline G4double GetRange(G4double& kineticEnergy, const G4MaterialCutsCouple*);
    271   inline G4double GetCSDARange(G4double& kineticEnergy, const G4MaterialCutsCouple*);
    272   inline G4double GetRangeForLoss(G4double& kineticEnergy, const G4MaterialCutsCouple*);
    273   inline G4double GetKineticEnergy(G4double& range, const G4MaterialCutsCouple*);
    274   inline G4double GetLambda(G4double& kineticEnergy, const G4MaterialCutsCouple*);
    275 
    276   inline G4bool TablesAreBuilt() const;
    277 
    278   //------------------------------------------------------------------------
    279   // Define and access particle type
    280   //------------------------------------------------------------------------
    281 
    282   inline void SetBaseParticle(const G4ParticleDefinition* p);
    283   inline const G4ParticleDefinition* Particle() const;
    284   inline const G4ParticleDefinition* BaseParticle() const;
    285   inline const G4ParticleDefinition* SecondaryParticle() const;
    286 
    287   //------------------------------------------------------------------------
    288   // Specific methods to set, access, modify models
    289   //------------------------------------------------------------------------
    290 
    291   // Add EM model coupled with fluctuation model for the region
    292   inline void AddEmModel(G4int, G4VEmModel*,
    293                          G4VEmFluctuationModel* fluc = 0,
    294                          const G4Region* region = 0);
    295 
    296   // Assign a model to a process
    297   inline void SetEmModel(G4VEmModel*, G4int index=1);
    298  
    299   // return the assigned model
    300   inline G4VEmModel* EmModel(G4int index=1);
    301  
    302   // Assign a fluctuation model to a process
    303   inline void SetFluctModel(G4VEmFluctuationModel*);
    304  
    305   // return the assigned fluctuation model
    306   inline G4VEmFluctuationModel* FluctModel();
    307    
    308   // Define new energy range for the model identified by the name
    309   inline void UpdateEmModel(const G4String&, G4double, G4double);
    310 
    311   // Access to models
    312   inline G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false);
    313 
    314   inline G4int NumberOfModels();
    315 
    316   //------------------------------------------------------------------------
    317   // Get/set parameters used for simulation of energy loss
    318   //------------------------------------------------------------------------
    319 
    320   inline void SetLossFluctuations(G4bool val);
    321   inline void SetRandomStep(G4bool val);
    322   inline void SetIntegral(G4bool val);
    323   inline G4bool IsIntegral() const;
    324 
    325   // Set/Get flag "isIonisation"
    326   inline void SetIonisation(G4bool val);
    327   inline G4bool IsIonisationProcess() const;
    328 
    329   // Redefine parameteters for stepping control
    330   //
    331   inline void SetLinearLossLimit(G4double val);
    332   inline void SetMinSubRange(G4double val);
    333   inline void SetStepFunction(G4double v1, G4double v2);
    334   inline void SetLambdaFactor(G4double val);
    335 
    336 
    337   // Add subcutoff option for the region
    338   void ActivateSubCutoff(G4bool val, const G4Region* region = 0);
    339 
    340   inline G4int NumberOfSubCutoffRegions() const;
    341 
    342   // Activate deexcitation code
    343   virtual void ActivateDeexcitation(G4bool, const G4Region* region = 0);
    344 
    345   //------------------------------------------------------------------------
    346   // Public interface to helper functions
    347   //------------------------------------------------------------------------
    348 
    349   inline
    350   G4VEmModel* SelectModelForMaterial(G4double kinEnergy, size_t& idx) const;
    351 
    352   inline G4double MeanFreePath(const G4Track& track);
    353 
    354   inline G4double ContinuousStepLimit(const G4Track& track,
    355                                       G4double previousStepSize,
    356                                       G4double currentMinimumStep,
    357                                       G4double& currentSafety);
    358 
    359416  //------------------------------------------------------------------------
    360417  // Run time method for simulation of ionisation
    … …  
    367424  inline void SetDynamicMassCharge(G4double massratio, G4double charge2ratio);
    368425
    369   // Access to cross section table
    370   G4double CrossSectionPerVolume(G4double kineticEnergy,
    371                                  const G4MaterialCutsCouple* couple);
    372 
    373 protected:
    374 
    375   G4PhysicsVector* LambdaPhysicsVector(const G4MaterialCutsCouple*,
    376                                        G4double cut);
    377 
    378   inline G4ParticleChangeForLoss* GetParticleChange();
    379 
    380   inline void SetParticle(const G4ParticleDefinition* p);
    381 
    382   inline void SetSecondaryParticle(const G4ParticleDefinition* p);
    383 
    384   inline void SelectModel(G4double kinEnergy);
    385 
    386   inline size_t CurrentMaterialCutsCoupleIndex() const;
    387 
    388   inline G4double GetCurrentRange() const;
    389 
    390426private:
    391 
    392   //------------------------------------------------------------------------
    393   // Management of tables
    394   //------------------------------------------------------------------------
    395 
    396   void Clear();
    397 
    398   G4bool StoreTable(const G4ParticleDefinition* p,
    399                     G4PhysicsTable*, G4bool ascii,
    400                     const G4String& directory,
    401                     const G4String& tname);
    402 
    403   G4bool RetrieveTable(const G4ParticleDefinition* p,
    404                        G4PhysicsTable*, G4bool ascii,
    405                        const G4String& directory,
    406                        const G4String& tname,
    407                        G4bool mandatory);
    408427
    409428  // define material and indexes
    410429  inline void DefineMaterial(const G4MaterialCutsCouple* couple);
    411430
    412   // Returnd values for scaled energy using mass of the base particle
    413   //
     431  //------------------------------------------------------------------------
     432  // Compute values using scaling relation, mass and charge of based particle
     433  //------------------------------------------------------------------------
     434
    414435  inline G4double GetDEDXForScaledEnergy(G4double scaledKinEnergy);
    415436  inline G4double GetSubDEDXForScaledEnergy(G4double scaledKinEnergy);
    … …  
    418439  inline G4double GetScaledRangeForScaledEnergy(G4double scaledKinEnergy);
    419440  inline G4double GetLimitScaledRangeForScaledEnergy(G4double scaledKinEnergy);
     441  inline G4double ScaledKinEnergyForLoss(G4double range);
    420442  inline G4double GetLambdaForScaledEnergy(G4double scaledKinEnergy);
    421   inline G4double ScaledKinEnergyForLoss(G4double range);
    422443  inline void ComputeLambdaForScaledEnergy(G4double scaledKinEnergy);
    423444
    424445  // hide  assignment operator
    425 
    426446  G4VEnergyLossProcess(G4VEnergyLossProcess &);
    427447  G4VEnergyLossProcess & operator=(const G4VEnergyLossProcess &right);
    … …  
    444464  G4VEmFluctuationModel*                fluctModel;
    445465  std::vector<const G4Region*>          scoffRegions;
     466  std::vector<const G4Region*>          deRegions;
    446467  G4int                                 nSCoffRegions;
    447   G4int*                                idxSCoffRegions;
     468  G4int                                 nDERegions;
     469  G4bool*                               idxSCoffRegions;
     470  G4bool*                               idxDERegions;
    448471
    449472  std::vector<G4VEnergyLossProcess*>    scProcesses;
    … …  
    493516  G4bool   isIonisation;
    494517  G4bool   useSubCutoff;
     518  G4bool   useDeexcitation;
    495519
    496520protected:
    … …  
    531555//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    532556
    533 inline void G4VEnergyLossProcess::DefineMaterial(
    534             const G4MaterialCutsCouple* couple)
    535 {
    536   if(couple != currentCouple) {
    537     currentCouple   = couple;
    538     currentMaterial = couple->GetMaterial();
    539     currentMaterialIndex = couple->GetIndex();
    540     mfpKinEnergy = DBL_MAX;
    541   }
     557inline size_t G4VEnergyLossProcess::CurrentMaterialCutsCoupleIndex() const
     558{
     559  return currentMaterialIndex;
     560}
     561
     562//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     563 
     564inline G4double G4VEnergyLossProcess::GetCurrentRange() const
     565{
     566  return fRange;
     567}
     568
     569//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     570
     571inline void G4VEnergyLossProcess::SelectModel(G4double kinEnergy)
     572{
     573  currentModel = modelManager->SelectModel(kinEnergy, currentMaterialIndex);
     574  currentModel->SetCurrentCouple(currentCouple);
     575}
     576
     577//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     578
     579inline G4VEmModel* G4VEnergyLossProcess::SelectModelForMaterial(
     580                   G4double kinEnergy, size_t& idx) const
     581{
     582  return modelManager->SelectModel(kinEnergy, idx);
     583}
     584
     585//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     586
     587inline void G4VEnergyLossProcess::SetFluctModel(G4VEmFluctuationModel* p)
     588{
     589  fluctModel = p;
     590}
     591
     592//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     593
     594inline G4VEmFluctuationModel* G4VEnergyLossProcess::FluctModel()
     595{
     596  return fluctModel;
     597}
     598
     599//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     600
     601inline void G4VEnergyLossProcess::SetParticle(const G4ParticleDefinition* p)
     602{
     603  particle = p;
     604}
     605
     606//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     607
     608inline void G4VEnergyLossProcess::SetSecondaryParticle(const G4ParticleDefinition* p)
     609{
     610  secondaryParticle = p;
     611}
     612
     613//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     614
     615inline void G4VEnergyLossProcess::SetBaseParticle(const G4ParticleDefinition* p)
     616{
     617  baseParticle = p;
     618}
     619
     620//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     621
     622inline const G4ParticleDefinition* G4VEnergyLossProcess::Particle() const
     623{
     624  return particle;
     625}
     626
     627//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     628
     629inline const G4ParticleDefinition* G4VEnergyLossProcess::BaseParticle() const
     630{
     631  return baseParticle;
     632}
     633
     634//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     635
     636inline const G4ParticleDefinition* G4VEnergyLossProcess::SecondaryParticle() const
     637{
     638  return secondaryParticle;
     639}
     640
     641//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     642
     643inline void G4VEnergyLossProcess::SetLossFluctuations(G4bool val)
     644{
     645  lossFluctuationFlag = val;
     646}
     647
     648//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     649
     650inline void G4VEnergyLossProcess::SetRandomStep(G4bool val)
     651{
     652  rndmStepFlag = val;
     653}
     654
     655//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     656
     657inline void G4VEnergyLossProcess::SetIntegral(G4bool val)
     658{
     659  integral = val;
     660}
     661
     662//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     663 
     664inline G4bool G4VEnergyLossProcess::IsIntegral() const
     665{
     666  return integral;
     667}
     668
     669//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     670
     671inline void G4VEnergyLossProcess::SetIonisation(G4bool val)
     672{
     673  isIonisation = val;
     674  if(val) aGPILSelection = CandidateForSelection;
     675  else    aGPILSelection = NotCandidateForSelection;
     676}
     677
     678//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     679
     680inline G4bool G4VEnergyLossProcess::IsIonisationProcess() const
     681{
     682  return isIonisation;
     683}
     684
     685//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     686
     687inline void G4VEnergyLossProcess::SetLinearLossLimit(G4double val)
     688{
     689  linLossLimit = val;
     690}
     691
     692//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     693
     694inline void G4VEnergyLossProcess::SetMinSubRange(G4double val)
     695{
     696  minSubRange = val;
     697}
     698
     699//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     700
     701inline void G4VEnergyLossProcess::SetLambdaFactor(G4double val)
     702{
     703  if(val > 0.0 && val <= 1.0) lambdaFactor = val;
     704}
     705
     706//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     707
     708void G4VEnergyLossProcess::SetStepFunction(G4double v1, G4double v2)
     709{
     710  dRoverRange = v1;
     711  finalRange = v2;
     712  if (dRoverRange > 0.999) dRoverRange = 1.0;
     713  currentCouple = 0;
     714  mfpKinEnergy  = DBL_MAX;
     715}
     716
     717//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     718
     719inline G4int G4VEnergyLossProcess::NumberOfSubCutoffRegions() const
     720{
     721  return nSCoffRegions;
     722}
     723
     724//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     725
     726inline G4int G4VEnergyLossProcess::NumberOfDERegions() const
     727{
     728  return nDERegions;
     729}
     730
     731//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     732
     733inline void G4VEnergyLossProcess::SetDEDXBinning(G4int nbins)
     734{
     735  nBins = nbins;
     736}
     737
     738//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     739
     740inline void G4VEnergyLossProcess::SetLambdaBinning(G4int nbins)
     741{
     742  nBins = nbins;
     743}
     744
     745//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     746
     747inline void G4VEnergyLossProcess::SetDEDXBinningForCSDARange(G4int nbins)
     748{
     749  nBinsCSDA = nbins;
     750}
     751
     752//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     753
     754inline void G4VEnergyLossProcess::SetMinKinEnergy(G4double e)
     755{
     756  minKinEnergy = e;
     757}
     758
     759//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     760
     761inline G4double G4VEnergyLossProcess::MinKinEnergy() const
     762{
     763  return minKinEnergy;
     764}
     765
     766//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     767
     768inline void G4VEnergyLossProcess::SetMaxKinEnergy(G4double e)
     769{
     770  maxKinEnergy = e;
     771  if(e < maxKinEnergyCSDA) maxKinEnergyCSDA = e;
     772}
     773
     774//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     775
     776inline G4double G4VEnergyLossProcess::MaxKinEnergy() const
     777{
     778  return maxKinEnergy;
     779}
     780
     781//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     782
     783inline void G4VEnergyLossProcess::SetMaxKinEnergyForCSDARange(G4double e)
     784{
     785  maxKinEnergyCSDA = e;
    542786}
    543787
    … …  
    558802  DefineMaterial(couple);
    559803  return GetSubDEDXForScaledEnergy(kineticEnergy*massRatio);
    560 }
    561 
    562 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    563 
    564 inline G4double G4VEnergyLossProcess::GetDEDXForScaledEnergy(G4double e)
    565 {
    566   G4bool b;
    567   G4double x =
    568     ((*theDEDXTable)[currentMaterialIndex]->GetValue(e, b))*chargeSqRatio;
    569   if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
    570   return x;
    571 }
    572 
    573 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    574 
    575 inline G4double G4VEnergyLossProcess::GetSubDEDXForScaledEnergy(G4double e)
    576 {
    577   G4bool b;
    578   G4double x =
    579     ((*theDEDXSubTable)[currentMaterialIndex]->GetValue(e, b))*chargeSqRatio;
    580   if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
    581   return x;
    582 }
    583 
    584 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    585 
    586 inline G4double G4VEnergyLossProcess::GetIonisationForScaledEnergy(G4double e)
    587 {
    588   G4bool b;
    589   G4double x = 0.0;
    590   //  if(theIonisationTable) {
    591   x = ((*theIonisationTable)[currentMaterialIndex]->GetValue(e, b))
    592     *chargeSqRatio;
    593   if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
    594   //}
    595   return x;
    596 }
    597 
    598 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    599 
    600 inline
    601 G4double G4VEnergyLossProcess::GetSubIonisationForScaledEnergy(G4double e)
    602 {
    603   G4bool b;
    604   G4double x = 0.0;
    605   //if(theIonisationSubTable) {
    606   x = ((*theIonisationSubTable)[currentMaterialIndex]->GetValue(e, b))
    607     *chargeSqRatio;
    608   if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
    609   //}
    610   return x;
    611804}
    612805
    … …  
    643836//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    644837
     838inline G4double G4VEnergyLossProcess::GetRangeForLoss(
     839                G4double& kineticEnergy,
     840                const G4MaterialCutsCouple* couple)
     841{
     842  DefineMaterial(couple);
     843  G4double x = DBL_MAX;
     844  if(theRangeTableForLoss)
     845    x = GetScaledRangeForScaledEnergy(kineticEnergy*massRatio)*reduceFactor;
     846  //  G4cout << "Range from " << GetProcessName()
     847  //         << "  e= " << kineticEnergy << " r= " << x << G4endl;
     848  return x;
     849}
     850
     851//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     852
     853inline G4double G4VEnergyLossProcess::GetKineticEnergy(
     854                G4double& range,
     855                const G4MaterialCutsCouple* couple)
     856{
     857  DefineMaterial(couple);
     858  G4double r = range/reduceFactor;
     859  G4double e = ScaledKinEnergyForLoss(r)/massRatio;
     860  return e;
     861}
     862
     863//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     864
     865inline G4double G4VEnergyLossProcess::GetLambda(G4double& kineticEnergy,
     866                                          const G4MaterialCutsCouple* couple)
     867{
     868  DefineMaterial(couple);
     869  G4double x = 0.0;
     870  if(theLambdaTable) x = GetLambdaForScaledEnergy(kineticEnergy*massRatio);
     871  return x;
     872}
     873
     874//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     875
     876inline G4bool G4VEnergyLossProcess::TablesAreBuilt() const
     877{
     878  return  tablesAreBuilt;
     879}
     880
     881//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     882
     883inline G4PhysicsTable* G4VEnergyLossProcess::DEDXTable() const
     884{
     885  return theDEDXTable;
     886}
     887
     888//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     889
     890inline G4PhysicsTable* G4VEnergyLossProcess::DEDXTableForSubsec() const
     891{
     892  return theDEDXSubTable;
     893}
     894
     895//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     896
     897inline G4PhysicsTable* G4VEnergyLossProcess::DEDXunRestrictedTable() const
     898{
     899  return theDEDXunRestrictedTable;
     900}
     901
     902//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     903
     904inline G4PhysicsTable* G4VEnergyLossProcess::IonisationTable() const
     905{
     906  G4PhysicsTable* t = theDEDXTable;
     907  if(theIonisationTable) t = theIonisationTable;
     908  return t;
     909}
     910
     911//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     912
     913inline G4PhysicsTable* G4VEnergyLossProcess::IonisationTableForSubsec() const
     914{
     915  G4PhysicsTable* t = theDEDXSubTable;
     916  if(theIonisationSubTable) t = theIonisationSubTable;
     917  return t;
     918}
     919
     920//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     921
     922inline G4PhysicsTable* G4VEnergyLossProcess::CSDARangeTable() const
     923{
     924  return theCSDARangeTable;
     925}
     926
     927//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     928
     929inline G4PhysicsTable* G4VEnergyLossProcess::RangeTableForLoss() const
     930{
     931  return theRangeTableForLoss;
     932}
     933
     934//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     935
     936inline G4PhysicsTable* G4VEnergyLossProcess::InverseRangeTable() const
     937{
     938  return theInverseRangeTable;
     939}
     940
     941//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     942
     943inline G4PhysicsTable* G4VEnergyLossProcess::LambdaTable()
     944{
     945  return theLambdaTable;
     946}
     947
     948//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     949
     950inline G4PhysicsTable* G4VEnergyLossProcess::SubLambdaTable()
     951{
     952  return theSubLambdaTable;
     953}
     954
     955//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     956
     957inline G4double G4VEnergyLossProcess::SampleRange()
     958{
     959  G4double e = amu_c2*preStepKinEnergy/particle->GetPDGMass();
     960  G4bool b;
     961  G4double s = fRange*std::pow(10.,vstrag->GetValue(e,b));
     962  G4double x = fRange + G4RandGauss::shoot(0.0,s);
     963  if(x > 0.0) fRange = x;
     964  return fRange;
     965}
     966
     967//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     968
     969inline void G4VEnergyLossProcess::SetDynamicMassCharge(G4double massratio,
     970                                                       G4double charge2ratio)
     971{
     972  massRatio     = massratio;
     973  chargeSqRatio = charge2ratio;
     974  reduceFactor  = 1.0/(chargeSqRatio*massRatio);
     975}
     976
     977//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     978
     979inline void G4VEnergyLossProcess::DefineMaterial(
     980            const G4MaterialCutsCouple* couple)
     981{
     982  if(couple != currentCouple) {
     983    currentCouple   = couple;
     984    currentMaterial = couple->GetMaterial();
     985    currentMaterialIndex = couple->GetIndex();
     986    mfpKinEnergy = DBL_MAX;
     987  }
     988}
     989
     990//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     991
     992inline G4double G4VEnergyLossProcess::GetDEDXForScaledEnergy(G4double e)
     993{
     994  G4bool b;
     995  G4double x =
     996    ((*theDEDXTable)[currentMaterialIndex]->GetValue(e, b))*chargeSqRatio;
     997  if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
     998  return x;
     999}
     1000
     1001//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1002
     1003inline G4double G4VEnergyLossProcess::GetSubDEDXForScaledEnergy(G4double e)
     1004{
     1005  G4bool b;
     1006  G4double x =
     1007    ((*theDEDXSubTable)[currentMaterialIndex]->GetValue(e, b))*chargeSqRatio;
     1008  if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
     1009  return x;
     1010}
     1011
     1012//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1013
     1014inline G4double G4VEnergyLossProcess::GetIonisationForScaledEnergy(G4double e)
     1015{
     1016  G4bool b;
     1017  G4double x = 0.0;
     1018  //  if(theIonisationTable) {
     1019  x = ((*theIonisationTable)[currentMaterialIndex]->GetValue(e, b))
     1020    *chargeSqRatio;
     1021  if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
     1022  //}
     1023  return x;
     1024}
     1025
     1026//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1027
     1028inline
     1029G4double G4VEnergyLossProcess::GetSubIonisationForScaledEnergy(G4double e)
     1030{
     1031  G4bool b;
     1032  G4double x = 0.0;
     1033  //if(theIonisationSubTable) {
     1034  x = ((*theIonisationSubTable)[currentMaterialIndex]->GetValue(e, b))
     1035    *chargeSqRatio;
     1036  if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
     1037  //}
     1038  return x;
     1039}
     1040
     1041//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1042
     1043inline G4double G4VEnergyLossProcess::GetScaledRangeForScaledEnergy(G4double e)
     1044{
     1045  G4bool b;
     1046  G4double x = ((*theRangeTableForLoss)[currentMaterialIndex])->GetValue(e, b);
     1047  if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
     1048  return x;
     1049}
     1050
     1051//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1052
    6451053inline G4double G4VEnergyLossProcess::GetLimitScaledRangeForScaledEnergy(
    6461054                G4double e)
    … …  
    6571065  }
    6581066  return x;
    659 }
    660 
    661 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    662 
    663 inline G4double G4VEnergyLossProcess::GetRangeForLoss(
    664                 G4double& kineticEnergy,
    665                 const G4MaterialCutsCouple* couple)
    666 {
    667   DefineMaterial(couple);
    668   G4double x = DBL_MAX;
    669   if(theRangeTableForLoss)
    670     x = GetScaledRangeForScaledEnergy(kineticEnergy*massRatio)*reduceFactor;
    671   //  G4cout << "Range from " << GetProcessName()
    672   //         << "  e= " << kineticEnergy << " r= " << x << G4endl;
    673   return x;
    674 }
    675 
    676 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    677 
    678 inline G4double G4VEnergyLossProcess::GetScaledRangeForScaledEnergy(G4double e)
    679 {
    680   G4bool b;
    681   G4double x = ((*theRangeTableForLoss)[currentMaterialIndex])->GetValue(e, b);
    682   if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
    683   return x;
    684 }
    685 
    686 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    687 
    688 inline G4double G4VEnergyLossProcess::GetKineticEnergy(
    689                 G4double& range,
    690                 const G4MaterialCutsCouple* couple)
    691 {
    692   DefineMaterial(couple);
    693   G4double r = range/reduceFactor;
    694   G4double e = ScaledKinEnergyForLoss(r)/massRatio;
    695   return e;
    6961067}
    6971068
    … …  
    7111082  }
    7121083  return e;
    713 }
    714 
    715 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    716 
    717 inline G4double G4VEnergyLossProcess::GetLambda(G4double& kineticEnergy,
    718                                           const G4MaterialCutsCouple* couple)
    719 {
    720   DefineMaterial(couple);
    721   G4double x = 0.0;
    722   if(theLambdaTable) x = GetLambdaForScaledEnergy(kineticEnergy*massRatio);
    723   return x;
    7241084}
    7251085
    … …  
    7581118//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    7591119
    760 inline G4double G4VEnergyLossProcess::ContinuousStepLimit(
    761          const G4Track& track, G4double x, G4double y, G4double& z)
    762 {
    763   G4GPILSelection sel;
    764   return AlongStepGetPhysicalInteractionLength(track, x, y, z, &sel);
    765 }
    766 
    767 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    768 
    769 inline G4double G4VEnergyLossProcess::SampleRange()
    770 {
    771   G4double e = amu_c2*preStepKinEnergy/particle->GetPDGMass();
    772   G4bool b;
    773   G4double s = fRange*std::pow(10.,vstrag->GetValue(e,b));
    774   G4double x = fRange + G4RandGauss::shoot(0.0,s);
    775   if(x > 0.0) fRange = x;
    776   return fRange;
    777 }
    778 
    779 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    780 
    781 inline G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track)
    782 {
    783   DefineMaterial(track.GetMaterialCutsCouple());
    784   preStepLambda = GetLambdaForScaledEnergy(track.GetKineticEnergy()*massRatio);
    785   G4double x = DBL_MAX;
    786   if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
    787   return x;
    788 }
    789 
    790 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    791 
    792 inline G4double G4VEnergyLossProcess::MinPrimaryEnergy(
    793                 const G4ParticleDefinition*, const G4Material*, G4double cut)
    794 {
    795   return cut;
    796 }
    797 
    798 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    799 
    800 inline void G4VEnergyLossProcess::SelectModel(G4double kinEnergy)
    801 {
    802   currentModel = modelManager->SelectModel(kinEnergy, currentMaterialIndex);
    803 }
    804 
    805 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    806 
    807 inline G4VEmModel* G4VEnergyLossProcess::SelectModelForMaterial(
    808                    G4double kinEnergy, size_t& idx) const
    809 {
    810   return modelManager->SelectModel(kinEnergy, idx);
    811 }
    812 
    813 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    814 
    815 inline G4ParticleChangeForLoss* G4VEnergyLossProcess::GetParticleChange()
    816 {
    817   return &fParticleChange;
    818 }
    819 
    820 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    821 
    822 inline const G4ParticleDefinition* G4VEnergyLossProcess::Particle() const
    823 {
    824   return particle;
    825 }
    826 
    827 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    828 
    829 inline const G4ParticleDefinition* G4VEnergyLossProcess::BaseParticle() const
    830 {
    831   return baseParticle;
    832 }
    833 
    834 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    835 
    836 inline const G4ParticleDefinition* G4VEnergyLossProcess::SecondaryParticle() const
    837 {
    838   return secondaryParticle;
    839 }
    840 
    841 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    842 
    843 inline G4PhysicsTable* G4VEnergyLossProcess::DEDXTable() const
    844 {
    845   return theDEDXTable;
    846 }
    847 
    848 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    849 
    850 inline G4PhysicsTable* G4VEnergyLossProcess::DEDXTableForSubsec() const
    851 {
    852   return theDEDXSubTable;
    853 }
    854 
    855 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    856 
    857 inline G4PhysicsTable* G4VEnergyLossProcess::DEDXunRestrictedTable() const
    858 {
    859   return theDEDXunRestrictedTable;
    860 }
    861 
    862 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    863 
    864 inline G4PhysicsTable* G4VEnergyLossProcess::IonisationTable() const
    865 {
    866   G4PhysicsTable* t = theDEDXTable;
    867   if(theIonisationTable) t = theIonisationTable;
    868   return t;
    869 }
    870 
    871 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    872 
    873 inline G4PhysicsTable* G4VEnergyLossProcess::IonisationTableForSubsec() const
    874 {
    875   G4PhysicsTable* t = theDEDXSubTable;
    876   if(theIonisationSubTable) t = theIonisationSubTable;
    877   return t;
    878 }
    879 
    880 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    881 
    882 inline G4PhysicsTable* G4VEnergyLossProcess::CSDARangeTable() const
    883 {
    884   return theCSDARangeTable;
    885 }
    886 
    887 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    888 
    889 inline G4PhysicsTable* G4VEnergyLossProcess::RangeTableForLoss() const
    890 {
    891   return theRangeTableForLoss;
    892 }
    893 
    894 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    895 
    896 inline G4PhysicsTable* G4VEnergyLossProcess::InverseRangeTable() const
    897 {
    898   return theInverseRangeTable;
    899 }
    900 
    901 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    902 
    903 inline G4PhysicsTable* G4VEnergyLossProcess::LambdaTable()
    904 {
    905   return theLambdaTable;
    906 }
    907 
    908 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    909 
    910 inline G4PhysicsTable* G4VEnergyLossProcess::SubLambdaTable()
    911 {
    912   return theSubLambdaTable;
    913 }
    914 
    915 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    916  
    917 inline G4bool G4VEnergyLossProcess::IsIntegral() const
    918 {
    919   return integral;
    920 }
    921 
    922 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    923 
    924 inline size_t G4VEnergyLossProcess::CurrentMaterialCutsCoupleIndex() const
    925 {
    926   return currentMaterialIndex;
    927 }
    928 
    929 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    930 
    931 inline void G4VEnergyLossProcess::SetDynamicMassCharge(G4double massratio,
    932                                                        G4double charge2ratio)
    933 {
    934   massRatio     = massratio;
    935   chargeSqRatio = charge2ratio;
    936   reduceFactor  = 1.0/(chargeSqRatio*massRatio);
    937 }
    938 
    939 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    940  
    941 inline G4double G4VEnergyLossProcess::GetCurrentRange() const
    942 {
    943   return fRange;
    944 }
    945 
    946 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    947 
    948 inline
    949 void G4VEnergyLossProcess::AddEmModel(G4int order, G4VEmModel* p,
    950                                       G4VEmFluctuationModel* fluc,
    951                                       const G4Region* region)
    952 {
    953   modelManager->AddEmModel(order, p, fluc, region);
    954   if(p) p->SetParticleChange(pParticleChange, fluc);
    955 }
    956 
    957 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    958 
    959 inline
    960 G4VEmModel* G4VEnergyLossProcess::GetModelByIndex(G4int idx, G4bool ver)
    961 {
    962   return modelManager->GetModel(idx, ver);
    963 }
    964 
    965 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    966 
    967 inline G4int G4VEnergyLossProcess::NumberOfModels()
    968 {
    969   return modelManager->NumberOfModels();
    970 }
    971 
    972 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    973 
    974 inline void G4VEnergyLossProcess::SetEmModel(G4VEmModel* p, G4int index)
    975 {
    976   G4int n = emModels.size();
    977   if(index >= n) for(G4int i=n; i<index+1; i++) {emModels.push_back(0);}
    978   emModels[index] = p;
    979 }
    980 
    981 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    982 
    983 inline G4VEmModel* G4VEnergyLossProcess::EmModel(G4int index)
    984 {
    985   G4VEmModel* p = 0;
    986   if(index >= 0 && index <  G4int(emModels.size())) p = emModels[index];
    987   return p;
    988 }
    989 
    990 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    991 
    992 inline void G4VEnergyLossProcess::SetFluctModel(G4VEmFluctuationModel* p)
    993 {
    994   fluctModel = p;
    995 }
    996 
    997 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    998 
    999 inline G4VEmFluctuationModel* G4VEnergyLossProcess::FluctModel()
    1000 {
    1001   return fluctModel;
    1002 }
    1003 
    1004 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1005 
    1006 inline void G4VEnergyLossProcess::UpdateEmModel(const G4String& nam,
    1007                                                 G4double emin, G4double emax)
    1008 {
    1009   modelManager->UpdateEmModel(nam, emin, emax);
    1010 }
    1011 
    1012 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1013 
    1014 inline void G4VEnergyLossProcess::SetIntegral(G4bool val)
    1015 {
    1016   integral = val;
    1017 }
    1018 
    1019 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1020 
    1021 inline void G4VEnergyLossProcess::SetParticle(const G4ParticleDefinition* p)
    1022 {
    1023   particle = p;
    1024 }
    1025 
    1026 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1027 
    1028 inline void G4VEnergyLossProcess::SetBaseParticle(const G4ParticleDefinition* p)
    1029 {
    1030   baseParticle = p;
    1031 }
    1032 
    1033 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1034 
    1035 inline void G4VEnergyLossProcess::SetSecondaryParticle(const G4ParticleDefinition* p)
    1036 {
    1037   secondaryParticle = p;
    1038 }
    1039 
    1040 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1041 
    1042 inline void G4VEnergyLossProcess::SetLinearLossLimit(G4double val)
    1043 {
    1044   linLossLimit = val;
    1045 }
    1046 
    1047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1048 
    1049 inline void G4VEnergyLossProcess::SetLossFluctuations(G4bool val)
    1050 {
    1051   lossFluctuationFlag = val;
    1052 }
    1053 
    1054 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1055 
    1056 inline void G4VEnergyLossProcess::SetRandomStep(G4bool val)
    1057 {
    1058   rndmStepFlag = val;
    1059 }
    1060 
    1061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1062 
    1063 inline void G4VEnergyLossProcess::SetMinSubRange(G4double val)
    1064 {
    1065   minSubRange = val;
    1066 }
    1067 
    1068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1069 
    1070 inline G4bool G4VEnergyLossProcess::TablesAreBuilt() const
    1071 {
    1072   return  tablesAreBuilt;
    1073 }
    1074 
    1075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1076 
    1077 inline G4int G4VEnergyLossProcess::NumberOfSubCutoffRegions() const
    1078 {
    1079   return nSCoffRegions;
    1080 }
    1081 
    1082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1083 
    1084 inline void G4VEnergyLossProcess::SetDEDXBinning(G4int nbins)
    1085 {
    1086   nBins = nbins;
    1087 }
    1088 
    1089 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1090 
    1091 inline void G4VEnergyLossProcess::SetLambdaBinning(G4int nbins)
    1092 {
    1093   nBins = nbins;
    1094 }
    1095 
    1096 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1097 
    1098 inline void G4VEnergyLossProcess::SetDEDXBinningForCSDARange(G4int nbins)
    1099 {
    1100   nBinsCSDA = nbins;
    1101 }
    1102 
    1103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1104 
    1105 inline G4double G4VEnergyLossProcess::MinKinEnergy() const
    1106 {
    1107   return minKinEnergy;
    1108 }
    1109 
    1110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1111 
    1112 inline void G4VEnergyLossProcess::SetMinKinEnergy(G4double e)
    1113 {
    1114   minKinEnergy = e;
    1115 }
    1116 
    1117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1118 
    1119 inline void G4VEnergyLossProcess::SetMaxKinEnergy(G4double e)
    1120 {
    1121   maxKinEnergy = e;
    1122   if(e < maxKinEnergyCSDA) maxKinEnergyCSDA = e;
    1123 }
    1124 
    1125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1126 
    1127 inline void G4VEnergyLossProcess::SetMaxKinEnergyForCSDARange(G4double e)
    1128 {
    1129   maxKinEnergyCSDA = e;
    1130 }
    1131 
    1132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1133 
    1134 inline G4double G4VEnergyLossProcess::MaxKinEnergy() const
    1135 {
    1136   return maxKinEnergy;
    1137 }
    1138 
    1139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1140 
    1141 inline void G4VEnergyLossProcess::SetLambdaFactor(G4double val)
    1142 {
    1143   if(val > 0.0 && val <= 1.0) lambdaFactor = val;
    1144 }
    1145 
    1146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1147 
    1148 inline void G4VEnergyLossProcess::SetIonisation(G4bool val)
    1149 {
    1150   isIonisation = val;
    1151   if(val) aGPILSelection = CandidateForSelection;
    1152   else    aGPILSelection = NotCandidateForSelection;
    1153 }
    1154 
    1155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1156 
    1157 inline G4bool G4VEnergyLossProcess::IsIonisationProcess() const
    1158 {
    1159   return isIonisation;
    1160 }
    1161 
    1162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1163 
    1164 void G4VEnergyLossProcess::SetStepFunction(G4double v1, G4double v2)
    1165 {
    1166   dRoverRange = v1;
    1167   finalRange = v2;
    1168   if (dRoverRange > 0.999) dRoverRange = 1.0;
    1169   currentCouple = 0;
    1170   mfpKinEnergy  = DBL_MAX;
    1171 }
    1172 
    1173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1174 
    11751120#endif

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMscModel.hh,v 1.4 2008/03/10 10:39:28 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VMscModel.hh,v 1.9 2009/04/07 18:39:47 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    3939//
    4040// Modifications:
    41 //
     41// 07.04.2009 V.Ivanchenko moved msc methods from G4VEmModel to G4VMscModel
    4242//
    4343// Class Description:
    … …  
    5454#include "G4MscStepLimitType.hh"
    5555#include "globals.hh"
     56#include "G4ThreeVector.hh"
     57#include "G4Track.hh"
     58#include "G4SafetyHelper.hh"
     59
     60class G4ParticleChangeForMSC;
    5661
    5762class G4VMscModel : public G4VEmModel
    … …  
    6469  virtual ~G4VMscModel();
    6570
     71  virtual G4double ComputeTruePathLengthLimit(const G4Track& track,
     72                                              G4PhysicsTable* theLambdaTable,
     73                                              G4double currentMinimalStep);
     74
     75  virtual G4double ComputeGeomPathLength(G4double truePathLength);
     76
     77  virtual G4double ComputeTrueStepLength(G4double geomPathLength);
     78
     79  virtual void SampleScattering(const G4DynamicParticle*,
     80                                G4double safety);
     81
     82  // empty
     83  virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
     84                                 const G4MaterialCutsCouple*,
     85                                 const G4DynamicParticle*,
     86                                 G4double tmin,
     87                                 G4double tmax);
     88
     89  //================================================================
     90  //  Set parameters of multiple scattering models
     91  //================================================================
     92 
    6693  inline void SetStepLimitType(G4MscStepLimitType);
    6794
    … …  
    73100
    74101  inline void SetSkin(G4double);
     102
     103  inline void SetSampleZ(G4bool);
     104
     105protected:
     106
     107  // initialisation of the ParticleChange for the model
     108  G4ParticleChangeForMSC* GetParticleChangeForMSC();
     109
     110  // initialisation of interface with geometry
     111  void InitialiseSafetyHelper();
     112
     113  // shift point of the track PostStep
     114  void ComputeDisplacement(G4ParticleChangeForMSC*, 
     115                           const G4ThreeVector& displDir,
     116                           G4double displacement,
     117                           G4double postsafety);
     118
     119  // compute safety
     120  inline G4double ComputeSafety(const G4ThreeVector& position, G4double limit);
     121
     122  // compute linear distance to a geometry boundary
     123  inline G4double ComputeGeomLimit(const G4Track& position, G4double& presafety,
     124                                   G4double limit);
    75125
    76126private:
    … …  
    79129  G4VMscModel & operator=(const  G4VMscModel &right);
    80130  G4VMscModel(const  G4VMscModel&);
     131
     132  G4SafetyHelper* safetyHelper;
    81133
    82134protected:
    … …  
    88140  G4double dtrl;
    89141  G4double lambdalimit;
     142  G4double geommax;
    90143
    91144  G4MscStepLimitType steppingAlgorithm;
    … …  
    134187//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    135188
     189inline void G4VMscModel::SetSampleZ(G4bool val)
     190{
     191  samplez = val;
     192}
     193
     194//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     195
     196inline G4double G4VMscModel::ComputeSafety(const G4ThreeVector& position,
     197                                           G4double)
     198{
     199  return safetyHelper->ComputeSafety(position);
     200}
     201
     202//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     203
     204inline G4double G4VMscModel::ComputeGeomLimit(const G4Track& track,
     205                                              G4double& presafety,
     206                                              G4double limit)
     207{
     208  G4double res = geommax;
     209  if(track.GetVolume() != safetyHelper->GetWorldVolume()) {
     210    res = safetyHelper->CheckNextStep(
     211          track.GetStep()->GetPreStepPoint()->GetPosition(),
     212          track.GetMomentumDirection(),
     213          limit, presafety);
     214  }
     215  return res;
     216}
     217
     218//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     219
    136220#endif
    137221

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMultipleScattering.hh,v 1.54 2008/07/31 13:01:26 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VMultipleScattering.hh,v 1.56 2009/04/07 18:39:47 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6464// 27-10-07 Virtual functions moved to source (V.Ivanchenko)
    6565// 15-07-08 Reorder class members for further multi-thread development (VI)
     66// 07-04-09 Moved msc methods from G4VEmModel to G4VMscModel (VI)
    6667//
    6768
    … …  
    8081#include "G4Step.hh"
    8182#include "G4EmModelManager.hh"
    82 #include "G4VEmModel.hh"
     83#include "G4VMscModel.hh"
    8384#include "G4MscStepLimitType.hh"
    8485
    … …  
    145146                              G4bool ascii);
    146147
    147   //------------------------------------------------------------------------
    148   // Specific methods for msc processes
    149   //------------------------------------------------------------------------
    150 
    151148  // The function overloads the corresponding function of the base
    152149  // class.It limits the step near to boundaries only
    153150  // and invokes the method GetMscContinuousStepLimit at every step.
    154   virtual G4double AlongStepGetPhysicalInteractionLength(
     151  G4double AlongStepGetPhysicalInteractionLength(
    155152                                            const G4Track&,
    156153                                            G4double  previousStepSize,
    … …  
    192189  inline G4PhysicsTable* LambdaTable() const;
    193190
    194   //------------------------------------------------------------------------
    195   // Define and access particle type
    196   //------------------------------------------------------------------------
    197 
     191  // access particle type
    198192  inline const G4ParticleDefinition* Particle() const;
    199   inline void SetParticle(const G4ParticleDefinition*);
    200193
    201194  //------------------------------------------------------------------------
    … …  
    203196  //------------------------------------------------------------------------
    204197
    205   inline void AddEmModel(G4int, G4VEmModel*, const G4Region* region = 0);
    206 
     198protected:
     199  // Select model in run time
     200  inline G4VEmModel* SelectModel(G4double kinEnergy);
     201
     202public:
     203  // Select model in run time
    207204  inline G4VEmModel* SelectModelForMaterial(G4double kinEnergy,
    208205                                            size_t& idxRegion) const;
    209206
    210   // Access to models
    211   inline G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false);
    212 
    213   //------------------------------------------------------------------------
    214   // Set parameters for simulation of multiple scattering
    215   //------------------------------------------------------------------------
    216 
     207  // Add model for region, smaller value of order defines which
     208  // model will be selected for a given energy interval 
     209  void AddEmModel(G4int order, G4VEmModel*, const G4Region* region = 0);
     210
     211  // Access to models by index
     212  G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false);
     213
     214  //------------------------------------------------------------------------
     215  // Get/Set parameters for simulation of multiple scattering
     216  //------------------------------------------------------------------------
     217
     218  inline G4bool LateralDisplasmentFlag() const;
    217219  inline void SetLateralDisplasmentFlag(G4bool val);
    218220
     221  inline G4double Skin() const;
    219222  inline void SetSkin(G4double val);
    220223
     224  inline G4double RangeFactor() const;
    221225  inline void SetRangeFactor(G4double val);
    222226
     227  inline G4double GeomFactor() const;
    223228  inline void SetGeomFactor(G4double val);
    224229
     230  inline G4double PolarAngleLimit() const;
    225231  inline void SetPolarAngleLimit(G4double val);
    226232
     233  inline G4MscStepLimitType StepLimitType() const;
    227234  inline void SetStepLimitType(G4MscStepLimitType val);
    228235
     236  //------------------------------------------------------------------------
     237  // Run time methods
     238  //------------------------------------------------------------------------
     239
    229240protected:
    230241
    231   // This method is used for tracking, it returns mean free path value
     242  // This method is not used for tracking, it returns mean free path value
    232243  G4double GetMeanFreePath(const G4Track& track,
    233244                           G4double,
    234245                           G4ForceCondition* condition);
    235 
    236   //------------------------------------------------------------------------
    237   // Run time methods
    238   //------------------------------------------------------------------------
    239246
    240247  // This method is not used for tracking, it returns step limit
    … …  
    244251                                  G4double& currentSafety);
    245252
     253  // This method returns inversed transport cross section
    246254  inline G4double GetLambda(const G4ParticleDefinition* p,
    247255                            G4double& kineticEnergy);
    … …  
    253261                                            G4double& currentSafety);
    254262
    255   inline G4VEmModel* SelectModel(G4double kinEnergy);
    256   // Select concrete model
     263  // defines current material in run time
     264  inline void DefineMaterial(const G4MaterialCutsCouple* couple);
    257265
    258266  inline const G4MaterialCutsCouple* CurrentMaterialCutsCouple() const;
    259267
    260   // define current material
    261   inline void DefineMaterial(const G4MaterialCutsCouple* couple);
    262 
    263   //------------------------------------------------------------------------
    264   // Access parameters of multiple scattering
    265   //------------------------------------------------------------------------
    266 
    267   inline G4ParticleChangeForMSC* GetParticleChange();
    268 
    269   inline G4double Skin() const;
    270 
    271   inline G4double RangeFactor() const;
    272 
    273   inline G4double GeomFactor() const;
    274 
    275   inline G4double PolarAngleLimit() const;
    276 
    277   inline G4MscStepLimitType StepLimitType() const;
    278 
    279   inline G4bool LateralDisplasmentFlag() const;
    280 
    281268private:
    282269
    283270  // hide  assignment operator
    284 
    285271  G4VMultipleScattering(G4VMultipleScattering &);
    286272  G4VMultipleScattering & operator=(const G4VMultipleScattering &right);
    … …  
    318304private:
    319305
    320   G4VEmModel*                 currentModel;
     306  G4VMscModel*                currentModel;
    321307
    322308  // cache
    … …  
    330316//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    331317
     318inline G4double G4VMultipleScattering::ContinuousStepLimit(
     319                                       const G4Track& track,
     320                                       G4double previousStepSize,
     321                                       G4double currentMinimalStep,
     322                                       G4double& currentSafety)
     323{
     324  return GetMscContinuousStepLimit(track,previousStepSize,currentMinimalStep,
     325                                   currentSafety);
     326}
     327
     328//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     329
     330inline void G4VMultipleScattering::SetBinning(G4int nbins)
     331{
     332  nBins = nbins;
     333}
     334
     335//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     336
     337inline G4int G4VMultipleScattering::Binning() const
     338{
     339  return nBins;
     340}
     341
     342//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     343
     344inline void G4VMultipleScattering::SetMinKinEnergy(G4double e)
     345{
     346  minKinEnergy = e;
     347}
     348
     349//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     350
     351inline G4double G4VMultipleScattering::MinKinEnergy() const
     352{
     353  return minKinEnergy;
     354}
     355
     356//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     357
     358inline void G4VMultipleScattering::SetMaxKinEnergy(G4double e)
     359{
     360  maxKinEnergy = e;
     361}
     362
     363//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     364
     365inline G4double G4VMultipleScattering::MaxKinEnergy() const
     366{
     367  return maxKinEnergy;
     368}
     369
     370//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     371
     372inline  void G4VMultipleScattering::SetBuildLambdaTable(G4bool val)
     373{
     374  buildLambdaTable = val;
     375}
     376
     377//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     378
     379inline G4PhysicsTable* G4VMultipleScattering::LambdaTable() const
     380{
     381  return theLambdaTable;
     382}
     383
     384//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     385
     386inline  const G4ParticleDefinition* G4VMultipleScattering::Particle() const
     387{
     388  return currentParticle;
     389}
     390
     391//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     392
     393inline G4VEmModel* G4VMultipleScattering::SelectModel(G4double kinEnergy)
     394{
     395  return modelManager->SelectModel(kinEnergy, currentMaterialIndex);
     396}
     397
     398//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     399
     400inline G4VEmModel* G4VMultipleScattering::SelectModelForMaterial(
     401                   G4double kinEnergy, size_t& idxRegion) const
     402{
     403  return modelManager->SelectModel(kinEnergy, idxRegion);
     404}
     405
     406//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     407
     408inline  G4bool G4VMultipleScattering::LateralDisplasmentFlag() const
     409{
     410  return latDisplasment;
     411}
     412
     413//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     414
     415inline  void G4VMultipleScattering::SetLateralDisplasmentFlag(G4bool val)
     416{
     417  latDisplasment = val;
     418}
     419
     420//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     421
     422inline  G4double G4VMultipleScattering::Skin() const
     423{
     424  return skin;
     425}
     426
     427//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     428
     429inline  void G4VMultipleScattering::SetSkin(G4double val)
     430{
     431  if(val < 1.0) skin = 0.0;
     432  else          skin = val;
     433}
     434
     435//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     436
     437inline  G4double G4VMultipleScattering::RangeFactor() const
     438{
     439  return facrange;
     440}
     441
     442//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     443
     444inline  void G4VMultipleScattering::SetRangeFactor(G4double val)
     445{
     446  if(val > 0.0) facrange = val;
     447}
     448
     449//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     450
     451inline  G4double G4VMultipleScattering::GeomFactor() const
     452{
     453  return facgeom;
     454}
     455
     456//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     457
     458inline  void G4VMultipleScattering::SetGeomFactor(G4double val)
     459{
     460  if(val > 0.0) facgeom = val;
     461}
     462
     463//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     464
     465inline  G4double G4VMultipleScattering::PolarAngleLimit() const
     466{
     467  return polarAngleLimit;
     468}
     469
     470//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     471
     472inline  void G4VMultipleScattering::SetPolarAngleLimit(G4double val)
     473{
     474  if(val < 0.0)     polarAngleLimit = 0.0;
     475  else if(val > pi) polarAngleLimit = pi;
     476  else              polarAngleLimit = val;
     477}
     478
     479//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     480
     481inline G4MscStepLimitType G4VMultipleScattering::StepLimitType() const
     482{
     483  return stepLimit;
     484}
     485
     486//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     487
     488inline void G4VMultipleScattering::SetStepLimitType(G4MscStepLimitType val)
     489{
     490  stepLimit = val;
     491  if(val == fMinimal) facrange = 0.2;
     492}
     493
     494//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     495
    332496inline
    333 void G4VMultipleScattering::DefineMaterial(const G4MaterialCutsCouple* couple)
    334 {
    335   if(couple != currentCouple) {
    336     currentCouple   = couple;
    337     currentMaterialIndex = couple->GetIndex();
     497G4double G4VMultipleScattering::GetLambda(const G4ParticleDefinition* p,
     498                                          G4double& e)
     499{
     500  G4double x;
     501  if(theLambdaTable) {
     502    G4bool b;
     503    x = ((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b);
     504  } else {
     505    x = currentModel->CrossSection(currentCouple,p,e);
    338506  }
     507  if(x > DBL_MIN) x = 1./x;
     508  else            x = DBL_MAX;
     509  return x;
    339510}
    340511
    … …  
    349520  G4double x = currentMinimalStep;
    350521  DefineMaterial(track.GetMaterialCutsCouple());
    351   currentModel = SelectModel(scaledKinEnergy);
     522  currentModel = static_cast<G4VMscModel*>(SelectModel(scaledKinEnergy));
    352523  if(x > 0.0 && scaledKinEnergy > 0.0) {
    353524    G4double tPathLength =
    … …  
    364535//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    365536
    366 inline G4double G4VMultipleScattering::ContinuousStepLimit(
    367                                        const G4Track& track,
    368                                        G4double previousStepSize,
    369                                        G4double currentMinimalStep,
    370                                        G4double& currentSafety)
    371 {
    372   return GetMscContinuousStepLimit(track,previousStepSize,currentMinimalStep,
    373                                    currentSafety);
    374 }
    375 
    376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    377 
    378537inline
    379 G4double G4VMultipleScattering::GetLambda(const G4ParticleDefinition* p,
    380                                           G4double& e)
    381 {
    382   G4double x;
    383   if(theLambdaTable) {
    384     G4bool b;
    385     x = ((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b);
    386   } else {
    387     x = currentModel->CrossSection(currentCouple,p,e);
     538void G4VMultipleScattering::DefineMaterial(const G4MaterialCutsCouple* couple)
     539{
     540  if(couple != currentCouple) {
     541    currentCouple   = couple;
     542    currentMaterialIndex = couple->GetIndex();
    388543  }
    389   if(x > DBL_MIN) x = 1./x;
    390   else            x = DBL_MAX;
    391   return x;
    392 }
    393 
    394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    395 
    396 inline G4VEmModel* G4VMultipleScattering::SelectModel(G4double kinEnergy)
    397 {
    398   return modelManager->SelectModel(kinEnergy, currentMaterialIndex);
    399 }
    400 
    401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    402 
    403 inline G4VEmModel* G4VMultipleScattering::SelectModelForMaterial(
    404                    G4double kinEnergy, size_t& idxRegion) const
    405 {
    406   return modelManager->SelectModel(kinEnergy, idxRegion);
    407 }
    408 
    409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    410 
    411 inline void G4VMultipleScattering::SetBinning(G4int nbins)
    412 {
    413   nBins = nbins;
    414 }
    415 
    416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    417 
    418 inline G4int G4VMultipleScattering::Binning() const
    419 {
    420   return nBins;
    421 }
    422 
    423 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    424 
    425 inline void G4VMultipleScattering::SetMinKinEnergy(G4double e)
    426 {
    427   minKinEnergy = e;
    428 }
    429 
    430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    431 
    432 inline G4double G4VMultipleScattering::MinKinEnergy() const
    433 {
    434   return minKinEnergy;
    435 }
    436 
    437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    438 
    439 inline void G4VMultipleScattering::SetMaxKinEnergy(G4double e)
    440 {
    441   maxKinEnergy = e;
    442 }
    443 
    444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    445 
    446 inline G4double G4VMultipleScattering::MaxKinEnergy() const
    447 {
    448   return maxKinEnergy;
    449 }
    450 
    451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    452 
    453 inline  G4bool G4VMultipleScattering::LateralDisplasmentFlag() const
    454 {
    455   return latDisplasment;
    456 }
    457 
    458 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    459 
    460 inline  void G4VMultipleScattering::SetLateralDisplasmentFlag(G4bool val)
    461 {
    462   latDisplasment = val;
    463 }
    464 
    465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    466 
    467 inline G4ParticleChangeForMSC* G4VMultipleScattering::GetParticleChange()
    468 {
    469   return &fParticleChange;
    470 }
    471 
    472 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    473 
    474 inline  G4double G4VMultipleScattering::Skin() const
    475 {
    476   return skin;
    477 }
    478 
    479 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    480 
    481 inline  void G4VMultipleScattering::SetSkin(G4double val)
    482 {
    483   if(val < 1.0) skin = 0.0;
    484   else          skin = val;
    485 }
    486 
    487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    488 
    489 inline  G4double G4VMultipleScattering::RangeFactor() const
    490 {
    491   return facrange;
    492 }
    493 
    494 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    495 
    496 inline  void G4VMultipleScattering::SetRangeFactor(G4double val)
    497 {
    498   if(val > 0.0) facrange = val;
    499 }
    500 
    501 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    502 
    503 inline  G4double G4VMultipleScattering::GeomFactor() const
    504 {
    505   return facgeom;
    506 }
    507 
    508 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    509 
    510 inline  void G4VMultipleScattering::SetGeomFactor(G4double val)
    511 {
    512   if(val > 0.0) facgeom = val;
    513 }
    514 
    515 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    516 
    517 inline  G4double G4VMultipleScattering::PolarAngleLimit() const
    518 {
    519   return polarAngleLimit;
    520 }
    521 
    522 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    523 
    524 inline  void G4VMultipleScattering::SetPolarAngleLimit(G4double val)
    525 {
    526   if(val < 0.0)     polarAngleLimit = 0.0;
    527   else if(val > pi) polarAngleLimit = pi;
    528   else              polarAngleLimit = val;
    529 }
    530 
    531 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    532 
    533 inline G4MscStepLimitType G4VMultipleScattering::StepLimitType() const
    534 {
    535   return stepLimit;
    536 }
    537 
    538 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    539 
    540 inline void G4VMultipleScattering::SetStepLimitType(G4MscStepLimitType val)
    541 {
    542   stepLimit = val;
    543   if(val == fMinimal) facrange = 0.2;
    544 }
    545 
    546 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    547 
    548 inline  void G4VMultipleScattering::SetBuildLambdaTable(G4bool val)
    549 {
    550   buildLambdaTable = val;
    551 }
    552 
    553 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    554 
    555 inline  const G4ParticleDefinition* G4VMultipleScattering::Particle() const
    556 {
    557   return currentParticle;
    558 }
    559 
    560 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    561 
    562 inline G4PhysicsTable* G4VMultipleScattering::LambdaTable() const
    563 {
    564   return theLambdaTable;
    565 }
    566 
    567 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    568 
    569 inline
    570 const G4MaterialCutsCouple* G4VMultipleScattering::CurrentMaterialCutsCouple() const
     544}
     545
     546//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     547
     548inline const G4MaterialCutsCouple*
     549G4VMultipleScattering::CurrentMaterialCutsCouple() const
    571550{
    572551  return currentCouple;
    … …  
    575554//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    576555
    577 inline void G4VMultipleScattering::AddEmModel(G4int order, G4VEmModel* p,
    578                                               const G4Region* region)
    579 {
    580   G4VEmFluctuationModel* fm = 0;
    581   modelManager->AddEmModel(order, p, fm, region);
    582   if(p) p->SetParticleChange(pParticleChange);
    583 }
    584 
    585 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    586 
    587 inline
    588 G4VEmModel* G4VMultipleScattering::GetModelByIndex(G4int idx, G4bool ver)
    589 {
    590   return modelManager->GetModel(idx, ver);
    591 }
    592 
    593 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    594 
    595556#endif

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DummyModel.cc,v 1.3 2007/05/22 17:31:58 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4DummyModel.cc,v 1.4 2009/04/07 18:39:47 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5252
    5353G4DummyModel::G4DummyModel(const G4String& nam)
    54   : G4VEmModel(nam)
     54  : G4VMscModel(nam)
    5555{}
    5656

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmCalculator.cc,v 1.44 2008/08/03 18:47:15 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4EmCalculator.cc,v 1.46 2009/02/24 09:56:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    755755    if(currentProcess) currentProcessName = currentProcess->GetProcessName();
    756756
    757     if(p->GetParticleType() == "nucleus" &&
    758        currentParticleName  != "deuteron" && currentParticleName != "triton") {
     757    if(p->GetParticleType() == "nucleus"
     758       && currentParticleName != "deuteron" 
     759       && currentParticleName != "triton"
     760       && currentParticleName != "alpha+"
     761       && currentParticleName != "helium"
     762       && currentParticleName != "hydrogen"
     763      ) {
    759764      baseParticle = theGenericIon;
    760765      massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmConfigurator.cc,v 1.3 2008/11/21 12:30:29 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4EmConfigurator.cc,v 1.4 2009/02/26 11:33:33 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    117117  G4String fname = "";
    118118  if(fm) fname = fm->GetName();
    119   AddModelForRegion(particleName, processName, mod->GetName(), regionName,
     119  G4String mname = "";
     120  if(mod) mname = mod->GetName();
     121  AddModelForRegion(particleName, processName, mname, regionName,
    120122                    emin, emax, fname);
    121123}
    … …  
    203205
    204206        for(G4int i=0; i<nm; i++) {
    205           if(modelName == modelList[i]->GetName() &&
     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 &&
    206212             (particleList[i] == "" || particleList[i] == particleName) ) {
    207213            mod  = modelList[i];
    … …  
    214220
    215221        if(!mod) {
    216           G4cout << "### G4EmConfigurator WARNING: fails to find a model <"
    217                  << modelName << "> for process <"
    218                  << processName << "> and " << particleName
    219                  << G4endl;
    220           if(flucModelName != "") 
    221             G4cout << "                            fluctuation model <"
    222                    << flucModelName << G4endl;
     222
     223          // set fluctuation model for ionisation processes
     224          if(fluc && ptype == eloss) {
     225            G4VEnergyLossProcess* p = reinterpret_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          }
    223238        } else {
    224239

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmElementSelector.cc,v 1.4 2008/08/21 18:53:32 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4EmElementSelector.cc,v 1.10 2009/05/26 16:59:35 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-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)
    … …  
    6565  nElmMinusOne = n - 1;
    6666  theElementVector = material->GetElementVector();
     67  element = (*theElementVector)[0];
    6768  if(nElmMinusOne > 0) {
    68     for(G4int i=0; i<nElmMinusOne; i++) {
     69    xSections.reserve(n);
     70    for(G4int i=0; i<n; i++) {
    6971      G4PhysicsLogVector* v = new G4PhysicsLogVector(lowEnergy,highEnergy,nbins);
    70       v->SetSpline(spline);
     72      //v->SetSpline(spline);
    7173      xSections.push_back(v);
    7274    }
    … …  
    8082{
    8183  if(nElmMinusOne > 0) {
    82     for(G4int i=0; i<nElmMinusOne; i++) {
     84    for(G4int i=0; i<=nElmMinusOne; i++) {
    8385      delete xSections[i];
    8486    }
    … …  
    101103
    102104  G4int i;
    103   G4int n = nElmMinusOne + 1;
    104   G4double* xsec = new G4double[n]; 
    105105
    106106  // loop over bins
    … …  
    110110    cross = 0.0;
    111111    //G4cout << "j= " << j << " e(MeV)= " << e/MeV << G4endl;
    112     for (i=0; i<n; i++) {
     112    for (i=0; i<=nElmMinusOne; i++) {
    113113      cross += theAtomNumDensityVector[i]*     
    114114        model->ComputeCrossSectionPerAtom(part, (*theElementVector)[i], e,
    115115                                          cutEnergy, e);
    116       xsec[i] = cross;
    117     }
    118     if(DBL_MIN >= cross) cross = 1.0;
    119     // normalise cross section sum
    120     for (i=0; i<nElmMinusOne; i++) {
    121       xSections[i]->PutValue(j, xsec[i]/cross);
    122       //G4cout << "i= " << i << " xs= " << xsec[i]/cross << G4endl;
     116      xSections[i]->PutValue(j, cross);
    123117    }
    124118  }
    125   delete [] xsec;
     119
     120  // 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++) {
     123      xSections[i]->PutValue(0, (*xSections[i])[1]);
     124    }
     125  }
     126  // xSections ends with null, so use probabilities from the previous bin
     127  if(DBL_MIN >= (*xSections[nElmMinusOne])[nbins-1]) {
     128    for (i=0; i<=nElmMinusOne; i++) {
     129      xSections[i]->PutValue(nbins-1, (*xSections[i])[nbins-2]);
     130    }
     131  }
     132  // perform normalization
     133  for(G4int j=0; j<nbins; j++) {
     134    cross = (*xSections[nElmMinusOne])[j];
     135    // only for positive X-section
     136    if(cross > DBL_MIN) {
     137      for (i=0; i<nElmMinusOne; i++) {
     138        xSections[i]->PutValue(j, (*xSections[i])[j]/cross);
     139      }
     140    }
     141  }
    126142}
    127143
    … …  
    139155    }
    140156  } 
    141   G4cout << "Last Element in element vector"
     157  G4cout << "Last Element in element vector "
    142158         << (*theElementVector)[nElmMinusOne]->GetName()
    143159         << G4endl;

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmModelManager.cc,v 1.46 2008/10/13 14:56:56 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4EmModelManager.cc,v 1.49 2009/04/17 10:35:32 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6161// 15-03-07 Add maxCutInRange (V.Ivanchenko)
    6262// 12-04-07 Add verbosity at destruction (V.Ivanchenko)
     63// 08-04-08 Fixed and simplified initialisation of G4RegionModel (VI)
    6364//
    6465// Class Description:
    … …  
    134135  theGamma = G4Gamma::Gamma();
    135136  thePositron = G4Positron::Positron();
     137  models.reserve(4);
     138  flucModels.reserve(4);
     139  regions.reserve(4);
     140  orderOfModels.reserve(4);
     141  isUsed.reserve(4);
    136142}
    137143
    … …  
    178184  regions.push_back(r);
    179185  orderOfModels.push_back(num);
     186  isUsed.push_back(0);
    180187  p->DefineForRegion(r);
    181188  nEmModels++;
    … …  
    187194                                     G4double emin, G4double emax)
    188195{
    189   if (nEmModels) {
     196  if (nEmModels > 0) {
    190197    for(G4int i=0; i<nEmModels; i++) {
    191198      if(nam == models[i]->GetName()) {
    … …  
    225232{
    226233  verboseLevel = val;
     234  G4String partname = p->GetParticleName();
    227235  if(1 < verboseLevel) {
    228236    G4cout << "G4EmModelManager::Initialise() for "
    229            << p->GetParticleName()
    230            << G4endl;
     237           << partname << G4endl;
    231238  }
    232239  // Are models defined?
    233240  if(!nEmModels) {
    234     G4Exception("G4EmModelManager::Initialise without any model defined");
     241    G4Exception("G4EmModelManager::Initialise without any model defined for "+partname);
    235242  }
    236243  particle = p;
    … …  
    271278    G4ProductionCutsTable::GetProductionCutsTable();
    272279  G4int numOfCouples = theCoupleTable->GetTableSize();
    273   idxOfRegionModels = new G4int[numOfCouples+1];
    274   idxOfRegionModels[numOfCouples] = 0;
     280  if(nRegions > 1) idxOfRegionModels = new G4int[numOfCouples];
    275281  setOfRegionModels = new G4RegionModels*[nRegions];
    276282
    277283  std::vector<G4int>    modelAtRegion(nEmModels);
    278284  std::vector<G4int>    modelOrd(nEmModels);
    279   G4DataVector          eLow(nEmModels);
     285  G4DataVector          eLow(nEmModels+1);
    280286  G4DataVector          eHigh(nEmModels);
    281   G4int nmax = nEmModels;
    282287
    283288  // Order models for regions
    … …  
    305310            if (region) G4cout << region->GetName();
    306311            G4cout << ">  "
    307                  << " tmin(MeV)= " << tmin/MeV
    308                  << "; tmax(MeV)= " << tmax/MeV
    309                  << "; order= " << ord
    310                  << G4endl;
     312                   << " tmin(MeV)= " << tmin/MeV
     313                   << "; tmax(MeV)= " << tmax/MeV
     314                   << "; order= " << ord
     315                   << G4endl;
    311316          }
    312317       
    313           if (n == 0) n++;
    314           else {
     318          if(n > 0) {
     319
     320            // extend energy range to previous models
    315321            tmin = std::min(tmin, eHigh[n-1]);
    316             if(tmin >= tmax) push = false;
    317             else {
    318 
    319               // high energy model
    320               if(tmin == eHigh[n-1] && tmax > eHigh[n-1]) n++;
    321               else if (tmax > eHigh[n-1]) {
    322 
    323                 // compare order of models
    324                 for(G4int k = n-1; k>=0; k--) {
    325        
    326                   if (ord >= modelOrd[k]) {
    327                     tmin = std::max(tmin, eHigh[k]);           
    328                     if(k < n-1) n = k + 2;
    329                     break;
    330                   } else if (tmin > eLow[k]) {
    331                     eHigh[k] = tmin;
    332                     n = k + 2;
    333                     break;
    334                   } else if (tmin == eLow[k]) {
    335                     n = k + 1;
    336                     break;
     322            tmax = std::max(tmax, eLow[0]);
     323            //G4cout << "tmin= " << tmin << "  tmax= "
     324            //     << tmax << "  ord= " << ord <<G4endl;
     325            // empty energy range
     326            if( tmax - tmin <= eV) push = false;
     327            // low-energy model
     328            else if (tmax == eLow[0]) {
     329              push = false;
     330              insert = true;
     331              idx = 0;
     332              // resolve intersections
     333            } else if(tmin < eHigh[n-1]) {
     334              // compare order
     335              for(G4int k=0; k<n; k++) {
     336                // new model has lower application
     337                if(ord >= modelOrd[k]) {
     338                  if(tmin < eHigh[k]  && tmin >= eLow[k]) tmin = eHigh[k];
     339                  if(tmax <= eHigh[k] && tmax >  eLow[k]) tmax = eLow[k];
     340                  if(tmax > eHigh[k] && tmin < eLow[k]) {
     341                    if(tmax - eHigh[k] > eLow[k] - tmin) tmin = eHigh[k];
     342                    else tmax = eLow[k];
     343                  }
     344                  if( tmax - tmin <= eV) {
     345                    push = false;
     346                    break;
    337347                  }
    338348                }
    339                 if(tmin < eLow[0]) n = 1;
    340                 idx = n - 1;
    341 
    342                 // low energy model
    343               } else {
    344 
    345                 tmax = std::max(tmax, eLow[0]);
    346                 insert = true;
    347                 push = false;
    348                 idx = 0;
    349                 if(tmax <= eLow[0]) tmax = eLow[0];
    350                 else {
    351 
    352                   for(G4int k=0; k<n; k++) {
    353        
    354                     if (ord >= modelOrd[k]) {
    355                       if(k == 0) {
    356                         if(tmin < eLow[0]) tmax = eLow[0];
    357                         else insert = false;
    358                         break;
    359                       } else {
    360                         insert = false;
    361                         break;
    362                       }
    363                     } else if(tmax < eHigh[k]) {
    364                       idx = k;
    365                       if(k > 0) tmin = eLow[k];
    366                       eLow[k] = tmax;               
    367                       break;
    368                     } else if(tmax == eHigh[k]) {           
    369                       insert = false;
    370                       push = true;
    371                       idx = k;
    372                       if(k > 0) tmin = eLow[k];
    373                       else tmin = std::min(tmin,eLow[0]);
    374                       break;
    375                     } else {
    376                       modelAtRegion[k] = ii;
    377                       modelOrd[k] = ord;
    378                       if(k == 0) eLow[idx] = std::min(tmin,eLow[0]);
     349              }
     350              //G4cout << "tmin= " << tmin << "  tmax= "
     351              //     << tmax << "  push= " << push << " idx= " << idx <<G4endl;
     352              if(push) {
     353                if (tmax == eLow[0]) {
     354                  push = false;
     355                  insert = true;
     356                  idx = 0;
     357                  // continue resolve intersections
     358                } else if(tmin < eHigh[n-1]) {
     359                  // last energy interval
     360                  if(tmin > eLow[n-1] && tmax >= eHigh[n-1]) {
     361                    eHigh[n-1] = tmin;
     362                    // first energy interval
     363                  } else if(tmin <= eLow[0] && tmax < eHigh[0]) {
     364                    eLow[0] = tmax;
     365                    push = false;
     366                    insert = true;
     367                    idx = 0;
     368                  } else {
     369                    // find energy interval to replace
     370                    for(G4int k=0; k<n; k++) {
     371                      if(tmin <= eLow[k] && tmax >= eHigh[k]) {
     372                        push = false;
     373                        modelAtRegion[k] = ii;
     374                        modelOrd[k] = ord;
     375                        isUsed[ii] = 1;
     376                      }
    379377                    }
    380378                  }
    381379                }
    382                 if(insert && idx < n) n++;
    383                 else insert = false;
    384380              }
    385381            }
    386382          }
    387           if(n > nmax) {
    388             nmax = n;
    389             modelAtRegion.resize(nmax);
    390             modelOrd.resize(nmax);
    391             eLow.resize(nmax);
    392             eHigh.resize(nmax);
    393           }
    394383          if(insert) {
    395             for(G4int k=n-2; k>=idx; k--) {           
     384            for(G4int k=n-1; k>=idx; k--) {           
    396385              modelAtRegion[k+1] = modelAtRegion[k];
    397386              modelOrd[k+1] = modelOrd[k];
    … …  
    400389            }
    401390          }
     391          //G4cout << "push= " << push << " insert= " << insert
     392          //<< " idx= " << idx <<G4endl;
    402393          if (push || insert) {
     394            n++;
    403395            modelAtRegion[idx] = ii;
    404396            modelOrd[idx] = ord;
    405397            eLow[idx]  = tmin;
    406398            eHigh[idx] = tmax;
     399            isUsed[ii] = 1;
    407400          }
    408401        }
    … …  
    416409    }
    417410    eLow[0] = 0.0;
    418     if(n >= nmax) eLow.resize(nmax+1);
    419411    eLow[n] = eHigh[n-1];
    420412
    … …  
    430422  }
    431423
     424  currRegionModel = setOfRegionModels[0];
     425
    432426  // Access to materials and build cuts
    433 
    434427  for(G4int i=0; i<numOfCouples; i++) {
    435428
    … …  
    439432    const G4ProductionCuts* pcuts = couple->GetProductionCuts();
    440433 
    441     G4int reg = nRegions;
    442     do {reg--;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));
    443     idxOfRegionModels[i] = reg;
    444 
     434    G4int reg = 0;
     435    if(nRegions > 1) {
     436      reg = nRegions;
     437      do {reg--;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));
     438      idxOfRegionModels[i] = reg;
     439    }
    445440    if(1 < verboseLevel) {
    446441      G4cout << "G4EmModelManager::Initialise() for "
    447              << material->GetName()
    448              << " indexOfCouple= " << i
    449              << " indexOfRegion= " << reg
    450              << G4endl;
     442             << material->GetName()
     443             << " indexOfCouple= " << i
     444             << " indexOfRegion= " << reg
     445             << G4endl;
    451446    }
    452447
    … …  
    494489
    495490  for(G4int jj=0; jj<nEmModels; jj++) {
    496     models[jj]->Initialise(particle, theCuts);
    497     if(flucModels[jj]) flucModels[jj]->InitialiseMe(particle);
    498   }
    499 
     491    if(1 == isUsed[jj]) {
     492      models[jj]->Initialise(particle, theCuts);
     493      if(flucModels[jj]) flucModels[jj]->InitialiseMe(particle);
     494    }
     495  }
    500496
    501497  if(1 < verboseLevel) {
    … …  
    534530  }
    535531
    536   G4int reg  = idxOfRegionModels[i];
     532  G4int reg  = 0;
     533  if(nRegions > 1) reg = idxOfRegionModels[i];
    537534  const G4RegionModels* regModels = setOfRegionModels[reg];
    538535  G4int nmod = regModels->NumberOfModels();
    … …  
    683680  }
    684681
    685   G4int reg  = idxOfRegionModels[i];
     682  G4int reg  = 0;
     683  if(nRegions > 1) reg  = idxOfRegionModels[i];
    686684  const G4RegionModels* regModels = setOfRegionModels[reg];
    687685  G4int nmod = regModels->NumberOfModels();
    … …  
    799797    G4RegionModels* r = setOfRegionModels[i];
    800798    const G4Region* reg = r->Region();
    801     if(verb > 1 || nRegions > 1) {
    802     }
     799    //    if(verb > -1 || nRegions > 1) {
     800    // }
    803801    G4int n = r->NumberOfModels(); 
    804     if(verb > 1 || n > 0) {
     802    if(verb > -1 || n > 0) {
    805803      G4cout << "      ===== EM models for the G4Region  " << reg->GetName()
    806804             << " ======" << G4endl;;

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmProcessOptions.cc,v 1.24 2008/04/17 10:33:27 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4EmProcessOptions.cc,v 1.26 2009/02/18 14:43:27 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5959#include "G4VMultipleScattering.hh"
    6060#include "G4Region.hh"
     61#include "G4RegionStore.hh"
    6162
    6263//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    392393//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    393394
    394 void G4EmProcessOptions::ActivateDeexcitation(G4bool val, const G4Region* r)
    395 {
    396   const std::vector<G4VEnergyLossProcess*>& v =
    397         theManager->GetEnergyLossProcessVector();
    398   std::vector<G4VEnergyLossProcess*>::const_iterator itr;
    399   for(itr = v.begin(); itr != v.end(); itr++) {
    400     G4VEnergyLossProcess* p = *itr;
    401     if(p) p->ActivateDeexcitation(val,r);
    402   }
    403   const std::vector<G4VEmProcess*>& w =
    404         theManager->GetEmProcessVector();
    405   std::vector<G4VEmProcess*>::const_iterator itp;
    406   for(itp = w.begin(); itp != w.end(); itp++) {
    407     G4VEmProcess* q = *itp;
    408     if(q) q->ActivateDeexcitation(val,r);
     395void 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    }
    409429  }
    410430}

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

    r1007 r1055  
    2525//
    2626//
    27 // $Id: G4EnergyLossMessenger.cc,v 1.35 2008/10/20 13:27:45 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4EnergyLossMessenger.cc,v 1.37 2009/02/18 14:43:27 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    178178  aplCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
    179179
     180  deexCmd = new G4UIcommand("/process/em/deexcitation",this);
     181  deexCmd->SetGuidance("Set deexcitation flag per process and G4Region.");
     182  deexCmd->SetGuidance("  procName  : process name");
     183  deexCmd->SetGuidance("  flag      : flag");
     184  deexCmd->SetGuidance("  regName   : G4Region name");
     185
     186  G4UIparameter* pName = new G4UIparameter("pName",'s',false);
     187  deexCmd->SetParameter(pName);
     188
     189  G4UIparameter* flag = new G4UIparameter("flag",'s',false);
     190  deexCmd->SetParameter(flag);
     191
     192  G4UIparameter* regName = new G4UIparameter("regName",'s',false);
     193  deexCmd->SetParameter(regName);
     194
    180195  dedxCmd = new G4UIcmdWithAnInteger("/process/eLoss/binsDEDX",this);
    181196  dedxCmd->SetGuidance("Set number of bins for DEDX tables");
    … …  
    259274  delete MinSubSecCmd;
    260275  delete StepFuncCmd;
     276  delete deexCmd;
    261277  delete eLossDirectory;
    262278  delete mscDirectory;
    … …  
    320336  }
    321337
     338  if (command == deexCmd) {
     339    G4String s1 (""), s2(""), s3("");
     340    G4bool b = false;
     341    std::istringstream is(newValue);
     342    is >> s1 >> s2 >> s3;
     343    if(s2 == "true") b = true;
     344    opt->ActivateDeexcitation(s1,b,s3);
     345  }
     346
    322347  if (command == mscCmd) {
    323348    if(newValue == "Minimal")

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LossTableBuilder.cc,v 1.27 2008/07/22 15:55:15 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4LossTableBuilder.cc,v 1.28 2009/02/18 16:24:47 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    132132      G4double dedx1  = pv->GetValue(elow, b);
    133133
     134      //G4cout << "nbins= " << nbins << " dedx1= " << dedx1 << G4endl;
     135
    134136      // protection for specific cases dedx=0
    135137      if(dedx1 == 0.0) {
    136         for (size_t k=1; k<nbins; k++) {
     138        for (size_t k=1; k<nbins-1; k++) {
    137139          bin0++;
    138140          elow  = pv->GetLowEdgeEnergy(k);
    … …  
    143145      }
    144146
     147      //G4cout << "nbins= " << nbins << " elow= " << elow << " ehigh= " << ehigh << G4endl;
    145148      // initialisation of a new vector
    146149      G4PhysicsLogVector* v = new G4PhysicsLogVector(elow, ehigh, nbins);
     150      // dedx is exect zero
     151      if(nbins <= 2) {
     152        v->PutValue(0,1000.);
     153        v->PutValue(1,2000.);
     154        G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
     155        return;
     156      }
    147157      v->SetSpline(splineFlag);
    148158

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LossTableManager.cc,v 1.95 2008/11/13 18:23:39 schaelic Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4LossTableManager.cc,v 1.96 2009/04/09 16:10:57 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9393#include "G4EmCorrections.hh"
    9494#include "G4EmSaturation.hh"
     95#include "G4EmConfigurator.hh"
    9596#include "G4EmTableType.hh"
    9697#include "G4LossTableBuilder.hh"
    … …  
    164165  emCorrections= new G4EmCorrections();
    165166  emSaturation = new G4EmSaturation();
     167  emConfigurator = new G4EmConfigurator();
    166168  integral = true;
    167169  integralActive = false;
    … …  
    932934}
    933935
     936//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     937
     938G4EmConfigurator* G4LossTableManager::EmConfigurator()
     939{
     940  return emConfigurator;
     941}
     942
    934943//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmFluctuationModel.cc,v 1.3 2008/07/15 16:56:39 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VEmFluctuationModel.cc,v 1.4 2009/02/19 11:25:50 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6565}
    6666
     67void G4VEmFluctuationModel::InitialiseMe(const G4ParticleDefinition*)
     68{}
    6769
     70void G4VEmFluctuationModel::SetParticleAndCharge(const G4ParticleDefinition*,
     71                                                 G4double)
     72{}
     73
     74

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmModel.cc,v 1.20 2008/11/13 23:13:18 schaelic Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VEmModel.cc,v 1.27 2009/05/26 15:00:49 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4141// 25.10.2005 Set default highLimit=100.TeV (V.Ivanchenko)
    4242// 06.02.2006 add method ComputeMeanFreePath() (mma)
     43// 16.02.2009 Move implementations of virtual methods to source (VI)
    4344//
    4445//
    … …  
    5354#include "G4LossTableManager.hh"
    5455#include "G4ProductionCutsTable.hh"
     56#include "G4ParticleChangeForLoss.hh"
     57#include "G4ParticleChangeForGamma.hh"
    5558
    5659//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    6063  fluc(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV),
    6164  polarAngleLimit(0.0),secondaryThreshold(DBL_MAX),theLPMflag(false),
    62   pParticleChange(0),nuclearStopping(false),nsec(5)
     65  pParticleChange(0),nuclearStopping(false),
     66  currentCouple(0),currentElement(0),
     67  nsec(5),flagDeexcitation(false)
    6368{
    6469  xsec.resize(nsec);
    … …  
    7883    }
    7984  }
     85}
     86
     87//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     88
     89G4ParticleChangeForLoss* G4VEmModel::GetParticleChangeForLoss()
     90{
     91  G4ParticleChangeForLoss* p = 0;
     92  if (pParticleChange) {
     93    p = static_cast<G4ParticleChangeForLoss*>(pParticleChange);
     94  } else {
     95    p = new G4ParticleChangeForLoss();
     96  }
     97  return p;
     98}
     99
     100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     101
     102G4ParticleChangeForGamma* G4VEmModel::GetParticleChangeForGamma()
     103{
     104  G4ParticleChangeForGamma* p = 0;
     105  if (pParticleChange) {
     106    p = static_cast<G4ParticleChangeForGamma*>(pParticleChange);
     107  } else {
     108    p = new G4ParticleChangeForGamma();
     109  }
     110  return p;
     111}
     112
     113//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     114
     115void G4VEmModel::InitialiseElementSelectors(const G4ParticleDefinition* p,
     116                                            const G4DataVector& cuts)
     117{
     118  // initialise before run
     119  flagDeexcitation = false;
     120
     121  G4int nbins = G4int(std::log10(highLimit/lowLimit) + 0.5);
     122  if(nbins < 3) nbins = 3;
     123  G4bool spline = G4LossTableManager::Instance()->SplineFlag();
     124
     125  G4ProductionCutsTable* theCoupleTable=
     126    G4ProductionCutsTable::GetProductionCutsTable();
     127  G4int numOfCouples = theCoupleTable->GetTableSize();
     128
     129  // prepare vector
     130  if(numOfCouples > nSelectors) elmSelectors.reserve(numOfCouples);
     131
     132  // initialise vector
     133  for(G4int i=0; i<numOfCouples; i++) {
     134    const G4MaterialCutsCouple* couple =
     135      theCoupleTable->GetMaterialCutsCouple(i);
     136    const G4Material* material = couple->GetMaterial();
     137    G4int idx = couple->GetIndex();
     138
     139    // selector already exist check if should be deleted
     140    G4bool create = true;
     141    if(i < nSelectors) {
     142      if(elmSelectors[i]) {
     143        if(material == elmSelectors[i]->GetMaterial()) create = false;
     144        else delete elmSelectors[i];
     145      }
     146    } else {
     147      nSelectors++;
     148      elmSelectors.push_back(0);
     149    }
     150    if(create) {
     151      elmSelectors[i] = new G4EmElementSelector(this,material,nbins,
     152                                                lowLimit,highLimit,spline);
     153    }
     154    elmSelectors[i]->Initialise(p, cuts[idx]);
     155    //elmSelectors[i]->Dump(p);
     156  }
     157}
     158
     159//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     160
     161G4double G4VEmModel::ComputeDEDXPerVolume(const G4Material*,
     162                                          const G4ParticleDefinition*,
     163                                          G4double,G4double)
     164{
     165  return 0.0;
    80166}
    81167
    … …  
    107193//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    108194
    109 G4double G4VEmModel::ComputeMeanFreePath(const G4ParticleDefinition* p,
    110                                          G4double ekin,
    111                                          const G4Material* material,     
    112                                          G4double emin,
    113                                          G4double emax)
    114 {
    115   G4double mfp = DBL_MAX;
    116   G4double cross = CrossSectionPerVolume(material,p,ekin,emin,emax);
    117   if (cross > DBL_MIN) mfp = 1./cross;
    118   return mfp;
    119 }
    120 
    121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    122 
    123 void G4VEmModel::InitialiseElementSelectors(const G4ParticleDefinition* p,
    124                                             const G4DataVector& cuts)
    125 {
    126   G4int nbins = G4int(std::log10(highLimit/lowLimit) + 0.5);
    127   if(nbins < 3) nbins = 3;
    128   G4bool spline = G4LossTableManager::Instance()->SplineFlag();
    129 
    130   G4ProductionCutsTable* theCoupleTable=
    131     G4ProductionCutsTable::GetProductionCutsTable();
    132   G4int numOfCouples = theCoupleTable->GetTableSize();
    133 
    134   // prepare vector
    135   if(numOfCouples > nSelectors) {
    136     elmSelectors.resize(numOfCouples);
    137     nSelectors = numOfCouples;
    138   }
    139 
    140   // initialise vector
    141   for(G4int i=0; i<numOfCouples; i++) {
    142     const G4MaterialCutsCouple* couple =
    143       theCoupleTable->GetMaterialCutsCouple(i);
    144     const G4Material* material = couple->GetMaterial();
    145     G4int idx = couple->GetIndex();
    146 
    147     // selector already exist check if should be deleted
    148     G4bool create = true;
    149     if(elmSelectors[i]) {
    150       if(material == elmSelectors[i]->GetMaterial()) create = false;
    151       else delete elmSelectors[i];
    152     }
    153     if(create) {
    154       elmSelectors[i] = new G4EmElementSelector(this,material,nbins,
    155                                                 lowLimit,highLimit,spline);
    156     }
    157     elmSelectors[i]->Initialise(p, cuts[idx]);
    158     //elmSelectors[i]->Dump(p);
    159   }
    160 }
    161 
    162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    163 
    164 
     195G4double G4VEmModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
     196                                                G4double, G4double, G4double,
     197                                                G4double, G4double)
     198{
     199  return 0.0;
     200}
     201
     202//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     203
     204G4double G4VEmModel::MinEnergyCut(const G4ParticleDefinition*,
     205                                  const G4MaterialCutsCouple*)
     206{
     207  return 0.0;
     208}
     209
     210//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     211
     212G4double G4VEmModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
     213                                          const G4Material*, G4double)
     214{
     215  G4double q = p->GetPDGCharge()/CLHEP::eplus;
     216  return q*q;
     217}
     218
     219//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     220
     221G4double G4VEmModel::GetParticleCharge(const G4ParticleDefinition* p,
     222                                       const G4Material*, G4double)
     223{
     224  return p->GetPDGCharge();
     225}
     226
     227//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     228
     229void G4VEmModel::CorrectionsAlongStep(const G4MaterialCutsCouple*,
     230                                      const G4DynamicParticle*,
     231                                      G4double&,G4double&,G4double)
     232{}
     233
     234//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     235
     236void G4VEmModel::SampleDeexcitationAlongStep(const G4Material*,
     237                                             const G4Track&,
     238                                             G4double& )
     239{}
     240
     241//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     242
     243G4double G4VEmModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
     244                                        G4double kineticEnergy)
     245{
     246  return kineticEnergy;
     247}
     248
     249//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     250
     251void G4VEmModel::DefineForRegion(const G4Region*)
     252{}
     253
     254//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     255
     256void G4VEmModel::SetupForMaterial(const G4ParticleDefinition*,
     257                                  const G4Material*, G4double)
     258{}
     259
     260//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmProcess.cc,v 1.60 2008/10/17 14:46:16 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VEmProcess.cc,v 1.66 2009/04/17 10:35:32 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    7878#include "G4Positron.hh"
    7979#include "G4PhysicsTableHelper.hh"
     80#include "G4EmConfigurator.hh"
    8081
    8182//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    9192  applyCuts(false),
    9293  startFromNull(true),
    93   nRegions(0),
    94   selectedModel(0),
     94  useDeexcitation(false),
     95  nDERegions(0),
     96  idxDERegions(0),
     97  currentModel(0),
    9598  particle(0),
    9699  currentCouple(0)
    … …  
    100103  // Size of tables assuming spline
    101104  minKinEnergy = 0.1*keV;
    102   maxKinEnergy = 100.0*TeV;
    103   nLambdaBins  = 84;
     105  maxKinEnergy = 10.0*TeV;
     106  nLambdaBins  = 77;
    104107
    105108  // default lambda factor
    … …  
    139142//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    140143
     144void G4VEmProcess::Clear()
     145{
     146  delete [] theEnergyOfCrossSectionMax;
     147  delete [] theCrossSectionMax;
     148  delete [] idxDERegions;
     149  theEnergyOfCrossSectionMax = 0;
     150  theCrossSectionMax = 0;
     151  idxDERegions = 0;
     152  currentCouple = 0;
     153  preStepLambda = 0.0;
     154  mfpKinEnergy  = DBL_MAX;
     155  deRegions.clear();
     156  nDERegions = 0;
     157}
     158
     159//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     160
     161void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* p,
     162                              const G4Region* region)
     163{
     164  G4VEmFluctuationModel* fm = 0;
     165  modelManager->AddEmModel(order, p, fm, region);
     166  if(p) p->SetParticleChange(pParticleChange);
     167}
     168
     169//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     170
     171void G4VEmProcess::SetModel(G4VEmModel* p, G4int index)
     172{
     173  G4int n = emModels.size();
     174  if(index >= n) for(G4int i=n; i<index+1; i++) {emModels.push_back(0);}
     175  emModels[index] = p;
     176}
     177
     178//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     179
     180G4VEmModel* G4VEmProcess::Model(G4int index)
     181{
     182  G4VEmModel* p = 0;
     183  if(index >= 0 && index <  G4int(emModels.size())) p = emModels[index];
     184  return p;
     185}
     186
     187//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     188
     189void G4VEmProcess::UpdateEmModel(const G4String& nam,
     190                                 G4double emin, G4double emax)
     191{
     192  modelManager->UpdateEmModel(nam, emin, emax);
     193}
     194
     195//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     196
     197G4VEmModel* G4VEmProcess::GetModelByIndex(G4int idx, G4bool ver)
     198{
     199  return modelManager->GetModel(idx, ver);
     200}
     201
     202//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     203
    141204void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
    142205{
    … …  
    149212           << G4endl;
    150213  }
     214
     215  (G4LossTableManager::Instance())->EmConfigurator()->AddModels();
    151216
    152217  if(particle == &part) {
    … …  
    159224    theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
    160225    theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
    161     if(buildLambdaTable)
     226    if(buildLambdaTable){
    162227      theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
    163   }
    164 }
    165 
    166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    167 
    168 void G4VEmProcess::Clear()
    169 {
    170   if(theEnergyOfCrossSectionMax) delete [] theEnergyOfCrossSectionMax;
    171   if(theCrossSectionMax) delete [] theCrossSectionMax;
    172   theEnergyOfCrossSectionMax = 0;
    173   theCrossSectionMax = 0;
    174   currentCouple = 0;
    175   preStepLambda = 0.0;
    176   mfpKinEnergy  = DBL_MAX;
     228    }
     229  }
     230  // Sub Cutoff and Deexcitation
     231  if (nDERegions>0) {
     232
     233    const G4ProductionCutsTable* theCoupleTable=
     234          G4ProductionCutsTable::GetProductionCutsTable();
     235    size_t numOfCouples = theCoupleTable->GetTableSize();
     236
     237    idxDERegions = new G4bool[numOfCouples];
     238
     239    for (size_t j=0; j<numOfCouples; j++) {
     240
     241      const G4MaterialCutsCouple* couple =
     242        theCoupleTable->GetMaterialCutsCouple(j);
     243      const G4ProductionCuts* pcuts = couple->GetProductionCuts();
     244      G4bool reg = false;
     245      for(G4int i=0; i<nDERegions; i++) {
     246        if(deRegions[i]) {
     247          if(pcuts == deRegions[i]->GetProductionCuts()) reg = true;
     248        }
     249      }
     250      idxDERegions[j] = reg;
     251    }
     252  }
     253  if (1 < verboseLevel && nDERegions>0) {
     254    G4cout << " Deexcitation is activated for regions: " << G4endl;
     255    for (G4int i=0; i<nDERegions; i++) {
     256      const G4Region* r = deRegions[i];
     257      G4cout << "           " << r->GetName() << G4endl;
     258    }
     259  }
    177260}
    178261
    … …  
    237320      G4cout << *theLambdaTable << G4endl;
    238321    }
     322  }
     323}
     324
     325//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     326
     327void G4VEmProcess::PrintInfoDefinition()
     328{
     329  if(verboseLevel > 0) {
     330    G4cout << G4endl << GetProcessName() << ":   for  "
     331           << particle->GetParticleName();
     332    if(integral) G4cout << ", integral: 1 ";
     333    if(applyCuts) G4cout << ", applyCuts: 1 ";
     334    G4cout << "    SubType= " << GetProcessSubType() << G4endl;
     335    if(buildLambdaTable) {
     336      G4cout << "      Lambda tables from "
     337             << G4BestUnit(minKinEnergy,"Energy")
     338             << " to "
     339             << G4BestUnit(maxKinEnergy,"Energy")
     340             << " in " << nLambdaBins << " bins, spline: "
     341             << (G4LossTableManager::Instance())->SplineFlag()
     342             << G4endl;
     343    }
     344    PrintInfo();
     345    modelManager->DumpModelList(verboseLevel);
     346  }
     347
     348  if(verboseLevel > 2 && buildLambdaTable) {
     349    G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
     350    if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
    239351  }
    240352}
    … …  
    304416//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    305417
    306 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
    307                                        G4double,
    308                                        G4ForceCondition* condition)
    309 {
    310   *condition = NotForced;
    311   return G4VEmProcess::MeanFreePath(track);
    312 }
    313 
    314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    315 
    316418G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
    317419                                              const G4Step&)
    … …  
    342444  }
    343445
    344   G4VEmModel* currentModel = SelectModel(finalT);
    345 
     446  SelectModel(finalT);
     447  if(useDeexcitation) {
     448    currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
     449  }
    346450  /* 
    347451  if(0 < verboseLevel) {
    … …  
    404508//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    405509
    406 void G4VEmProcess::PrintInfoDefinition()
    407 {
    408   if(verboseLevel > 0) {
    409     G4cout << G4endl << GetProcessName() << ":   for  "
    410            << particle->GetParticleName();
    411     if(integral) G4cout << ", integral: 1 ";
    412     if(applyCuts) G4cout << ", applyCuts: 1 ";
    413     G4cout << "    SubType= " << GetProcessSubType() << G4endl;
    414     if(buildLambdaTable) {
    415       G4cout << "      Lambda tables from "
    416              << G4BestUnit(minKinEnergy,"Energy")
    417              << " to "
    418              << G4BestUnit(maxKinEnergy,"Energy")
    419              << " in " << nLambdaBins << " bins, spline: "
    420              << (G4LossTableManager::Instance())->SplineFlag()
    421              << G4endl;
    422     }
    423     PrintInfo();
    424     modelManager->DumpModelList(verboseLevel);
    425   }
    426 
    427   if(verboseLevel > 2 && buildLambdaTable) {
    428     G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
    429     if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
    430   }
    431 }
    432 
    433 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    434 
    435 G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
    436                                              const G4MaterialCutsCouple* couple)
    437 {
    438   // Cross section per atom is calculated
    439   DefineMaterial(couple);
    440   G4double cross = 0.0;
    441   G4bool b;
    442   if(theLambdaTable) {
    443     cross = (((*theLambdaTable)[currentMaterialIndex])->
    444                            GetValue(kineticEnergy, b));
    445   } else {
    446     G4VEmModel* model = SelectModel(kineticEnergy);
    447     cross =
    448       model->CrossSectionPerVolume(currentMaterial,particle,kineticEnergy);
    449   }
    450 
    451   return cross;
    452 }
    453 
    454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    455 
    456510G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
    457511                                       const G4String& directory,
    … …  
    521575
    522576  return yes;
     577}
     578
     579//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     580
     581void G4VEmProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
     582{
     583  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     584  const G4Region* reg = r;
     585  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
     586
     587  // the region is in the list
     588  if (nDERegions) {
     589    for (G4int i=0; i<nDERegions; i++) {
     590      if (reg == deRegions[i]) {
     591        if(!val) deRegions[i] = 0;
     592        return;
     593      }
     594    }
     595  }
     596
     597  // new region
     598  if(val) {
     599    useDeexcitation = true;
     600    deRegions.push_back(reg);
     601    nDERegions++;
     602  } else {
     603    useDeexcitation = false;
     604  }
     605}
     606
     607//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     608
     609G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
     610                                             const G4MaterialCutsCouple* couple)
     611{
     612  // Cross section per atom is calculated
     613  DefineMaterial(couple);
     614  G4double cross = 0.0;
     615  G4bool b;
     616  if(theLambdaTable) {
     617    cross = (((*theLambdaTable)[currentMaterialIndex])->
     618                           GetValue(kineticEnergy, b));
     619  } else {
     620    SelectModel(kineticEnergy);
     621    cross = currentModel->CrossSectionPerVolume(currentMaterial,
     622                                                particle,kineticEnergy);
     623  }
     624
     625  return cross;
     626}
     627
     628//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     629
     630G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
     631                                       G4double,
     632                                       G4ForceCondition* condition)
     633{
     634  *condition = NotForced;
     635  return G4VEmProcess::MeanFreePath(track);
    523636}
    524637

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEnergyLossProcess.cc,v 1.143 2008/10/17 14:46:16 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VEnergyLossProcess.cc,v 1.149 2009/04/17 10:35:32 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    142142#include "G4SafetyHelper.hh"
    143143#include "G4TransportationManager.hh"
     144#include "G4EmConfigurator.hh"
    144145
    145146//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    150151  secondaryParticle(0),
    151152  nSCoffRegions(0),
     153  nDERegions(0),
    152154  idxSCoffRegions(0),
     155  idxDERegions(0),
    153156  nProcesses(0),
    154157  theDEDXTable(0),
    … …  
    176179  isIonisation(true),
    177180  useSubCutoff(false),
     181  useDeexcitation(false),
    178182  particle(0),
    179183  currentCouple(0),
    … …  
    188192  // Size of tables assuming spline
    189193  minKinEnergy     = 0.1*keV;
    190   maxKinEnergy     = 100.0*TeV;
    191   nBins            = 84;
     194  maxKinEnergy     = 10.0*TeV;
     195  nBins            = 77;
    192196  maxKinEnergyCSDA = 1.0*GeV;
    193197  nBinsCSDA        = 35;
    … …  
    237241           << G4endl;
    238242  delete vstrag;
    239   Clear();
     243  Clean();
    240244
    241245  if ( !baseParticle ) {
    … …  
    291295//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    292296
    293 void G4VEnergyLossProcess::Clear()
     297void G4VEnergyLossProcess::Clean()
    294298{
    295299  if(1 < verboseLevel) {
    … …  
    302306  delete [] theCrossSectionMax;
    303307  delete [] idxSCoffRegions;
     308  delete [] idxDERegions;
    304309
    305310  theDEDXAtMaxEnergy = 0;
    … …  
    312317  scProcesses.clear();
    313318  nProcesses = 0;
     319}
     320
     321//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     322
     323G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
     324                                                const G4Material*,
     325                                                G4double cut)
     326{
     327  return cut;
     328}
     329
     330//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     331
     332void G4VEnergyLossProcess::AddEmModel(G4int order, G4VEmModel* p,
     333                                      G4VEmFluctuationModel* fluc,
     334                                      const G4Region* region)
     335{
     336  modelManager->AddEmModel(order, p, fluc, region);
     337  if(p) p->SetParticleChange(pParticleChange, fluc);
     338}
     339
     340//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     341
     342void G4VEnergyLossProcess::UpdateEmModel(const G4String& nam,
     343                                         G4double emin, G4double emax)
     344{
     345  modelManager->UpdateEmModel(nam, emin, emax);
     346}
     347
     348//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     349
     350void G4VEnergyLossProcess::SetEmModel(G4VEmModel* p, G4int index)
     351{
     352  G4int n = emModels.size();
     353  if(index >= n) for(G4int i=n; i<index+1; i++) {emModels.push_back(0);}
     354  emModels[index] = p;
     355}
     356
     357//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     358
     359G4VEmModel* G4VEnergyLossProcess::EmModel(G4int index)
     360{
     361  G4VEmModel* p = 0;
     362  if(index >= 0 && index <  G4int(emModels.size())) p = emModels[index];
     363  return p;
     364}
     365
     366//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     367
     368G4VEmModel* G4VEnergyLossProcess::GetModelByIndex(G4int idx, G4bool ver)
     369{
     370  return modelManager->GetModel(idx, ver);
     371}
     372
     373//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     374
     375G4int G4VEnergyLossProcess::NumberOfModels()
     376{
     377  return modelManager->NumberOfModels();
    314378}
    315379
    … …  
    364428  }
    365429
    366   Clear();
     430  Clean();
     431  lManager->EmConfigurator()->AddModels();
    367432
    368433  // Base particle and set of models can be defined here
    … …  
    407472                                     minSubRange, verboseLevel);
    408473
    409   // Sub Cutoff Regime
    410   if (nSCoffRegions>0) {
     474  // Sub Cutoff and Deexcitation
     475  if (nSCoffRegions>0 || nDERegions>0) {
    411476    theSubCuts = modelManager->SubCutoff();
    412477
    … …  
    414479          G4ProductionCutsTable::GetProductionCutsTable();
    415480    size_t numOfCouples = theCoupleTable->GetTableSize();
    416     idxSCoffRegions = new G4int[numOfCouples];
     481
     482    if(nSCoffRegions>0) idxSCoffRegions = new G4bool[numOfCouples];
     483    if(nDERegions>0) idxDERegions = new G4bool[numOfCouples];
    417484 
    418485    for (size_t j=0; j<numOfCouples; j++) {
    … …  
    421488        theCoupleTable->GetMaterialCutsCouple(j);
    422489      const G4ProductionCuts* pcuts = couple->GetProductionCuts();
    423       G4int reg = 0;
    424       for(G4int i=0; i<nSCoffRegions; i++) {
    425         if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = 1;
    426       }
    427       idxSCoffRegions[j] = reg;
     490     
     491      if(nSCoffRegions>0) {
     492        G4bool reg = false;
     493        for(G4int i=0; i<nSCoffRegions; i++) {
     494          if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = true;
     495        }
     496        idxSCoffRegions[j] = reg;
     497      }
     498      if(nDERegions>0) {
     499        G4bool reg = false;
     500        for(G4int i=0; i<nDERegions; i++) {
     501          if( pcuts == deRegions[i]->GetProductionCuts()) reg = true;
     502        }
     503        idxDERegions[j] = reg;
     504      }
    428505    }
    429506  }
    … …  
    445522      }
    446523    }
     524    if (nDERegions) {
     525      G4cout << " Deexcitation is ON for regions: " << G4endl;
     526      for (G4int i=0; i<nDERegions; i++) {
     527        const G4Region* r = deRegions[i];
     528        G4cout << "           " << r->GetName() << G4endl;
     529      }
     530    }
    447531  }
    448532}
    … …  
    479563    if(isIonisation) G4cout << "  isIonisation  flag = 1";
    480564    G4cout << G4endl;
    481   }
    482 }
    483 
    484 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    485 
    486 void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
    487 {
    488   G4RegionStore* regionStore = G4RegionStore::GetInstance();
    489   if(val) {
    490     useSubCutoff = true;
    491     if (!r) {r = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
    492     if (nSCoffRegions) {
    493       for (G4int i=0; i<nSCoffRegions; i++) {
    494         if (r == scoffRegions[i]) return;
    495       }
    496     }
    497     scoffRegions.push_back(r);
    498     nSCoffRegions++;
    499   } else {
    500     useSubCutoff = false;
    501565  }
    502566}
    … …  
    640704//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    641705
    642 G4double G4VEnergyLossProcess::GetContinuousStepLimit(
    643                 const G4Track&,
    644                 G4double, G4double, G4double&)
    645 {
    646   return DBL_MAX;
     706void G4VEnergyLossProcess::PrintInfoDefinition()
     707{
     708  if(0 < verboseLevel) {
     709    G4cout << G4endl << GetProcessName() << ":   for  "
     710           << particle->GetParticleName()
     711           << "    SubType= " << GetProcessSubType()
     712           << G4endl
     713           << "      dE/dx and range tables from "
     714           << G4BestUnit(minKinEnergy,"Energy")
     715           << " to " << G4BestUnit(maxKinEnergy,"Energy")
     716           << " in " << nBins << " bins" << G4endl
     717           << "      Lambda tables from threshold to "
     718           << G4BestUnit(maxKinEnergy,"Energy")
     719           << " in " << nBins << " bins, spline: "
     720           << (G4LossTableManager::Instance())->SplineFlag()
     721           << G4endl;
     722    if(theRangeTableForLoss && isIonisation) {
     723      G4cout << "      finalRange(mm)= " << finalRange/mm
     724             << ", dRoverRange= " << dRoverRange
     725             << ", integral: " << integral
     726             << ", fluct: " << lossFluctuationFlag
     727             << ", linLossLimit= " << linLossLimit
     728             << G4endl;
     729    }
     730    PrintInfo();
     731    modelManager->DumpModelList(verboseLevel);
     732    if(theCSDARangeTable && isIonisation) {
     733      G4cout << "      CSDA range table up"
     734             << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
     735             << " in " << nBinsCSDA << " bins" << G4endl;
     736    }
     737    if(nSCoffRegions>0 && isIonisation) {
     738      G4cout << "      Subcutoff sampling in " << nSCoffRegions
     739             << " regions" << G4endl;
     740    }
     741    if(2 < verboseLevel) {
     742      G4cout << "      DEDXTable address= " << theDEDXTable << G4endl;
     743      if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
     744      G4cout << "non restricted DEDXTable address= "
     745             << theDEDXunRestrictedTable << G4endl;
     746      if(theDEDXunRestrictedTable && isIonisation) {
     747           G4cout << (*theDEDXunRestrictedTable) << G4endl;
     748      }
     749      if(theDEDXSubTable && isIonisation) {
     750        G4cout << (*theDEDXSubTable) << G4endl;
     751      }
     752      G4cout << "      CSDARangeTable address= " << theCSDARangeTable
     753             << G4endl;
     754      if(theCSDARangeTable && isIonisation) {
     755        G4cout << (*theCSDARangeTable) << G4endl;
     756      }
     757      G4cout << "      RangeTableForLoss address= " << theRangeTableForLoss
     758             << G4endl;
     759      if(theRangeTableForLoss && isIonisation) {
     760             G4cout << (*theRangeTableForLoss) << G4endl;
     761      }
     762      G4cout << "      InverseRangeTable address= " << theInverseRangeTable
     763             << G4endl;
     764      if(theInverseRangeTable && isIonisation) {
     765             G4cout << (*theInverseRangeTable) << G4endl;
     766      }
     767      G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
     768      if(theLambdaTable && isIonisation) {
     769        G4cout << (*theLambdaTable) << G4endl;
     770      }
     771      G4cout << "      SubLambdaTable address= " << theSubLambdaTable << G4endl;
     772      if(theSubLambdaTable && isIonisation) {
     773        G4cout << (*theSubLambdaTable) << G4endl;
     774      }
     775    }
     776  }
     777}
     778
     779//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     780
     781void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
     782{
     783  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     784  const G4Region* reg = r;
     785  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
     786
     787  // the region is in the list
     788  if (nSCoffRegions) {
     789    for (G4int i=0; i<nSCoffRegions; i++) {
     790      if (reg == scoffRegions[i]) {
     791        if(!val) deRegions[i] = 0;
     792        return;
     793      }
     794    }
     795  }
     796
     797  // new region
     798  if(val) {
     799    useSubCutoff = true;
     800    scoffRegions.push_back(reg);
     801    nSCoffRegions++;
     802  } else {
     803    useSubCutoff = false;
     804  }
     805}
     806
     807//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     808
     809void G4VEnergyLossProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
     810{
     811  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     812  const G4Region* reg = r;
     813  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
     814
     815  // the region is in the list
     816  if (nDERegions) {
     817    for (G4int i=0; i<nDERegions; i++) {
     818      if (reg == deRegions[i]) {
     819        if(!val) deRegions[i] = 0;
     820        return;
     821      }
     822    }
     823  }
     824
     825  // new region
     826  if(val) {
     827    useDeexcitation = true;
     828    deRegions.push_back(reg);
     829    nDERegions++;
     830  } else {
     831    useDeexcitation = false;
     832  }
    647833}
    648834
    … …  
    674860  //  <<" stepLimit= "<<x<<G4endl;
    675861  return x;
    676 }
    677 
    678 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    679 
    680 G4double G4VEnergyLossProcess::GetMeanFreePath(
    681                              const G4Track& track,
    682                              G4double,
    683                              G4ForceCondition* condition)
    684 
    685 {
    686   *condition = NotForced;
    687   return MeanFreePath(track);
    688862}
    689863
    … …  
    806980  if (length >= fRange) {
    807981    eloss = preStepKinEnergy;
    808     currentModel->CorrectionsAlongStep(currentCouple, dynParticle,
    809                                        eloss, esecdep, length);
     982    if (useDeexcitation) {
     983      if(idxDERegions[currentMaterialIndex]) {
     984        currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
     985        if(eloss < 0.0) eloss = 0.0;
     986      }
     987    }
    810988    fParticleChange.SetProposedKineticEnergy(0.0);
    811     fParticleChange.ProposeLocalEnergyDeposit(preStepKinEnergy);
     989    fParticleChange.ProposeLocalEnergyDeposit(eloss);
    812990    return &fParticleChange;
    813991  }
    … …  
    9341112      G4double tmax =
    9351113        std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut);
     1114      G4double emean = eloss;
    9361115      eloss = fluc->SampleFluctuations(currentMaterial,dynParticle,
    937                                        tmax,length,eloss);
     1116                                       tmax,length,emean);
    9381117      /*                           
    9391118      if(-1 < verboseLevel)
    … …  
    9501129  eloss += esecdep;
    9511130  if(eloss < 0.0) eloss = 0.0;
     1131
     1132  // deexcitation
     1133  else if (useDeexcitation) {
     1134    if(idxDERegions[currentMaterialIndex]) {
     1135      currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
     1136      if(eloss < 0.0) eloss = 0.0;
     1137    }
     1138  }
    9521139
    9531140  // Energy balanse
    … …  
    11061293
    11071294  SelectModel(postStepScaledEnergy);
     1295  if(useDeexcitation) {
     1296    currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
     1297  }
     1298
    11081299  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
    11091300  G4double tcut = (*theCuts)[currentMaterialIndex];
    … …  
    11401331//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    11411332
    1142 void G4VEnergyLossProcess::PrintInfoDefinition()
    1143 {
    1144   if(0 < verboseLevel) {
    1145     G4cout << G4endl << GetProcessName() << ":   for  "
    1146            << particle->GetParticleName()
    1147            << "    SubType= " << GetProcessSubType()
    1148            << G4endl
    1149            << "      dE/dx and range tables from "
    1150            << G4BestUnit(minKinEnergy,"Energy")
    1151            << " to " << G4BestUnit(maxKinEnergy,"Energy")
    1152            << " in " << nBins << " bins" << G4endl
    1153            << "      Lambda tables from threshold to "
    1154            << G4BestUnit(maxKinEnergy,"Energy")
    1155            << " in " << nBins << " bins, spline: "
    1156            << (G4LossTableManager::Instance())->SplineFlag()
     1333G4bool G4VEnergyLossProcess::StorePhysicsTable(
     1334       const G4ParticleDefinition* part, const G4String& directory,
     1335       G4bool ascii)
     1336{
     1337  G4bool res = true;
     1338  if ( baseParticle || part != particle ) return res;
     1339
     1340  if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX"))
     1341    {res = false;}
     1342
     1343  if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr"))
     1344    {res = false;}
     1345
     1346  if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX"))
     1347    {res = false;}
     1348
     1349  if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation"))
     1350    {res = false;}
     1351
     1352  if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation"))
     1353    {res = false;}
     1354
     1355  if(isIonisation &&
     1356     !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange"))
     1357    {res = false;}
     1358
     1359  if(isIonisation &&
     1360     !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range"))
     1361    {res = false;}
     1362 
     1363  if(isIonisation &&
     1364     !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange"))
     1365    {res = false;}
     1366 
     1367  if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda"))
     1368    {res = false;}
     1369
     1370  if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda"))
     1371    {res = false;}
     1372
     1373  if ( res ) {
     1374    if(0 < verboseLevel) {
     1375      G4cout << "Physics tables are stored for " << particle->GetParticleName()
     1376             << " and process " << GetProcessName()
     1377             << " in the directory <" << directory
     1378             << "> " << G4endl;
     1379    }
     1380  } else {
     1381    G4cout << "Fail to store Physics Tables for "
     1382           << particle->GetParticleName()
     1383           << " and process " << GetProcessName()
     1384           << " in the directory <" << directory
     1385           << "> " << G4endl;
     1386  }
     1387  return res;
     1388}
     1389
     1390//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1391
     1392G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
     1393       const G4ParticleDefinition* part, const G4String& directory,
     1394       G4bool ascii)
     1395{
     1396  G4bool res = true;
     1397  const G4String particleName = part->GetParticleName();
     1398
     1399  if(1 < verboseLevel) {
     1400    G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
     1401           << particleName << " and process " << GetProcessName()
     1402           << "; tables_are_built= " << tablesAreBuilt
    11571403           << G4endl;
    1158     if(theRangeTableForLoss && isIonisation) {
    1159       G4cout << "      finalRange(mm)= " << finalRange/mm
    1160              << ", dRoverRange= " << dRoverRange
    1161              << ", integral: " << integral
    1162              << ", fluct: " << lossFluctuationFlag
    1163              << ", linLossLimit= " << linLossLimit
    1164              << G4endl;
    1165     }
    1166     PrintInfo();
    1167     modelManager->DumpModelList(verboseLevel);
    1168     if(theCSDARangeTable && isIonisation) {
    1169       G4cout << "      CSDA range table up"
    1170              << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
    1171              << " in " << nBinsCSDA << " bins" << G4endl;
    1172     }
    1173     if(nSCoffRegions>0 && isIonisation) {
    1174       G4cout << "      Subcutoff sampling in " << nSCoffRegions
    1175              << " regions" << G4endl;
    1176     }
    1177     if(2 < verboseLevel) {
    1178       G4cout << "      DEDXTable address= " << theDEDXTable << G4endl;
    1179       if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
    1180       G4cout << "non restricted DEDXTable address= "
    1181              << theDEDXunRestrictedTable << G4endl;
    1182       if(theDEDXunRestrictedTable && isIonisation) {
    1183            G4cout << (*theDEDXunRestrictedTable) << G4endl;
    1184       }
    1185       if(theDEDXSubTable && isIonisation) {
    1186         G4cout << (*theDEDXSubTable) << G4endl;
    1187       }
    1188       G4cout << "      CSDARangeTable address= " << theCSDARangeTable
     1404  }
     1405  if(particle == part) {
     1406
     1407    //    G4bool yes = true;
     1408    if ( !baseParticle ) {
     1409
     1410      G4bool fpi = true;
     1411      if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi))
     1412        {fpi = false;}
     1413
     1414      if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false))
     1415        {fpi = false;}
     1416
     1417      if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi))
     1418        {res = false;}
     1419
     1420      if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false))
     1421        {res = false;}
     1422
     1423      if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false))
     1424        {res = false;}
     1425
     1426      if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi))
     1427        {res = false;}
     1428
     1429      if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true))
     1430        {res = false;}
     1431
     1432      G4bool yes = false;
     1433      if(nSCoffRegions > 0) {yes = true;}
     1434
     1435      if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes))
     1436        {res = false;}
     1437
     1438      if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes))
     1439        {res = false;}
     1440
     1441      if(!fpi) yes = false;
     1442      if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes))
     1443        {res = false;}
     1444    }
     1445  }
     1446
     1447  return res;
     1448}
     1449
     1450//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1451
     1452G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part,
     1453                                        G4PhysicsTable* aTable, G4bool ascii,
     1454                                        const G4String& directory,
     1455                                        const G4String& tname)
     1456{
     1457  G4bool res = true;
     1458  if ( aTable ) {
     1459    const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
     1460    if( !aTable->StorePhysicsTable(name,ascii)) res = false;
     1461  }
     1462  return res;
     1463}
     1464
     1465//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1466
     1467G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part,
     1468                                           G4PhysicsTable* aTable, G4bool ascii,
     1469                                           const G4String& directory,
     1470                                           const G4String& tname,
     1471                                           G4bool mandatory)
     1472{
     1473  G4bool res = true;
     1474  G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
     1475  G4bool yes = aTable->ExistPhysicsTable(filename);
     1476  if(yes) {
     1477    yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii);
     1478    if((G4LossTableManager::Instance())->SplineFlag()) {
     1479      size_t n = aTable->length();
     1480      for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);}
     1481    }
     1482  }
     1483  if(yes) {
     1484    if (0 < verboseLevel) {
     1485      G4cout << tname << " table for " << part->GetParticleName()
     1486             << " is Retrieved from <" << filename << ">"
    11891487             << G4endl;
    1190       if(theCSDARangeTable && isIonisation) {
    1191         G4cout << (*theCSDARangeTable) << G4endl;
    1192       }
    1193       G4cout << "      RangeTableForLoss address= " << theRangeTableForLoss
     1488    }
     1489  } else {
     1490    if(mandatory) res = false;
     1491    if(mandatory || 1 < verboseLevel) {
     1492      G4cout << tname << " table for " << part->GetParticleName()
     1493             << " from file <"
     1494             << filename << "> is not Retrieved"
    11941495             << G4endl;
    1195       if(theRangeTableForLoss && isIonisation) {
    1196              G4cout << (*theRangeTableForLoss) << G4endl;
    1197       }
    1198       G4cout << "      InverseRangeTable address= " << theInverseRangeTable
    1199              << G4endl;
    1200       if(theInverseRangeTable && isIonisation) {
    1201              G4cout << (*theInverseRangeTable) << G4endl;
    1202       }
    1203       G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
    1204       if(theLambdaTable && isIonisation) {
    1205         G4cout << (*theLambdaTable) << G4endl;
    1206       }
    1207       G4cout << "      SubLambdaTable address= " << theSubLambdaTable << G4endl;
    1208       if(theSubLambdaTable && isIonisation) {
    1209         G4cout << (*theSubLambdaTable) << G4endl;
    1210       }
     1496    }
     1497  }
     1498  return res;
     1499}
     1500
     1501//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1502
     1503G4double G4VEnergyLossProcess::GetDEDXDispersion(
     1504                                  const G4MaterialCutsCouple *couple,
     1505                                  const G4DynamicParticle* dp,
     1506                                        G4double length)
     1507{
     1508  DefineMaterial(couple);
     1509  G4double ekin = dp->GetKineticEnergy();
     1510  SelectModel(ekin*massRatio);
     1511  G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
     1512  tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]);
     1513  G4double d = 0.0;
     1514  G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
     1515  if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length);
     1516  return d;
     1517}
     1518
     1519//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1520
     1521G4double G4VEnergyLossProcess::CrossSectionPerVolume(
     1522         G4double kineticEnergy, const G4MaterialCutsCouple* couple)
     1523{
     1524  // Cross section per volume is calculated
     1525  DefineMaterial(couple);
     1526  G4double cross = 0.0;
     1527  G4bool b;
     1528  if(theLambdaTable) {
     1529    cross =
     1530      ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b);
     1531  } else {
     1532    SelectModel(kineticEnergy);
     1533    cross =
     1534      currentModel->CrossSectionPerVolume(currentMaterial,
     1535                                          particle, kineticEnergy,
     1536                                          (*theCuts)[currentMaterialIndex]);
     1537  }
     1538
     1539  return cross;
     1540}
     1541
     1542//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1543
     1544G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track)
     1545{
     1546  DefineMaterial(track.GetMaterialCutsCouple());
     1547  preStepLambda = GetLambdaForScaledEnergy(track.GetKineticEnergy()*massRatio);
     1548  G4double x = DBL_MAX;
     1549  if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
     1550  return x;
     1551}
     1552
     1553//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1554
     1555G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track,
     1556                                                   G4double x, G4double y,
     1557                                                   G4double& z)
     1558{
     1559  G4GPILSelection sel;
     1560  return AlongStepGetPhysicalInteractionLength(track, x, y, z, &sel);
     1561}
     1562
     1563//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1564
     1565G4double G4VEnergyLossProcess::GetMeanFreePath(
     1566                             const G4Track& track,
     1567                             G4double,
     1568                             G4ForceCondition* condition)
     1569
     1570{
     1571  *condition = NotForced;
     1572  return MeanFreePath(track);
     1573}
     1574
     1575//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1576
     1577G4double G4VEnergyLossProcess::GetContinuousStepLimit(
     1578                const G4Track&,
     1579                G4double, G4double, G4double&)
     1580{
     1581  return DBL_MAX;
     1582}
     1583
     1584//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1585
     1586G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
     1587                 const G4MaterialCutsCouple* couple, G4double cut)
     1588{
     1589  //  G4double cut  = (*theCuts)[couple->GetIndex()];
     1590  //  G4int nbins = nLambdaBins;
     1591  G4double tmin =
     1592    std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
     1593             minKinEnergy);
     1594  if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
     1595  G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
     1596  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
     1597
     1598  return v;
     1599}
     1600
     1601//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1602 
     1603void G4VEnergyLossProcess::AddCollaborativeProcess(
     1604            G4VEnergyLossProcess* p)
     1605{
     1606  G4bool add = true;
     1607  if(p->GetProcessName() != "eBrem") add = false;
     1608  if(add && nProcesses > 0) {
     1609    for(G4int i=0; i<nProcesses; i++) {
     1610      if(p == scProcesses[i]) {
     1611        add = false;
     1612        break;
     1613      }
     1614    }
     1615  }
     1616  if(add) {
     1617    scProcesses.push_back(p);
     1618    nProcesses++;
     1619    if (1 < verboseLevel) {
     1620      G4cout << "### The process " << p->GetProcessName()
     1621             << " is added to the list of collaborative processes of "
     1622             << GetProcessName() << G4endl;
    12111623    }
    12121624  }
    … …  
    13771789//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    13781790
    1379 G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
    1380                  const G4MaterialCutsCouple* couple, G4double cut)
    1381 {
    1382   //  G4double cut  = (*theCuts)[couple->GetIndex()];
    1383   //  G4int nbins = nLambdaBins;
    1384   G4double tmin =
    1385     std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
    1386              minKinEnergy);
    1387   if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
    1388   G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
    1389   v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
    1390 
    1391   return v;
    1392 }
    1393 
    1394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1395 
    1396 G4double G4VEnergyLossProcess::CrossSectionPerVolume(
    1397          G4double kineticEnergy, const G4MaterialCutsCouple* couple)
    1398 {
    1399   // Cross section per volume is calculated
    1400   DefineMaterial(couple);
    1401   G4double cross = 0.0;
    1402   G4bool b;
    1403   if(theLambdaTable) {
    1404     cross =
    1405       ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b);
    1406   } else {
    1407     SelectModel(kineticEnergy);
    1408     cross =
    1409       currentModel->CrossSectionPerVolume(currentMaterial,
    1410                                           particle, kineticEnergy,
    1411                                           (*theCuts)[currentMaterialIndex]);
    1412   }
    1413 
    1414   return cross;
    1415 }
    1416 
    1417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1418 
    1419 G4bool G4VEnergyLossProcess::StorePhysicsTable(
    1420        const G4ParticleDefinition* part, const G4String& directory,
    1421        G4bool ascii)
    1422 {
    1423   G4bool res = true;
    1424   if ( baseParticle || part != particle ) return res;
    1425 
    1426   if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX"))
    1427     {res = false;}
    1428 
    1429   if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr"))
    1430     {res = false;}
    1431 
    1432   if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX"))
    1433     {res = false;}
    1434 
    1435   if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation"))
    1436     {res = false;}
    1437 
    1438   if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation"))
    1439     {res = false;}
    1440 
    1441   if(isIonisation &&
    1442      !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange"))
    1443     {res = false;}
    1444 
    1445   if(isIonisation &&
    1446      !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range"))
    1447     {res = false;}
    1448  
    1449   if(isIonisation &&
    1450      !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange"))
    1451     {res = false;}
    1452  
    1453   if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda"))
    1454     {res = false;}
    1455 
    1456   if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda"))
    1457     {res = false;}
    1458 
    1459   if ( res ) {
    1460     if(0 < verboseLevel) {
    1461       G4cout << "Physics tables are stored for " << particle->GetParticleName()
    1462              << " and process " << GetProcessName()
    1463              << " in the directory <" << directory
    1464              << "> " << G4endl;
    1465     }
    1466   } else {
    1467     G4cout << "Fail to store Physics Tables for "
    1468            << particle->GetParticleName()
    1469            << " and process " << GetProcessName()
    1470            << " in the directory <" << directory
    1471            << "> " << G4endl;
    1472   }
    1473   return res;
    1474 }
    1475 
    1476 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    1477 
    1478 G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part,
    1479                                         G4PhysicsTable* aTable, G4bool ascii,
    1480                                         const G4String& directory,
    1481                                         const G4String& tname)
    1482 {
    1483   G4bool res = true;
    1484   if ( aTable ) {
    1485     const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
    1486     if( !aTable->StorePhysicsTable(name,ascii)) res = false;
    1487   }
    1488   return res;
    1489 }
    1490 
    1491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    1492 
    1493 G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
    1494        const G4ParticleDefinition* part, const G4String& directory,
    1495        G4bool ascii)
    1496 {
    1497   G4bool res = true;
    1498   const G4String particleName = part->GetParticleName();
    1499 
    1500   if(1 < verboseLevel) {
    1501     G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
    1502            << particleName << " and process " << GetProcessName()
    1503            << "; tables_are_built= " << tablesAreBuilt
    1504            << G4endl;
    1505   }
    1506   if(particle == part) {
    1507 
    1508     //    G4bool yes = true;
    1509     if ( !baseParticle ) {
    1510 
    1511       G4bool fpi = true;
    1512       if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi))
    1513         {fpi = false;}
    1514 
    1515       if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false))
    1516         {fpi = false;}
    1517 
    1518       if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi))
    1519         {res = false;}
    1520 
    1521       if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false))
    1522         {res = false;}
    1523 
    1524       if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false))
    1525         {res = false;}
    1526 
    1527       if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi))
    1528         {res = false;}
    1529 
    1530       if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true))
    1531         {res = false;}
    1532 
    1533       G4bool yes = false;
    1534       if(nSCoffRegions > 0) {yes = true;}
    1535 
    1536       if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes))
    1537         {res = false;}
    1538 
    1539       if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes))
    1540         {res = false;}
    1541 
    1542       if(!fpi) yes = false;
    1543       if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes))
    1544         {res = false;}
    1545     }
    1546   }
    1547 
    1548   return res;
    1549 }
    1550 
    1551 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    1552 
    1553 G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part,
    1554                                            G4PhysicsTable* aTable, G4bool ascii,
    1555                                            const G4String& directory,
    1556                                            const G4String& tname,
    1557                                            G4bool mandatory)
    1558 {
    1559   G4bool res = true;
    1560   G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
    1561   G4bool yes = aTable->ExistPhysicsTable(filename);
    1562   if(yes) {
    1563     yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii);
    1564     if((G4LossTableManager::Instance())->SplineFlag()) {
    1565       size_t n = aTable->length();
    1566       for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);}
    1567     }
    1568   }
    1569   if(yes) {
    1570     if (0 < verboseLevel) {
    1571       G4cout << tname << " table for " << part->GetParticleName()
    1572              << " is Retrieved from <" << filename << ">"
    1573              << G4endl;
    1574     }
    1575   } else {
    1576     if(mandatory) res = false;
    1577     if(mandatory || 1 < verboseLevel) {
    1578       G4cout << tname << " table for " << part->GetParticleName()
    1579              << " from file <"
    1580              << filename << "> is not Retrieved"
    1581              << G4endl;
    1582     }
    1583   }
    1584   return res;
    1585 }
    1586 
    1587 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1588  
    1589 void G4VEnergyLossProcess::AddCollaborativeProcess(
    1590             G4VEnergyLossProcess* p)
    1591 {
    1592   G4bool add = true;
    1593   if(p->GetProcessName() != "eBrem") add = false;
    1594   if(add && nProcesses > 0) {
    1595     for(G4int i=0; i<nProcesses; i++) {
    1596       if(p == scProcesses[i]) {
    1597         add = false;
    1598         break;
    1599       }
    1600     }
    1601   }
    1602   if(add) {
    1603     scProcesses.push_back(p);
    1604     nProcesses++;
    1605     if (1 < verboseLevel) {
    1606       G4cout << "### The process " << p->GetProcessName()
    1607              << " is added to the list of collaborative processes of "
    1608              << GetProcessName() << G4endl;
    1609     }
    1610   }
    1611 }
    1612 
    1613 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1614 
    1615 G4double G4VEnergyLossProcess::GetDEDXDispersion(
    1616                                   const G4MaterialCutsCouple *couple,
    1617                                   const G4DynamicParticle* dp,
    1618                                         G4double length)
    1619 {
    1620   DefineMaterial(couple);
    1621   G4double ekin = dp->GetKineticEnergy();
    1622   SelectModel(ekin*massRatio);
    1623   G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
    1624   tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]);
    1625   G4double d = 0.0;
    1626   G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
    1627   if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length);
    1628   return d;
    1629 }
    1630 
    1631 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1632 
    1633 void G4VEnergyLossProcess::ActivateDeexcitation(G4bool, const G4Region*)
    1634 {}
    1635 
    1636 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1637 

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMscModel.cc,v 1.4 2008/03/10 18:39:45 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VMscModel.cc,v 1.12 2009/05/10 19:36:19 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4949
    5050#include "G4VMscModel.hh"
     51#include "G4ParticleChangeForMSC.hh"
     52#include "G4TransportationManager.hh"
    5153
    5254//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    5557G4VMscModel::G4VMscModel(const G4String& nam):
    5658  G4VEmModel(nam),
    57   facrange(0.02),
     59  safetyHelper(0),
     60  facrange(0.04),
    5861  facgeom(2.5),
    5962  facsafety(0.25),
    … …  
    6164  dtrl(0.05),
    6265  lambdalimit(mm),
     66  geommax(1.e50*mm),
    6367  steppingAlgorithm(fUseSafety),
    6468  samplez(false),
    … …  
    7276
    7377//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     78
     79G4ParticleChangeForMSC* G4VMscModel::GetParticleChangeForMSC()
     80{
     81  G4ParticleChangeForMSC* p = 0;
     82  if (pParticleChange) {
     83    p = static_cast<G4ParticleChangeForMSC*>(pParticleChange);
     84  } else {
     85    p = new G4ParticleChangeForMSC();
     86  }
     87  return p;
     88}
     89
     90//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     91
     92void G4VMscModel::InitialiseSafetyHelper()
     93{
     94  if(!safetyHelper) {
     95    safetyHelper = G4TransportationManager::GetTransportationManager()
     96      ->GetSafetyHelper();
     97    safetyHelper->InitialiseHelper();
     98  }
     99}
     100
     101//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     102
     103void G4VMscModel::ComputeDisplacement(G4ParticleChangeForMSC* fParticleChange, 
     104                                      const G4ThreeVector& dir,
     105                                      G4double displacement,
     106                                      G4double postsafety)
     107{
     108  const G4ThreeVector* pos = fParticleChange->GetProposedPosition();
     109  G4double r = displacement;
     110  if(r >  postsafety) {
     111    G4double newsafety = safetyHelper->ComputeSafety(*pos);
     112    if(r > newsafety) r = newsafety;
     113  }
     114  if(r > 0.) {
     115
     116    // compute new endpoint of the Step
     117    G4ThreeVector newPosition = *pos + r*dir;
     118
     119    // definitely not on boundary
     120    if(displacement == r) {
     121      safetyHelper->ReLocateWithinVolume(newPosition);
     122
     123    } else {
     124      // check safety after displacement
     125      G4double postsafety = safetyHelper->ComputeSafety(newPosition);
     126
     127      // displacement to boundary
     128      if(postsafety <= 0.0) {
     129        safetyHelper->Locate(newPosition,
     130                             *fParticleChange->GetProposedMomentumDirection());
     131
     132        // not on the boundary
     133      } else {
     134        safetyHelper->ReLocateWithinVolume(newPosition);
     135      }
     136    }
     137    fParticleChange->ProposePosition(newPosition);
     138  }
     139}
     140
     141//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     142
     143void G4VMscModel::SampleScattering(const G4DynamicParticle*, G4double)
     144{}
     145
     146//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     147
     148G4double G4VMscModel::ComputeTruePathLengthLimit(const G4Track&,
     149                                                 G4PhysicsTable*,
     150                                                 G4double)
     151{
     152  return DBL_MAX;
     153}
     154
     155//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     156
     157G4double G4VMscModel::ComputeGeomPathLength(G4double truePathLength)
     158{
     159  return truePathLength;
     160}
     161
     162//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     163
     164G4double G4VMscModel::ComputeTrueStepLength(G4double geomPathLength)
     165{
     166  return geomPathLength;
     167}
     168
     169//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     170
     171void G4VMscModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
     172                                    const G4MaterialCutsCouple*,
     173                                    const G4DynamicParticle*,
     174                                    G4double, G4double)
     175{}
     176
     177//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMultipleScattering.cc,v 1.60 2008/11/20 20:32:40 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4VMultipleScattering.cc,v 1.66 2009/05/27 11:55:02 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8383#include "G4GenericIon.hh"
    8484#include "G4Electron.hh"
     85#include "G4EmConfigurator.hh"
    8586
    8687//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    9495  stepLimit(fUseSafety),
    9596  skin(3.0),
    96   facrange(0.02),
     97  facrange(0.04),
    9798  facgeom(2.5),
    9899  latDisplasment(true),
    … …  
    105106  // Size of tables assuming spline
    106107  minKinEnergy = 0.1*keV;
    107   maxKinEnergy = 100.0*TeV;
    108   nBins        = 84;
     108  maxKinEnergy = 10.0*TeV;
     109  nBins        = 77;
    109110
    110111  // default limit on polar angle
    … …  
    132133  }
    133134  (G4LossTableManager::Instance())->DeRegister(this);
     135}
     136
     137//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     138
     139void G4VMultipleScattering::AddEmModel(G4int order, G4VEmModel* p,
     140                                       const G4Region* region)
     141{
     142  G4VEmFluctuationModel* fm = 0;
     143  modelManager->AddEmModel(order, p, fm, region);
     144  if(p) p->SetParticleChange(pParticleChange);
     145}
     146
     147//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     148
     149G4VEmModel* G4VMultipleScattering::GetModelByIndex(G4int idx, G4bool ver)
     150{
     151  return modelManager->GetModel(idx, ver);
    134152}
    135153
    … …  
    208226           << G4endl;
    209227  }
     228
     229  (G4LossTableManager::Instance())->EmConfigurator()->AddModels();
    210230
    211231  if(firstParticle == &part) {
Note: See TracChangeset for help on using the changeset viewer.