Ignore:
Timestamp:
Apr 20, 2009, 4:53:50 PM (15 years ago)
Author:
garnier
Message:

fichier oublies

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

Legend:

Unmodified
Added
Removed
  • trunk/source/processes/electromagnetic/utils/History

    r991 r1005  
    1 $Id: History,v 1.364 2008/11/20 20:32:40 vnivanch Exp $
     1$Id: History,v 1.372 2009/02/26 11:33:33 vnivanch Exp $
    22-------------------------------------------------------------------
    33
     
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     2026 February 09: V.Ivant (emutils-V09-02-03)
     21G4EmConfigurator - fixed for the case if only fluctuation model is set
     22                   and main model is default
     23
     2422 February 09: V.Ivant (emutils-V09-02-02)
     25- G4VEmModel - make methods to access geometry protected, added new
     26               method SetSampleZ, added geommax private member
     27- G4EmCalculator - added possibility to be used by DNA processes:
     28                   take into account special DNA particles
     29
     3018 February 09: V.Ivant (emutils-V09-02-01)
     31G4VEmModel, G4VEmFluctuationModel, G4VEnegryLossProcess, G4VEmProcess,
     32G4VMultipleScattering - move all virtual methods to source, update comments
     33G4VEmModel - added flagDeexcitation and Get/Set methods
     34G4VEnegryLossProcess, G4VEmProcess - added calls to deexcitation PostStep
     35G4EmProcessOptions - added ActivateDeexcitation method
     36G4EnergyLossMessenger - added /process/em/deexcitation UI command
     37G4LossTableBuilder - added protection in BuildRangeTable against zero dedx
     38
     3927 January 09: V.Ivant (emutils-V09-02-00)
     40G4VEmModel - added method SampleDeexcitationAlongStep
     41G4VEnegryLossProcess - added deexcitation AlongStep per region
     42G4VMscModel - added methdos: InitialiseSafetyHelper, ComputeSafety,
     43              ComputeGeomLimit, ComputeDisplacement
     44G4VEmProcess - added possibility to set more than 1 model
    1945
    204620 November 08: V.Ivant (emutils-V09-01-37)
  • trunk/source/processes/electromagnetic/utils/include/G4EmProcessOptions.hh

    r991 r1005  
    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-02-ref-02 $
    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

    r991 r1005  
    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-02-ref-02 $
    2929//
    3030// -------------------------------------------------------------------
     
    107107  G4UIcmdWithADouble*        MinSubSecCmd;
    108108  G4UIcommand*               StepFuncCmd;
     109  G4UIcommand*               deexCmd;
    109110  G4UIcmdWithAString*        mscCmd;
    110111  G4UIcmdWithADoubleAndUnit* MinEnCmd;
  • trunk/source/processes/electromagnetic/utils/include/G4VEmFluctuationModel.hh

    r991 r1005  
    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-02-ref-02 $
    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

    r991 r1005  
    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.66 2009/02/19 09:57:36 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
     
    6565// 12-09-08 Added methods GetParticleCharge, GetChargeSquareRatio,
    6666//          CorrectionsAlongStep, ActivateNuclearStopping (VI)
     67// 16-02-09 Move implementations of virtual methods to source (VI)
    6768//
    6869// Class Description:
     
    120121  //------------------------------------------------------------------------
    121122
    122   // dEdx per unit length
    123   virtual G4double ComputeDEDX(const G4MaterialCutsCouple*,
    124                                const G4ParticleDefinition*,
    125                                G4double kineticEnergy,
    126                                G4double cutEnergy = DBL_MAX);
    127 
    128123  // main method to compute dEdx
    129124  virtual G4double ComputeDEDXPerVolume(const G4Material*,
     
    131126                                        G4double kineticEnergy,
    132127                                        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);
    140128
    141129  // main method to compute cross section per Volume
     
    175163                                    G4double length);
    176164
     165  // sample PIXE deexcitation
     166  virtual void SampleDeexcitationAlongStep(const G4Material*,
     167                                           const G4Track&,
     168                                           G4double& eloss);
     169
    177170protected:
    178171
     
    212205                                  const G4DataVector&);
    213206
     207  // dEdx per unit length
     208  inline G4double ComputeDEDX(const G4MaterialCutsCouple*,
     209                              const G4ParticleDefinition*,
     210                              G4double kineticEnergy,
     211                              G4double cutEnergy = DBL_MAX);
     212
     213  // cross section per volume
     214  inline G4double CrossSection(const G4MaterialCutsCouple*,
     215                                const G4ParticleDefinition*,
     216                                G4double kineticEnergy,
     217                                G4double cutEnergy = 0.0,
     218                                G4double maxEnergy = DBL_MAX);
     219
    214220  // 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);
     221  inline G4double ComputeMeanFreePath(const G4ParticleDefinition*,
     222                                      G4double kineticEnergy,
     223                                      const G4Material*,   
     224                                      G4double cutEnergy = 0.0,
     225                                      G4double maxEnergy = DBL_MAX);
    220226
    221227  // generic cross section per element
     
    233239                                           G4double maxEnergy = DBL_MAX);
    234240
    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
     241  // to select atom cross section per volume is recomputed for each element
    237242  inline const G4Element* SelectRandomAtom(const G4Material*,
    238243                                           const G4ParticleDefinition*,
     
    260265  inline G4bool LPMFlag() const;
    261266
     267  inline G4bool DeexcitationFlag() const;
     268
    262269  inline void SetHighEnergyLimit(G4double);
    263270
     
    270277  inline void SetLPMFlag(G4bool val);
    271278
     279  inline void SetDeexcitationFlag(G4bool val);
     280
    272281  inline void ActivateNuclearStopping(G4bool);
    273282
     
    278287  inline void SetParticleChange(G4VParticleChange*, G4VEmFluctuationModel*);
    279288
     289  inline void SetCurrentCouple(const G4MaterialCutsCouple*);
     290
    280291protected:
    281292
     293  inline const G4MaterialCutsCouple* CurrentCouple() const;
     294
     295  inline void SetCurrentElement(const G4Element*);
     296
    282297  inline const G4Element* GetCurrentElement() const;
    283 
    284   inline void SetCurrentElement(const G4Element*);
    285298
    286299private:
     
    315328private:
    316329
    317   const G4Element* currentElement;
     330  const G4MaterialCutsCouple* currentCouple;
     331  const G4Element*            currentElement;
     332
    318333  G4int                  nsec;
     334  G4bool                 flagDeexcitation;
    319335  std::vector<G4double>  xsec;
    320336
     
    324340//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    325341
    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;
     342inline G4double G4VEmModel::ComputeDEDX(const G4MaterialCutsCouple* c,
     343                                        const G4ParticleDefinition* p,
     344                                        G4double kinEnergy,
     345                                        G4double cutEnergy)
     346{
     347  currentCouple = c;
     348  return ComputeDEDXPerVolume(c->GetMaterial(),p,kinEnergy,cutEnergy);
     349}
     350
     351//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     352
     353inline G4double G4VEmModel::CrossSection(const G4MaterialCutsCouple* c,
     354                                         const G4ParticleDefinition* p,
     355                                         G4double kinEnergy,
     356                                         G4double cutEnergy,
     357                                         G4double maxEnergy)
     358{
     359  currentCouple = c;
     360  return CrossSectionPerVolume(c->GetMaterial(),p,kinEnergy,cutEnergy,maxEnergy);
     361}
     362
     363//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     364
     365inline G4double G4VEmModel::ComputeMeanFreePath(const G4ParticleDefinition* p,
     366                                                G4double ekin,
     367                                                const G4Material* material,     
     368                                                G4double emin,
     369                                                G4double emax)
     370{
     371  G4double mfp = DBL_MAX;
     372  G4double cross = CrossSectionPerVolume(material,p,ekin,emin,emax);
     373  if (cross > DBL_MIN) mfp = 1./cross;
     374  return mfp;
    399375}
    400376
     
    415391//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    416392
    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 
    507393inline
    508394const G4Element* G4VEmModel::SelectRandomAtom(const G4MaterialCutsCouple* couple,
     
    512398                                              G4double maxEnergy)
    513399{
     400  currentCouple = couple;
    514401  if(nSelectors > 0) {
    515402    currentElement =
     
    571458//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    572459
    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;
     460inline G4VEmFluctuationModel* G4VEmModel::GetModelOfFluctuations()
     461{
     462  return fluc;
     463}
     464
     465//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     466
     467inline G4double G4VEmModel::HighEnergyLimit() const
     468{
     469  return highLimit;
     470}
     471
     472//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     473
     474inline G4double G4VEmModel::LowEnergyLimit() const
     475{
     476  return lowLimit;
     477}
     478
     479//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     480
     481inline G4double G4VEmModel::PolarAngleLimit() const
     482{
     483  return polarAngleLimit;
     484}
     485
     486//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     487
     488inline G4double G4VEmModel::SecondaryThreshold() const
     489{
     490  return secondaryThreshold;
     491}
     492
     493//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     494
     495inline G4bool G4VEmModel::LPMFlag() const
     496{
     497  return theLPMflag;
     498}
     499
     500//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     501
     502inline G4bool G4VEmModel::DeexcitationFlag() const
     503{
     504  return flagDeexcitation;
     505}
     506
     507//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     508
     509inline void G4VEmModel::SetHighEnergyLimit(G4double val)
     510{
     511  highLimit = val;
     512}
     513
     514//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     515
     516inline void G4VEmModel::SetLowEnergyLimit(G4double val)
     517{
     518  lowLimit = val;
     519}
     520
     521//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     522
     523inline void G4VEmModel::SetPolarAngleLimit(G4double val)
     524{
     525  polarAngleLimit = val;
     526}
     527
     528//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     529
     530inline void G4VEmModel::SetSecondaryThreshold(G4double val)
     531{
     532  secondaryThreshold = val;
     533}
     534
     535//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     536
     537inline void G4VEmModel::SetLPMFlag(G4bool val)
     538{
     539  theLPMflag = val;
     540}
     541
     542//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     543
     544inline void G4VEmModel::SetDeexcitationFlag(G4bool val)
     545{
     546  flagDeexcitation = val;
     547}
     548
     549//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     550
     551inline void G4VEmModel::ActivateNuclearStopping(G4bool val)
     552{
     553  nuclearStopping = val;
    583554}
    584555
     
    594565//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    595566
    596 inline G4double G4VEmModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
    597                                                G4double kineticEnergy)
    598 {
    599   return kineticEnergy;
    600 }
    601 
    602 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    603 
    604567inline const G4String& G4VEmModel::GetName() const
    605568{
     
    608571
    609572//....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 {}
     573
     574inline void G4VEmModel::SetParticleChange(G4VParticleChange* p, 
     575                                          G4VEmFluctuationModel* f = 0)
     576{
     577  if(p && pParticleChange != p) pParticleChange = p;
     578  fluc = f;
     579}
     580
     581//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     582
     583inline void G4VEmModel::SetCurrentCouple(const G4MaterialCutsCouple* p)
     584{
     585  currentCouple = p;
     586}
     587
     588//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     589
     590inline const G4MaterialCutsCouple* G4VEmModel::CurrentCouple() const
     591{
     592  return currentCouple;
     593}
     594
     595//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     596
     597inline void G4VEmModel::SetCurrentElement(const G4Element* elm)
     598{
     599  currentElement = elm;
     600}
     601
     602//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     603
     604inline const G4Element* G4VEmModel::GetCurrentElement() const
     605{
     606  return currentElement;
     607}
    650608
    651609//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  • trunk/source/processes/electromagnetic/utils/include/G4VEmProcess.hh

    r991 r1005  
    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.50 2009/02/19 09:57:36 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    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
     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;
     210   
     211  // Add model for region, smaller value of order defines which
     212  // model will be selected for a given energy interval 
    209213  inline void AddEmModel(G4int, G4VEmModel*, const G4Region* region = 0);
    210    
     214
    211215  // Assign a model to a process
    212   inline void SetModel(G4VEmModel*);
     216  inline void SetModel(G4VEmModel*, G4int index = 1);
    213217 
    214218  // return the assigned model
    215   inline G4VEmModel* Model();
     219  inline G4VEmModel* Model(G4int index = 1);
    216220   
    217221  // Define new energy range for the model identified by the name
     
    221225  inline G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false);
    222226
    223   //------------------------------------------------------------------------
    224   // Get/set parameters used for simulation of energy loss
    225   //------------------------------------------------------------------------
    226 
    227   inline void ActivateDeexcitation(G4bool, const G4Region* r = 0);
    228 
    229227  inline void SetLambdaFactor(G4double val);
    230228
     
    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::AddEmModel(G4int order, G4VEmModel* p,
     487                                     const G4Region* region)
     488{
     489  G4VEmFluctuationModel* fm = 0;
     490  modelManager->AddEmModel(order, p, fm, region);
     491  if(p) p->SetParticleChange(pParticleChange);
     492}
     493
     494//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     495
     496inline void G4VEmProcess::SetModel(G4VEmModel* p, G4int index)
     497{
     498  G4int n = emModels.size();
     499  if(index >= n) for(G4int i=n; i<index+1; i++) {emModels.push_back(0);}
     500  emModels[index] = p;
     501}
     502
     503//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     504
     505inline G4VEmModel* G4VEmProcess::Model(G4int index)
     506{
     507  G4VEmModel* p = 0;
     508  if(index >= 0 && index <  G4int(emModels.size())) p = emModels[index];
     509  return p;
     510}
     511
     512//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     513
     514inline void G4VEmProcess::UpdateEmModel(const G4String& nam,
     515                                        G4double emin, G4double emax)
     516{
     517  modelManager->UpdateEmModel(nam, emin, emax);
     518}
     519
     520//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     521
     522inline G4VEmModel* G4VEmProcess::GetModelByIndex(G4int idx, G4bool ver)
     523{
     524  return modelManager->GetModel(idx, ver);
     525}
     526
     527//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     528
     529inline void G4VEmProcess::SetLambdaFactor(G4double val)
     530{
     531  if(val > 0.0 && val <= 1.0) lambdaFactor = val;
     532}
     533
     534//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     535
     536inline void G4VEmProcess::SetIntegral(G4bool val)
     537{
     538  if(particle && particle != theGamma) integral = val;
     539  if(integral) buildLambdaTable = true;
     540}
     541
     542//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     543
     544inline G4bool G4VEmProcess::IsIntegral() const
     545{
     546  return integral;
     547}
     548
     549//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     550
     551inline void G4VEmProcess::SetApplyCuts(G4bool val)
     552{
     553  applyCuts = val;
     554}
     555
     556//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     557
     558inline G4double G4VEmProcess::RecalculateLambda(G4double e,
     559                                                const G4MaterialCutsCouple* couple)
     560{
     561  DefineMaterial(couple);
     562  return ComputeCurrentLambda(e);
     563}
     564
     565//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     566
     567inline G4ParticleChangeForGamma* G4VEmProcess::GetParticleChange()
     568{
     569  return &fParticleChange;
     570}
     571
     572//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     573
     574inline void G4VEmProcess::SetParticle(const G4ParticleDefinition* p)
     575{
     576  particle = p;
     577}
     578
     579//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     580
     581inline void G4VEmProcess::SetSecondaryParticle(const G4ParticleDefinition* p)
     582{
     583  secondaryParticle = p;
     584}
     585
     586//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     587
     588inline size_t G4VEmProcess::CurrentMaterialCutsCoupleIndex() const
     589{
     590  return currentMaterialIndex;
     591}
     592
     593//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     594
     595inline G4double G4VEmProcess::GetGammaEnergyCut()
     596{
     597  return (*theCutsGamma)[currentMaterialIndex];
     598}
     599
     600//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     601
     602inline G4double G4VEmProcess::GetElectronEnergyCut()
     603{
     604  return (*theCutsElectron)[currentMaterialIndex];
     605}
     606
     607//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     608
     609inline void G4VEmProcess::SetBuildTableFlag(G4bool val)
     610{
     611  buildLambdaTable = val;
     612  if(!val) integral = false;
     613}
     614
     615//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     616
     617inline void G4VEmProcess::SetStartFromNullFlag(G4bool val)
     618{
     619  startFromNull = val;
     620}
     621
     622//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     623
     624inline void G4VEmProcess::InitialiseStep(const G4Track& track)
     625{
     626  preStepKinEnergy = track.GetKineticEnergy();
     627  DefineMaterial(track.GetMaterialCutsCouple());
     628  if (theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX;
     629}
     630
    350631//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    351632
     
    358639    mfpKinEnergy = DBL_MAX;
    359640  }
    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));
    417641}
    418642
     
    442666//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    443667
    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;
     668inline G4double G4VEmProcess::GetLambdaFromTable(G4double e)
     669{
     670  G4bool b;
     671  return (((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b));
     672}
     673
     674//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     675
     676inline G4double G4VEmProcess::GetCurrentLambda(G4double e)
     677{
     678  G4double x = 0.0;
     679  if(theLambdaTable) x = GetLambdaFromTable(e);
     680  else               x = ComputeCurrentLambda(e);
    450681  return x;
    451682}
     
    453684//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    454685
    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;
     686inline G4double G4VEmProcess::ComputeCurrentLambda(G4double e)
     687{
     688  SelectModel(e);
     689  G4double x = 0.0;
     690  if(currentModel) {
     691    x = currentModel->CrossSectionPerVolume(currentMaterial,particle,
     692                                            e,(*theCuts)[currentMaterialIndex]);
     693  }
     694  return x;
    684695}
    685696
  • trunk/source/processes/electromagnetic/utils/include/G4VEnergyLossProcess.hh

    r991 r1005  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEnergyLossProcess.hh,v 1.83 2008/09/12 16:19:01 vnivanch Exp $
     26// $Id: G4VEnergyLossProcess.hh,v 1.86 2009/02/19 09:57:36 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 G4ParticleChangeForLoss* GetParticleChange();
     277
     278  inline size_t CurrentMaterialCutsCoupleIndex() const;
     279
     280  inline G4double GetCurrentRange() const;
     281
     282  //------------------------------------------------------------------------
     283  // Specific methods to set, access, modify models
     284  //------------------------------------------------------------------------
     285
     286  // Select model in run time
     287  inline void SelectModel(G4double kinEnergy);
    206288
    207289public:
    208 
     290  // Select model by energy and region index
     291  inline G4VEmModel* SelectModelForMaterial(G4double kinEnergy,
     292                                            size_t& idx) const;
     293
     294  // Add EM model coupled with fluctuation model for region, smaller value
     295  // of order defines which pair of models will be selected for a given
     296  // energy interval 
     297  inline void AddEmModel(G4int, G4VEmModel*,
     298                         G4VEmFluctuationModel* fluc = 0,
     299                         const G4Region* region = 0);
     300
     301  // Define new energy range for the model identified by the name
     302  inline void UpdateEmModel(const G4String&, G4double, G4double);
     303
     304  // Assign a model to a process
     305  inline void SetEmModel(G4VEmModel*, G4int index=1);
     306 
     307  // return the assigned model
     308  inline G4VEmModel* EmModel(G4int index=1);
     309 
     310  // Access to models
     311  inline G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false);
     312
     313  inline G4int NumberOfModels();
     314
     315  // Assign a fluctuation model to a process
     316  inline void SetFluctModel(G4VEmFluctuationModel*);
     317 
     318  // return the assigned fluctuation model
     319  inline G4VEmFluctuationModel* FluctModel();
     320   
     321  //------------------------------------------------------------------------
     322  // Define and access particle type
     323  //------------------------------------------------------------------------
     324
     325protected:
     326  inline void SetParticle(const G4ParticleDefinition* p);
     327  inline void SetSecondaryParticle(const G4ParticleDefinition* p);
     328
     329public:
     330  inline void SetBaseParticle(const G4ParticleDefinition* p);
     331  inline const G4ParticleDefinition* Particle() const;
     332  inline const G4ParticleDefinition* BaseParticle() const;
     333  inline const G4ParticleDefinition* SecondaryParticle() const;
     334
     335  //------------------------------------------------------------------------
     336  // Get/set parameters to configure the process at initialisation time
     337  //------------------------------------------------------------------------
     338
     339  // Add subcutoff process (bremsstrahlung) to sample secondary
     340  // particle production in vicinity of the geometry boundary
    209341  void AddCollaborativeProcess(G4VEnergyLossProcess*);
    210342
    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);
     343  inline void SetLossFluctuations(G4bool val);
     344  inline void SetRandomStep(G4bool val);
     345
     346  inline void SetIntegral(G4bool val);
     347  inline G4bool IsIntegral() const;
     348
     349  // Set/Get flag "isIonisation"
     350  inline void SetIonisation(G4bool val);
     351  inline G4bool IsIonisationProcess() const;
     352
     353  // Redefine parameteters for stepping control
     354  //
     355  inline void SetLinearLossLimit(G4double val);
     356  inline void SetMinSubRange(G4double val);
     357  inline void SetLambdaFactor(G4double val);
     358  inline void SetStepFunction(G4double v1, G4double v2);
     359
     360  inline G4int NumberOfSubCutoffRegions() const;
     361  inline G4int NumberOfDERegions() const;
     362
     363  //------------------------------------------------------------------------
     364  // Specific methods to path Physics Tables to the process
     365  //------------------------------------------------------------------------
    226366
    227367  void SetDEDXTable(G4PhysicsTable* p, G4EmTableType tType);
    228368  void SetCSDARangeTable(G4PhysicsTable* pRange);
    229369  void SetRangeTableForLoss(G4PhysicsTable* p);
     370  void SetSecondaryRangeTable(G4PhysicsTable* p);
    230371  void SetInverseRangeTable(G4PhysicsTable* p);
    231   void SetSecondaryRangeTable(G4PhysicsTable* p);
    232372
    233373  void SetLambdaTable(G4PhysicsTable* p);
     
    251391  // Max kinetic energy for tables
    252392  inline void SetMaxKinEnergyForCSDARange(G4double e);
     393
     394  // Return values for given G4MaterialCutsCouple
     395  inline G4double GetDEDX(G4double& kineticEnergy, const G4MaterialCutsCouple*);
     396  inline G4double GetDEDXForSubsec(G4double& kineticEnergy,
     397                                   const G4MaterialCutsCouple*);
     398  inline G4double GetRange(G4double& kineticEnergy, const G4MaterialCutsCouple*);
     399  inline G4double GetCSDARange(G4double& kineticEnergy, const G4MaterialCutsCouple*);
     400  inline G4double GetRangeForLoss(G4double& kineticEnergy, const G4MaterialCutsCouple*);
     401  inline G4double GetKineticEnergy(G4double& range, const G4MaterialCutsCouple*);
     402  inline G4double GetLambda(G4double& kineticEnergy, const G4MaterialCutsCouple*);
     403
     404  inline G4bool TablesAreBuilt() const;
    253405
    254406  // Access to specific tables
     
    264416  inline G4PhysicsTable* SubLambdaTable();
    265417
    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 
    359418  //------------------------------------------------------------------------
    360419  // Run time method for simulation of ionisation
     
    367426  inline void SetDynamicMassCharge(G4double massratio, G4double charge2ratio);
    368427
    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 
    390428private:
    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);
    408429
    409430  // define material and indexes
    410431  inline void DefineMaterial(const G4MaterialCutsCouple* couple);
    411432
    412   // Returnd values for scaled energy using mass of the base particle
    413   //
     433  //------------------------------------------------------------------------
     434  // Compute values using scaling relation, mass and charge of based particle
     435  //------------------------------------------------------------------------
     436
    414437  inline G4double GetDEDXForScaledEnergy(G4double scaledKinEnergy);
    415438  inline G4double GetSubDEDXForScaledEnergy(G4double scaledKinEnergy);
     
    418441  inline G4double GetScaledRangeForScaledEnergy(G4double scaledKinEnergy);
    419442  inline G4double GetLimitScaledRangeForScaledEnergy(G4double scaledKinEnergy);
     443  inline G4double ScaledKinEnergyForLoss(G4double range);
    420444  inline G4double GetLambdaForScaledEnergy(G4double scaledKinEnergy);
    421   inline G4double ScaledKinEnergyForLoss(G4double range);
    422445  inline void ComputeLambdaForScaledEnergy(G4double scaledKinEnergy);
    423446
    424447  // hide  assignment operator
    425 
    426448  G4VEnergyLossProcess(G4VEnergyLossProcess &);
    427449  G4VEnergyLossProcess & operator=(const G4VEnergyLossProcess &right);
     
    444466  G4VEmFluctuationModel*                fluctModel;
    445467  std::vector<const G4Region*>          scoffRegions;
     468  std::vector<const G4Region*>          deRegions;
    446469  G4int                                 nSCoffRegions;
    447   G4int*                                idxSCoffRegions;
     470  G4int                                 nDERegions;
     471  G4bool*                               idxSCoffRegions;
     472  G4bool*                               idxDERegions;
    448473
    449474  std::vector<G4VEnergyLossProcess*>    scProcesses;
     
    493518  G4bool   isIonisation;
    494519  G4bool   useSubCutoff;
     520  G4bool   useDeexcitation;
    495521
    496522protected:
     
    531557//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    532558
    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   }
     559inline G4ParticleChangeForLoss* G4VEnergyLossProcess::GetParticleChange()
     560{
     561  return &fParticleChange;
     562}
     563
     564//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     565
     566inline size_t G4VEnergyLossProcess::CurrentMaterialCutsCoupleIndex() const
     567{
     568  return currentMaterialIndex;
     569}
     570
     571//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     572 
     573inline G4double G4VEnergyLossProcess::GetCurrentRange() const
     574{
     575  return fRange;
     576}
     577
     578//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     579
     580inline void G4VEnergyLossProcess::SelectModel(G4double kinEnergy)
     581{
     582  currentModel = modelManager->SelectModel(kinEnergy, currentMaterialIndex);
     583  currentModel->SetCurrentCouple(currentCouple);
     584}
     585
     586//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     587
     588inline G4VEmModel* G4VEnergyLossProcess::SelectModelForMaterial(
     589                   G4double kinEnergy, size_t& idx) const
     590{
     591  return modelManager->SelectModel(kinEnergy, idx);
     592}
     593
     594//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     595
     596inline
     597void G4VEnergyLossProcess::AddEmModel(G4int order, G4VEmModel* p,
     598                                      G4VEmFluctuationModel* fluc,
     599                                      const G4Region* region)
     600{
     601  modelManager->AddEmModel(order, p, fluc, region);
     602  if(p) p->SetParticleChange(pParticleChange, fluc);
     603}
     604
     605//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     606
     607inline void G4VEnergyLossProcess::UpdateEmModel(const G4String& nam,
     608                                                G4double emin, G4double emax)
     609{
     610  modelManager->UpdateEmModel(nam, emin, emax);
     611}
     612
     613//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     614
     615inline void G4VEnergyLossProcess::SetEmModel(G4VEmModel* p, G4int index)
     616{
     617  G4int n = emModels.size();
     618  if(index >= n) for(G4int i=n; i<index+1; i++) {emModels.push_back(0);}
     619  emModels[index] = p;
     620}
     621
     622//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     623
     624inline G4VEmModel* G4VEnergyLossProcess::EmModel(G4int index)
     625{
     626  G4VEmModel* p = 0;
     627  if(index >= 0 && index <  G4int(emModels.size())) p = emModels[index];
     628  return p;
     629}
     630
     631//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     632
     633inline
     634G4VEmModel* G4VEnergyLossProcess::GetModelByIndex(G4int idx, G4bool ver)
     635{
     636  return modelManager->GetModel(idx, ver);
     637}
     638
     639//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     640
     641inline G4int G4VEnergyLossProcess::NumberOfModels()
     642{
     643  return modelManager->NumberOfModels();
     644}
     645
     646//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     647
     648inline void G4VEnergyLossProcess::SetFluctModel(G4VEmFluctuationModel* p)
     649{
     650  fluctModel = p;
     651}
     652
     653//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     654
     655inline G4VEmFluctuationModel* G4VEnergyLossProcess::FluctModel()
     656{
     657  return fluctModel;
     658}
     659
     660//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     661
     662inline void G4VEnergyLossProcess::SetParticle(const G4ParticleDefinition* p)
     663{
     664  particle = p;
     665}
     666
     667//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     668
     669inline void G4VEnergyLossProcess::SetSecondaryParticle(const G4ParticleDefinition* p)
     670{
     671  secondaryParticle = p;
     672}
     673
     674//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     675
     676inline void G4VEnergyLossProcess::SetBaseParticle(const G4ParticleDefinition* p)
     677{
     678  baseParticle = p;
     679}
     680
     681//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     682
     683inline const G4ParticleDefinition* G4VEnergyLossProcess::Particle() const
     684{
     685  return particle;
     686}
     687
     688//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     689
     690inline const G4ParticleDefinition* G4VEnergyLossProcess::BaseParticle() const
     691{
     692  return baseParticle;
     693}
     694
     695//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     696
     697inline const G4ParticleDefinition* G4VEnergyLossProcess::SecondaryParticle() const
     698{
     699  return secondaryParticle;
     700}
     701
     702//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     703
     704inline void G4VEnergyLossProcess::SetLossFluctuations(G4bool val)
     705{
     706  lossFluctuationFlag = val;
     707}
     708
     709//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     710
     711inline void G4VEnergyLossProcess::SetRandomStep(G4bool val)
     712{
     713  rndmStepFlag = val;
     714}
     715
     716//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     717
     718inline void G4VEnergyLossProcess::SetIntegral(G4bool val)
     719{
     720  integral = val;
     721}
     722
     723//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     724 
     725inline G4bool G4VEnergyLossProcess::IsIntegral() const
     726{
     727  return integral;
     728}
     729
     730//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     731
     732inline void G4VEnergyLossProcess::SetIonisation(G4bool val)
     733{
     734  isIonisation = val;
     735  if(val) aGPILSelection = CandidateForSelection;
     736  else    aGPILSelection = NotCandidateForSelection;
     737}
     738
     739//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     740
     741inline G4bool G4VEnergyLossProcess::IsIonisationProcess() const
     742{
     743  return isIonisation;
     744}
     745
     746//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     747
     748inline void G4VEnergyLossProcess::SetLinearLossLimit(G4double val)
     749{
     750  linLossLimit = val;
     751}
     752
     753//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     754
     755inline void G4VEnergyLossProcess::SetMinSubRange(G4double val)
     756{
     757  minSubRange = val;
     758}
     759
     760//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     761
     762inline void G4VEnergyLossProcess::SetLambdaFactor(G4double val)
     763{
     764  if(val > 0.0 && val <= 1.0) lambdaFactor = val;
     765}
     766
     767//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     768
     769void G4VEnergyLossProcess::SetStepFunction(G4double v1, G4double v2)
     770{
     771  dRoverRange = v1;
     772  finalRange = v2;
     773  if (dRoverRange > 0.999) dRoverRange = 1.0;
     774  currentCouple = 0;
     775  mfpKinEnergy  = DBL_MAX;
     776}
     777
     778//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     779
     780inline G4int G4VEnergyLossProcess::NumberOfSubCutoffRegions() const
     781{
     782  return nSCoffRegions;
     783}
     784
     785//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     786
     787inline G4int G4VEnergyLossProcess::NumberOfDERegions() const
     788{
     789  return nDERegions;
     790}
     791
     792//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     793
     794inline void G4VEnergyLossProcess::SetDEDXBinning(G4int nbins)
     795{
     796  nBins = nbins;
     797}
     798
     799//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     800
     801inline void G4VEnergyLossProcess::SetLambdaBinning(G4int nbins)
     802{
     803  nBins = nbins;
     804}
     805
     806//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     807
     808inline void G4VEnergyLossProcess::SetDEDXBinningForCSDARange(G4int nbins)
     809{
     810  nBinsCSDA = nbins;
     811}
     812
     813//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     814
     815inline void G4VEnergyLossProcess::SetMinKinEnergy(G4double e)
     816{
     817  minKinEnergy = e;
     818}
     819
     820//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     821
     822inline G4double G4VEnergyLossProcess::MinKinEnergy() const
     823{
     824  return minKinEnergy;
     825}
     826
     827//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     828
     829inline void G4VEnergyLossProcess::SetMaxKinEnergy(G4double e)
     830{
     831  maxKinEnergy = e;
     832  if(e < maxKinEnergyCSDA) maxKinEnergyCSDA = e;
     833}
     834
     835//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     836
     837inline G4double G4VEnergyLossProcess::MaxKinEnergy() const
     838{
     839  return maxKinEnergy;
     840}
     841
     842//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     843
     844inline void G4VEnergyLossProcess::SetMaxKinEnergyForCSDARange(G4double e)
     845{
     846  maxKinEnergyCSDA = e;
    542847}
    543848
     
    558863  DefineMaterial(couple);
    559864  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;
    611865}
    612866
     
    643897//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    644898
     899inline G4double G4VEnergyLossProcess::GetRangeForLoss(
     900                G4double& kineticEnergy,
     901                const G4MaterialCutsCouple* couple)
     902{
     903  DefineMaterial(couple);
     904  G4double x = DBL_MAX;
     905  if(theRangeTableForLoss)
     906    x = GetScaledRangeForScaledEnergy(kineticEnergy*massRatio)*reduceFactor;
     907  //  G4cout << "Range from " << GetProcessName()
     908  //         << "  e= " << kineticEnergy << " r= " << x << G4endl;
     909  return x;
     910}
     911
     912//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     913
     914inline G4double G4VEnergyLossProcess::GetKineticEnergy(
     915                G4double& range,
     916                const G4MaterialCutsCouple* couple)
     917{
     918  DefineMaterial(couple);
     919  G4double r = range/reduceFactor;
     920  G4double e = ScaledKinEnergyForLoss(r)/massRatio;
     921  return e;
     922}
     923
     924//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     925
     926inline G4double G4VEnergyLossProcess::GetLambda(G4double& kineticEnergy,
     927                                          const G4MaterialCutsCouple* couple)
     928{
     929  DefineMaterial(couple);
     930  G4double x = 0.0;
     931  if(theLambdaTable) x = GetLambdaForScaledEnergy(kineticEnergy*massRatio);
     932  return x;
     933}
     934
     935//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     936
     937inline G4bool G4VEnergyLossProcess::TablesAreBuilt() const
     938{
     939  return  tablesAreBuilt;
     940}
     941
     942//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     943
     944inline G4PhysicsTable* G4VEnergyLossProcess::DEDXTable() const
     945{
     946  return theDEDXTable;
     947}
     948
     949//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     950
     951inline G4PhysicsTable* G4VEnergyLossProcess::DEDXTableForSubsec() const
     952{
     953  return theDEDXSubTable;
     954}
     955
     956//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     957
     958inline G4PhysicsTable* G4VEnergyLossProcess::DEDXunRestrictedTable() const
     959{
     960  return theDEDXunRestrictedTable;
     961}
     962
     963//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     964
     965inline G4PhysicsTable* G4VEnergyLossProcess::IonisationTable() const
     966{
     967  G4PhysicsTable* t = theDEDXTable;
     968  if(theIonisationTable) t = theIonisationTable;
     969  return t;
     970}
     971
     972//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     973
     974inline G4PhysicsTable* G4VEnergyLossProcess::IonisationTableForSubsec() const
     975{
     976  G4PhysicsTable* t = theDEDXSubTable;
     977  if(theIonisationSubTable) t = theIonisationSubTable;
     978  return t;
     979}
     980
     981//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     982
     983inline G4PhysicsTable* G4VEnergyLossProcess::CSDARangeTable() const
     984{
     985  return theCSDARangeTable;
     986}
     987
     988//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     989
     990inline G4PhysicsTable* G4VEnergyLossProcess::RangeTableForLoss() const
     991{
     992  return theRangeTableForLoss;
     993}
     994
     995//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     996
     997inline G4PhysicsTable* G4VEnergyLossProcess::InverseRangeTable() const
     998{
     999  return theInverseRangeTable;
     1000}
     1001
     1002//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1003
     1004inline G4PhysicsTable* G4VEnergyLossProcess::LambdaTable()
     1005{
     1006  return theLambdaTable;
     1007}
     1008
     1009//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1010
     1011inline G4PhysicsTable* G4VEnergyLossProcess::SubLambdaTable()
     1012{
     1013  return theSubLambdaTable;
     1014}
     1015
     1016//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1017
     1018inline G4double G4VEnergyLossProcess::SampleRange()
     1019{
     1020  G4double e = amu_c2*preStepKinEnergy/particle->GetPDGMass();
     1021  G4bool b;
     1022  G4double s = fRange*std::pow(10.,vstrag->GetValue(e,b));
     1023  G4double x = fRange + G4RandGauss::shoot(0.0,s);
     1024  if(x > 0.0) fRange = x;
     1025  return fRange;
     1026}
     1027
     1028//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1029
     1030inline void G4VEnergyLossProcess::SetDynamicMassCharge(G4double massratio,
     1031                                                       G4double charge2ratio)
     1032{
     1033  massRatio     = massratio;
     1034  chargeSqRatio = charge2ratio;
     1035  reduceFactor  = 1.0/(chargeSqRatio*massRatio);
     1036}
     1037
     1038//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1039
     1040inline void G4VEnergyLossProcess::DefineMaterial(
     1041            const G4MaterialCutsCouple* couple)
     1042{
     1043  if(couple != currentCouple) {
     1044    currentCouple   = couple;
     1045    currentMaterial = couple->GetMaterial();
     1046    currentMaterialIndex = couple->GetIndex();
     1047    mfpKinEnergy = DBL_MAX;
     1048  }
     1049}
     1050
     1051//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1052
     1053inline G4double G4VEnergyLossProcess::GetDEDXForScaledEnergy(G4double e)
     1054{
     1055  G4bool b;
     1056  G4double x =
     1057    ((*theDEDXTable)[currentMaterialIndex]->GetValue(e, b))*chargeSqRatio;
     1058  if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
     1059  return x;
     1060}
     1061
     1062//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1063
     1064inline G4double G4VEnergyLossProcess::GetSubDEDXForScaledEnergy(G4double e)
     1065{
     1066  G4bool b;
     1067  G4double x =
     1068    ((*theDEDXSubTable)[currentMaterialIndex]->GetValue(e, b))*chargeSqRatio;
     1069  if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
     1070  return x;
     1071}
     1072
     1073//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1074
     1075inline G4double G4VEnergyLossProcess::GetIonisationForScaledEnergy(G4double e)
     1076{
     1077  G4bool b;
     1078  G4double x = 0.0;
     1079  //  if(theIonisationTable) {
     1080  x = ((*theIonisationTable)[currentMaterialIndex]->GetValue(e, b))
     1081    *chargeSqRatio;
     1082  if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
     1083  //}
     1084  return x;
     1085}
     1086
     1087//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1088
     1089inline
     1090G4double G4VEnergyLossProcess::GetSubIonisationForScaledEnergy(G4double e)
     1091{
     1092  G4bool b;
     1093  G4double x = 0.0;
     1094  //if(theIonisationSubTable) {
     1095  x = ((*theIonisationSubTable)[currentMaterialIndex]->GetValue(e, b))
     1096    *chargeSqRatio;
     1097  if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
     1098  //}
     1099  return x;
     1100}
     1101
     1102//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1103
     1104inline G4double G4VEnergyLossProcess::GetScaledRangeForScaledEnergy(G4double e)
     1105{
     1106  G4bool b;
     1107  G4double x = ((*theRangeTableForLoss)[currentMaterialIndex])->GetValue(e, b);
     1108  if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy);
     1109  return x;
     1110}
     1111
     1112//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1113
    6451114inline G4double G4VEnergyLossProcess::GetLimitScaledRangeForScaledEnergy(
    6461115                G4double e)
     
    6571126  }
    6581127  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;
    6961128}
    6971129
     
    7111143  }
    7121144  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;
    7241145}
    7251146
     
    7581179//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    7591180
    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 
    11751181#endif
  • trunk/source/processes/electromagnetic/utils/include/G4VMscModel.hh

    r991 r1005  
    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.8 2009/02/24 09:56:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
     
    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  // empty
     72  virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
     73                                 const G4MaterialCutsCouple*,
     74                                 const G4DynamicParticle*,
     75                                 G4double tmin,
     76                                 G4double tmax);
     77
     78  //================================================================
     79  //  Set parameters of multiple scattering models
     80  //================================================================
     81 
    6682  inline void SetStepLimitType(G4MscStepLimitType);
    6783
     
    7389
    7490  inline void SetSkin(G4double);
     91
     92  inline void SetSampleZ(G4bool);
     93
     94protected:
     95
     96  // initialisation of interface with geometry
     97  void InitialiseSafetyHelper();
     98
     99  // shift point of the track PostStep
     100  void ComputeDisplacement(G4ParticleChangeForMSC*, 
     101                           const G4ThreeVector& displDir,
     102                           G4double displacement,
     103                           G4double postsafety);
     104
     105  // compute safety
     106  inline G4double ComputeSafety(const G4ThreeVector& position, G4double limit);
     107
     108  // compute linear distance to a geometry boundary
     109  inline G4double ComputeGeomLimit(const G4Track& position, G4double& presafety,
     110                                   G4double limit);
    75111
    76112private:
     
    79115  G4VMscModel & operator=(const  G4VMscModel &right);
    80116  G4VMscModel(const  G4VMscModel&);
     117
     118  G4SafetyHelper* safetyHelper;
    81119
    82120protected:
     
    88126  G4double dtrl;
    89127  G4double lambdalimit;
     128  G4double geommax;
    90129
    91130  G4MscStepLimitType steppingAlgorithm;
     
    134173//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    135174
     175inline void G4VMscModel::SetSampleZ(G4bool val)
     176{
     177  samplez = val;
     178}
     179
     180//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     181
     182inline G4double G4VMscModel::ComputeSafety(const G4ThreeVector& position,
     183                                           G4double)
     184{
     185  return safetyHelper->ComputeSafety(position);
     186}
     187
     188//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     189
     190inline G4double G4VMscModel::ComputeGeomLimit(const G4Track& track,
     191                                              G4double& presafety,
     192                                              G4double limit)
     193{
     194  G4double res = geommax;
     195  if(track.GetVolume() != safetyHelper->GetWorldVolume()) {
     196    res = safetyHelper->CheckNextStep(
     197          track.GetStep()->GetPreStepPoint()->GetPosition(),
     198          track.GetMomentumDirection(),
     199          limit, presafety);
     200  }
     201  return res;
     202}
     203
     204//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     205
    136206#endif
    137207
  • trunk/source/processes/electromagnetic/utils/include/G4VMultipleScattering.hh

    r991 r1005  
    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.55 2009/02/18 12:19:33 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
     
    145145                              G4bool ascii);
    146146
    147   //------------------------------------------------------------------------
    148   // Specific methods for msc processes
    149   //------------------------------------------------------------------------
    150 
    151147  // The function overloads the corresponding function of the base
    152148  // class.It limits the step near to boundaries only
    153149  // and invokes the method GetMscContinuousStepLimit at every step.
    154   virtual G4double AlongStepGetPhysicalInteractionLength(
     150  G4double AlongStepGetPhysicalInteractionLength(
    155151                                            const G4Track&,
    156152                                            G4double  previousStepSize,
     
    192188  inline G4PhysicsTable* LambdaTable() const;
    193189
    194   //------------------------------------------------------------------------
    195   // Define and access particle type
    196   //------------------------------------------------------------------------
    197 
     190  // access particle type
    198191  inline const G4ParticleDefinition* Particle() const;
    199   inline void SetParticle(const G4ParticleDefinition*);
    200192
    201193  //------------------------------------------------------------------------
     
    203195  //------------------------------------------------------------------------
    204196
    205   inline void AddEmModel(G4int, G4VEmModel*, const G4Region* region = 0);
    206 
     197protected:
     198  // Select model in run time
     199  inline G4VEmModel* SelectModel(G4double kinEnergy);
     200
     201public:
     202  // Select model in run time
    207203  inline G4VEmModel* SelectModelForMaterial(G4double kinEnergy,
    208204                                            size_t& idxRegion) const;
    209205
    210   // Access to models
     206  // Add model for region, smaller value of order defines which
     207  // model will be selected for a given energy interval 
     208  inline void AddEmModel(G4int order, G4VEmModel*, const G4Region* region = 0);
     209
     210  // Access to models by index
    211211  inline G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false);
    212212
    213213  //------------------------------------------------------------------------
    214   // Set parameters for simulation of multiple scattering
    215   //------------------------------------------------------------------------
    216 
     214  // Get/Set parameters for simulation of multiple scattering
     215  //------------------------------------------------------------------------
     216
     217  inline G4bool LateralDisplasmentFlag() const;
    217218  inline void SetLateralDisplasmentFlag(G4bool val);
    218219
     220  inline G4double Skin() const;
    219221  inline void SetSkin(G4double val);
    220222
     223  inline G4double RangeFactor() const;
    221224  inline void SetRangeFactor(G4double val);
    222225
     226  inline G4double GeomFactor() const;
    223227  inline void SetGeomFactor(G4double val);
    224228
     229  inline G4double PolarAngleLimit() const;
    225230  inline void SetPolarAngleLimit(G4double val);
    226231
     232  inline G4MscStepLimitType StepLimitType() const;
    227233  inline void SetStepLimitType(G4MscStepLimitType val);
     234
     235  //------------------------------------------------------------------------
     236  // Run time methods
     237  //------------------------------------------------------------------------
    228238
    229239protected:
     
    233243                           G4double,
    234244                           G4ForceCondition* condition);
    235 
    236   //------------------------------------------------------------------------
    237   // Run time methods
    238   //------------------------------------------------------------------------
    239245
    240246  // This method is not used for tracking, it returns step limit
     
    253259                                            G4double& currentSafety);
    254260
    255   inline G4VEmModel* SelectModel(G4double kinEnergy);
    256   // Select concrete model
    257 
    258   inline const G4MaterialCutsCouple* CurrentMaterialCutsCouple() const;
    259 
    260261  // define current material
    261262  inline void DefineMaterial(const G4MaterialCutsCouple* couple);
    262263
    263   //------------------------------------------------------------------------
    264   // Access parameters of multiple scattering
    265   //------------------------------------------------------------------------
     264  inline const G4MaterialCutsCouple* CurrentMaterialCutsCouple() const;
    266265
    267266  inline G4ParticleChangeForMSC* GetParticleChange();
    268267
    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;
    280268
    281269private:
    282270
    283271  // hide  assignment operator
    284 
    285272  G4VMultipleScattering(G4VMultipleScattering &);
    286273  G4VMultipleScattering & operator=(const G4VMultipleScattering &right);
     
    330317//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    331318
     319inline G4double G4VMultipleScattering::ContinuousStepLimit(
     320                                       const G4Track& track,
     321                                       G4double previousStepSize,
     322                                       G4double currentMinimalStep,
     323                                       G4double& currentSafety)
     324{
     325  return GetMscContinuousStepLimit(track,previousStepSize,currentMinimalStep,
     326                                   currentSafety);
     327}
     328
     329//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     330
     331inline void G4VMultipleScattering::SetBinning(G4int nbins)
     332{
     333  nBins = nbins;
     334}
     335
     336//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     337
     338inline G4int G4VMultipleScattering::Binning() const
     339{
     340  return nBins;
     341}
     342
     343//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     344
     345inline void G4VMultipleScattering::SetMinKinEnergy(G4double e)
     346{
     347  minKinEnergy = e;
     348}
     349
     350//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     351
     352inline G4double G4VMultipleScattering::MinKinEnergy() const
     353{
     354  return minKinEnergy;
     355}
     356
     357//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     358
     359inline void G4VMultipleScattering::SetMaxKinEnergy(G4double e)
     360{
     361  maxKinEnergy = e;
     362}
     363
     364//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     365
     366inline G4double G4VMultipleScattering::MaxKinEnergy() const
     367{
     368  return maxKinEnergy;
     369}
     370
     371//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     372
     373inline  void G4VMultipleScattering::SetBuildLambdaTable(G4bool val)
     374{
     375  buildLambdaTable = val;
     376}
     377
     378//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     379
     380inline G4PhysicsTable* G4VMultipleScattering::LambdaTable() const
     381{
     382  return theLambdaTable;
     383}
     384
     385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     386
     387inline  const G4ParticleDefinition* G4VMultipleScattering::Particle() const
     388{
     389  return currentParticle;
     390}
     391
     392//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     393
     394inline void G4VMultipleScattering::AddEmModel(G4int order, G4VEmModel* p,
     395                                              const G4Region* region)
     396{
     397  G4VEmFluctuationModel* fm = 0;
     398  modelManager->AddEmModel(order, p, fm, region);
     399  if(p) p->SetParticleChange(pParticleChange);
     400}
     401
     402//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     403
     404inline G4VEmModel* G4VMultipleScattering::SelectModel(G4double kinEnergy)
     405{
     406  return modelManager->SelectModel(kinEnergy, currentMaterialIndex);
     407}
     408
     409//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     410
     411inline G4VEmModel* G4VMultipleScattering::SelectModelForMaterial(
     412                   G4double kinEnergy, size_t& idxRegion) const
     413{
     414  return modelManager->SelectModel(kinEnergy, idxRegion);
     415}
     416
     417//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     418
    332419inline
    333 void G4VMultipleScattering::DefineMaterial(const G4MaterialCutsCouple* couple)
    334 {
    335   if(couple != currentCouple) {
    336     currentCouple   = couple;
    337     currentMaterialIndex = couple->GetIndex();
     420G4VEmModel* G4VMultipleScattering::GetModelByIndex(G4int idx, G4bool ver)
     421{
     422  return modelManager->GetModel(idx, ver);
     423}
     424
     425//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     426
     427inline  G4bool G4VMultipleScattering::LateralDisplasmentFlag() const
     428{
     429  return latDisplasment;
     430}
     431
     432//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     433
     434inline  void G4VMultipleScattering::SetLateralDisplasmentFlag(G4bool val)
     435{
     436  latDisplasment = val;
     437}
     438
     439//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     440
     441inline  G4double G4VMultipleScattering::Skin() const
     442{
     443  return skin;
     444}
     445
     446//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     447
     448inline  void G4VMultipleScattering::SetSkin(G4double val)
     449{
     450  if(val < 1.0) skin = 0.0;
     451  else          skin = val;
     452}
     453
     454//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     455
     456inline  G4double G4VMultipleScattering::RangeFactor() const
     457{
     458  return facrange;
     459}
     460
     461//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     462
     463inline  void G4VMultipleScattering::SetRangeFactor(G4double val)
     464{
     465  if(val > 0.0) facrange = val;
     466}
     467
     468//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     469
     470inline  G4double G4VMultipleScattering::GeomFactor() const
     471{
     472  return facgeom;
     473}
     474
     475//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     476
     477inline  void G4VMultipleScattering::SetGeomFactor(G4double val)
     478{
     479  if(val > 0.0) facgeom = val;
     480}
     481
     482//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     483
     484inline  G4double G4VMultipleScattering::PolarAngleLimit() const
     485{
     486  return polarAngleLimit;
     487}
     488
     489//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     490
     491inline  void G4VMultipleScattering::SetPolarAngleLimit(G4double val)
     492{
     493  if(val < 0.0)     polarAngleLimit = 0.0;
     494  else if(val > pi) polarAngleLimit = pi;
     495  else              polarAngleLimit = val;
     496}
     497
     498//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     499
     500inline G4MscStepLimitType G4VMultipleScattering::StepLimitType() const
     501{
     502  return stepLimit;
     503}
     504
     505//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     506
     507inline void G4VMultipleScattering::SetStepLimitType(G4MscStepLimitType val)
     508{
     509  stepLimit = val;
     510  if(val == fMinimal) facrange = 0.2;
     511}
     512
     513//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     514
     515inline
     516G4double G4VMultipleScattering::GetLambda(const G4ParticleDefinition* p,
     517                                          G4double& e)
     518{
     519  G4double x;
     520  if(theLambdaTable) {
     521    G4bool b;
     522    x = ((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b);
     523  } else {
     524    x = currentModel->CrossSection(currentCouple,p,e);
    338525  }
     526  if(x > DBL_MIN) x = 1./x;
     527  else            x = DBL_MAX;
     528  return x;
    339529}
    340530
     
    364554//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    365555
    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 
    378556inline
    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);
     557void G4VMultipleScattering::DefineMaterial(const G4MaterialCutsCouple* couple)
     558{
     559  if(couple != currentCouple) {
     560    currentCouple   = couple;
     561    currentMaterialIndex = couple->GetIndex();
    388562  }
    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
     563}
     564
     565//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     566
     567inline const G4MaterialCutsCouple*
     568G4VMultipleScattering::CurrentMaterialCutsCouple() const
    571569{
    572570  return currentCouple;
     
    575573//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    576574
    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);
     575inline G4ParticleChangeForMSC* G4VMultipleScattering::GetParticleChange()
     576{
     577  return &fParticleChange;
    591578}
    592579
  • trunk/source/processes/electromagnetic/utils/src/G4EmCalculator.cc

    r991 r1005  
    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-02-ref-02 $
    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

    r991 r1005  
    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-02-ref-02 $
    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/G4EmProcessOptions.cc

    r991 r1005  
    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-02-ref-02 $
    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

    r991 r1005  
    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-02-ref-02 $
    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

    r991 r1005  
    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-02-ref-02 $
    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/G4VEmFluctuationModel.cc

    r991 r1005  
    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-02-ref-02 $
    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

    r991 r1005  
    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.25 2009/02/19 09:57:36 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    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//
     
    6061  fluc(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV),
    6162  polarAngleLimit(0.0),secondaryThreshold(DBL_MAX),theLPMflag(false),
    62   pParticleChange(0),nuclearStopping(false),nsec(5)
     63  pParticleChange(0),nuclearStopping(false),
     64  currentCouple(0),currentElement(0),
     65  nsec(5),flagDeexcitation(false)
    6366{
    6467  xsec.resize(nsec);
     
    8285//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    8386
    84 G4double G4VEmModel::CrossSectionPerVolume(const G4Material* material,
    85                                            const G4ParticleDefinition* p,
    86                                            G4double ekin,
    87                                            G4double emin,
    88                                            G4double emax)
    89 {
    90   SetupForMaterial(p, material, ekin);
    91   G4double cross = 0.0;
    92   const G4ElementVector* theElementVector = material->GetElementVector();
    93   const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
    94   G4int nelm = material->GetNumberOfElements();
    95   if(nelm > nsec) {
    96     xsec.resize(nelm);
    97     nsec = nelm;
    98   }
    99   for (G4int i=0; i<nelm; i++) {
    100     cross += theAtomNumDensityVector[i]*
    101       ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax);
    102     xsec[i] = cross;
    103   }
    104   return cross;
    105 }
    106 
    107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    108 
    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 
    12387void G4VEmModel::InitialiseElementSelectors(const G4ParticleDefinition* p,
    12488                                            const G4DataVector& cuts)
    12589{
     90  // initialise before run
     91  flagDeexcitation = false;
     92
    12693  G4int nbins = G4int(std::log10(highLimit/lowLimit) + 0.5);
    12794  if(nbins < 3) nbins = 3;
     
    162129//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    163130
    164 
     131G4double G4VEmModel::ComputeDEDXPerVolume(const G4Material*,
     132                                          const G4ParticleDefinition*,
     133                                          G4double,G4double)
     134{
     135  return 0.0;
     136}
     137
     138//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     139
     140G4double G4VEmModel::CrossSectionPerVolume(const G4Material* material,
     141                                           const G4ParticleDefinition* p,
     142                                           G4double ekin,
     143                                           G4double emin,
     144                                           G4double emax)
     145{
     146  SetupForMaterial(p, material, ekin);
     147  G4double cross = 0.0;
     148  const G4ElementVector* theElementVector = material->GetElementVector();
     149  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
     150  G4int nelm = material->GetNumberOfElements();
     151  if(nelm > nsec) {
     152    xsec.resize(nelm);
     153    nsec = nelm;
     154  }
     155  for (G4int i=0; i<nelm; i++) {
     156    cross += theAtomNumDensityVector[i]*
     157      ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax);
     158    xsec[i] = cross;
     159  }
     160  return cross;
     161}
     162
     163//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     164
     165G4double G4VEmModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
     166                                                G4double, G4double, G4double,
     167                                                G4double, G4double)
     168{
     169  return 0.0;
     170}
     171
     172//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     173
     174G4double G4VEmModel::MinEnergyCut(const G4ParticleDefinition*,
     175                                  const G4MaterialCutsCouple*)
     176{
     177  return 0.0;
     178}
     179
     180//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     181
     182G4double G4VEmModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
     183                                          const G4Material*, G4double)
     184{
     185  G4double q = p->GetPDGCharge()/CLHEP::eplus;
     186  return q*q;
     187}
     188
     189//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     190
     191G4double G4VEmModel::GetParticleCharge(const G4ParticleDefinition* p,
     192                                       const G4Material*, G4double)
     193{
     194  return p->GetPDGCharge();
     195}
     196
     197//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     198
     199void G4VEmModel::CorrectionsAlongStep(const G4MaterialCutsCouple*,
     200                                      const G4DynamicParticle*,
     201                                      G4double&,G4double&,G4double)
     202{}
     203
     204//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     205
     206void G4VEmModel::SampleDeexcitationAlongStep(const G4Material*,
     207                                             const G4Track&,
     208                                             G4double& )
     209{}
     210
     211//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     212
     213G4double G4VEmModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
     214                                        G4double kineticEnergy)
     215{
     216  return kineticEnergy;
     217}
     218
     219//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     220
     221void G4VEmModel::SampleScattering(const G4DynamicParticle*, G4double)
     222{}
     223
     224//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     225
     226G4double G4VEmModel::ComputeTruePathLengthLimit(const G4Track&,
     227                                                G4PhysicsTable*,
     228                                                G4double)
     229{
     230  return DBL_MAX;
     231}
     232
     233//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     234
     235G4double G4VEmModel::ComputeGeomPathLength(G4double truePathLength)
     236{
     237  return truePathLength;
     238}
     239
     240//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     241
     242G4double G4VEmModel::ComputeTrueStepLength(G4double geomPathLength)
     243{
     244  return geomPathLength;
     245}
     246
     247//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     248
     249void G4VEmModel::DefineForRegion(const G4Region*)
     250{}
     251
     252//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     253
     254void G4VEmModel::SetupForMaterial(const G4ParticleDefinition*,
     255                                  const G4Material*, G4double)
     256{}
     257
     258//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  • trunk/source/processes/electromagnetic/utils/src/G4VEmProcess.cc

    r991 r1005  
    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.62 2009/02/19 09:57:36 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
     
    9191  applyCuts(false),
    9292  startFromNull(true),
    93   nRegions(0),
    94   selectedModel(0),
     93  useDeexcitation(false),
     94  nDERegions(0),
     95  idxDERegions(0),
     96  currentModel(0),
    9597  particle(0),
    9698  currentCouple(0)
     
    135137  delete modelManager;
    136138  (G4LossTableManager::Instance())->DeRegister(this);
     139}
     140
     141//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     142
     143void G4VEmProcess::Clear()
     144{
     145  delete [] theEnergyOfCrossSectionMax;
     146  delete [] theCrossSectionMax;
     147  delete [] idxDERegions;
     148  theEnergyOfCrossSectionMax = 0;
     149  theCrossSectionMax = 0;
     150  idxDERegions = 0;
     151  currentCouple = 0;
     152  preStepLambda = 0.0;
     153  mfpKinEnergy  = DBL_MAX;
     154  deRegions.clear();
     155  nDERegions = 0;
    137156}
    138157
     
    162181      theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
    163182  }
    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;
     183  // Sub Cutoff and Deexcitation
     184  if (nDERegions>0) {
     185
     186    const G4ProductionCutsTable* theCoupleTable=
     187          G4ProductionCutsTable::GetProductionCutsTable();
     188    size_t numOfCouples = theCoupleTable->GetTableSize();
     189
     190    idxDERegions = new G4bool[numOfCouples];
     191
     192    for (size_t j=0; j<numOfCouples; j++) {
     193
     194      const G4MaterialCutsCouple* couple =
     195        theCoupleTable->GetMaterialCutsCouple(j);
     196      const G4ProductionCuts* pcuts = couple->GetProductionCuts();
     197      G4bool reg = false;
     198      for(G4int i=0; i<nDERegions; i++) {
     199        if(deRegions[i]) {
     200          if(pcuts == deRegions[i]->GetProductionCuts()) reg = true;
     201        }
     202      }
     203      idxDERegions[j] = reg;
     204    }
     205  }
     206  if (1 < verboseLevel && nDERegions>0) {
     207    G4cout << " Deexcitation is activated for regions: " << G4endl;
     208    for (G4int i=0; i<nDERegions; i++) {
     209      const G4Region* r = deRegions[i];
     210      G4cout << "           " << r->GetName() << G4endl;
     211    }
     212  }
    177213}
    178214
     
    237273      G4cout << *theLambdaTable << G4endl;
    238274    }
     275  }
     276}
     277
     278//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     279
     280void G4VEmProcess::PrintInfoDefinition()
     281{
     282  if(verboseLevel > 0) {
     283    G4cout << G4endl << GetProcessName() << ":   for  "
     284           << particle->GetParticleName();
     285    if(integral) G4cout << ", integral: 1 ";
     286    if(applyCuts) G4cout << ", applyCuts: 1 ";
     287    G4cout << "    SubType= " << GetProcessSubType() << G4endl;
     288    if(buildLambdaTable) {
     289      G4cout << "      Lambda tables from "
     290             << G4BestUnit(minKinEnergy,"Energy")
     291             << " to "
     292             << G4BestUnit(maxKinEnergy,"Energy")
     293             << " in " << nLambdaBins << " bins, spline: "
     294             << (G4LossTableManager::Instance())->SplineFlag()
     295             << G4endl;
     296    }
     297    PrintInfo();
     298    modelManager->DumpModelList(verboseLevel);
     299  }
     300
     301  if(verboseLevel > 2 && buildLambdaTable) {
     302    G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
     303    if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
    239304  }
    240305}
     
    304369//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    305370
    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 
    316371G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
    317372                                              const G4Step&)
     
    342397  }
    343398
    344   G4VEmModel* currentModel = SelectModel(finalT);
    345 
     399  SelectModel(finalT);
     400  if(useDeexcitation) {
     401    currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
     402  }
    346403  /* 
    347404  if(0 < verboseLevel) {
     
    404461//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    405462
    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 
    456463G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
    457464                                       const G4String& directory,
     
    521528
    522529  return yes;
     530}
     531
     532//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     533
     534void G4VEmProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
     535{
     536  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     537  const G4Region* reg = r;
     538  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
     539
     540  // the region is in the list
     541  if (nDERegions) {
     542    for (G4int i=0; i<nDERegions; i++) {
     543      if (reg == deRegions[i]) {
     544        if(!val) deRegions[i] = 0;
     545        return;
     546      }
     547    }
     548  }
     549
     550  // new region
     551  if(val) {
     552    useDeexcitation = true;
     553    deRegions.push_back(reg);
     554    nDERegions++;
     555  } else {
     556    useDeexcitation = false;
     557  }
     558}
     559
     560//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     561
     562G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
     563                                             const G4MaterialCutsCouple* couple)
     564{
     565  // Cross section per atom is calculated
     566  DefineMaterial(couple);
     567  G4double cross = 0.0;
     568  G4bool b;
     569  if(theLambdaTable) {
     570    cross = (((*theLambdaTable)[currentMaterialIndex])->
     571                           GetValue(kineticEnergy, b));
     572  } else {
     573    SelectModel(kineticEnergy);
     574    cross = currentModel->CrossSectionPerVolume(currentMaterial,
     575                                                particle,kineticEnergy);
     576  }
     577
     578  return cross;
     579}
     580
     581//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     582
     583G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
     584                                       G4double,
     585                                       G4ForceCondition* condition)
     586{
     587  *condition = NotForced;
     588  return G4VEmProcess::MeanFreePath(track);
    523589}
    524590
  • trunk/source/processes/electromagnetic/utils/src/G4VEnergyLossProcess.cc

    r991 r1005  
    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.146 2009/02/19 11:25:50 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
     
    150150  secondaryParticle(0),
    151151  nSCoffRegions(0),
     152  nDERegions(0),
    152153  idxSCoffRegions(0),
     154  idxDERegions(0),
    153155  nProcesses(0),
    154156  theDEDXTable(0),
     
    176178  isIonisation(true),
    177179  useSubCutoff(false),
     180  useDeexcitation(false),
    178181  particle(0),
    179182  currentCouple(0),
     
    237240           << G4endl;
    238241  delete vstrag;
    239   Clear();
     242  Clean();
    240243
    241244  if ( !baseParticle ) {
     
    291294//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    292295
    293 void G4VEnergyLossProcess::Clear()
     296void G4VEnergyLossProcess::Clean()
    294297{
    295298  if(1 < verboseLevel) {
     
    302305  delete [] theCrossSectionMax;
    303306  delete [] idxSCoffRegions;
     307  delete [] idxDERegions;
    304308
    305309  theDEDXAtMaxEnergy = 0;
     
    312316  scProcesses.clear();
    313317  nProcesses = 0;
     318}
     319
     320//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     321
     322G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
     323                                                const G4Material*,
     324                                                G4double cut)
     325{
     326  return cut;
    314327}
    315328
     
    364377  }
    365378
    366   Clear();
     379  Clean();
    367380
    368381  // Base particle and set of models can be defined here
     
    407420                                     minSubRange, verboseLevel);
    408421
    409   // Sub Cutoff Regime
    410   if (nSCoffRegions>0) {
     422  // Sub Cutoff and Deexcitation
     423  if (nSCoffRegions>0 || nDERegions>0) {
    411424    theSubCuts = modelManager->SubCutoff();
    412425
     
    414427          G4ProductionCutsTable::GetProductionCutsTable();
    415428    size_t numOfCouples = theCoupleTable->GetTableSize();
    416     idxSCoffRegions = new G4int[numOfCouples];
     429
     430    if(nSCoffRegions>0) idxSCoffRegions = new G4bool[numOfCouples];
     431    if(nDERegions>0) idxDERegions = new G4bool[numOfCouples];
    417432 
    418433    for (size_t j=0; j<numOfCouples; j++) {
     
    421436        theCoupleTable->GetMaterialCutsCouple(j);
    422437      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;
     438     
     439      if(nSCoffRegions>0) {
     440        G4bool reg = false;
     441        for(G4int i=0; i<nSCoffRegions; i++) {
     442          if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = true;
     443        }
     444        idxSCoffRegions[j] = reg;
     445      }
     446      if(nDERegions>0) {
     447        G4bool reg = false;
     448        for(G4int i=0; i<nDERegions; i++) {
     449          if( pcuts == deRegions[i]->GetProductionCuts()) reg = true;
     450        }
     451        idxDERegions[j] = reg;
     452      }
    428453    }
    429454  }
     
    445470      }
    446471    }
     472    if (nDERegions) {
     473      G4cout << " Deexcitation is ON for regions: " << G4endl;
     474      for (G4int i=0; i<nDERegions; i++) {
     475        const G4Region* r = deRegions[i];
     476        G4cout << "           " << r->GetName() << G4endl;
     477      }
     478    }
    447479  }
    448480}
     
    479511    if(isIonisation) G4cout << "  isIonisation  flag = 1";
    480512    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;
    501513  }
    502514}
     
    640652//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    641653
    642 G4double G4VEnergyLossProcess::GetContinuousStepLimit(
    643                 const G4Track&,
    644                 G4double, G4double, G4double&)
    645 {
    646   return DBL_MAX;
     654void G4VEnergyLossProcess::PrintInfoDefinition()
     655{
     656  if(0 < verboseLevel) {
     657    G4cout << G4endl << GetProcessName() << ":   for  "
     658           << particle->GetParticleName()
     659           << "    SubType= " << GetProcessSubType()
     660           << G4endl
     661           << "      dE/dx and range tables from "
     662           << G4BestUnit(minKinEnergy,"Energy")
     663           << " to " << G4BestUnit(maxKinEnergy,"Energy")
     664           << " in " << nBins << " bins" << G4endl
     665           << "      Lambda tables from threshold to "
     666           << G4BestUnit(maxKinEnergy,"Energy")
     667           << " in " << nBins << " bins, spline: "
     668           << (G4LossTableManager::Instance())->SplineFlag()
     669           << G4endl;
     670    if(theRangeTableForLoss && isIonisation) {
     671      G4cout << "      finalRange(mm)= " << finalRange/mm
     672             << ", dRoverRange= " << dRoverRange
     673             << ", integral: " << integral
     674             << ", fluct: " << lossFluctuationFlag
     675             << ", linLossLimit= " << linLossLimit
     676             << G4endl;
     677    }
     678    PrintInfo();
     679    modelManager->DumpModelList(verboseLevel);
     680    if(theCSDARangeTable && isIonisation) {
     681      G4cout << "      CSDA range table up"
     682             << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
     683             << " in " << nBinsCSDA << " bins" << G4endl;
     684    }
     685    if(nSCoffRegions>0 && isIonisation) {
     686      G4cout << "      Subcutoff sampling in " << nSCoffRegions
     687             << " regions" << G4endl;
     688    }
     689    if(2 < verboseLevel) {
     690      G4cout << "      DEDXTable address= " << theDEDXTable << G4endl;
     691      if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
     692      G4cout << "non restricted DEDXTable address= "
     693             << theDEDXunRestrictedTable << G4endl;
     694      if(theDEDXunRestrictedTable && isIonisation) {
     695           G4cout << (*theDEDXunRestrictedTable) << G4endl;
     696      }
     697      if(theDEDXSubTable && isIonisation) {
     698        G4cout << (*theDEDXSubTable) << G4endl;
     699      }
     700      G4cout << "      CSDARangeTable address= " << theCSDARangeTable
     701             << G4endl;
     702      if(theCSDARangeTable && isIonisation) {
     703        G4cout << (*theCSDARangeTable) << G4endl;
     704      }
     705      G4cout << "      RangeTableForLoss address= " << theRangeTableForLoss
     706             << G4endl;
     707      if(theRangeTableForLoss && isIonisation) {
     708             G4cout << (*theRangeTableForLoss) << G4endl;
     709      }
     710      G4cout << "      InverseRangeTable address= " << theInverseRangeTable
     711             << G4endl;
     712      if(theInverseRangeTable && isIonisation) {
     713             G4cout << (*theInverseRangeTable) << G4endl;
     714      }
     715      G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
     716      if(theLambdaTable && isIonisation) {
     717        G4cout << (*theLambdaTable) << G4endl;
     718      }
     719      G4cout << "      SubLambdaTable address= " << theSubLambdaTable << G4endl;
     720      if(theSubLambdaTable && isIonisation) {
     721        G4cout << (*theSubLambdaTable) << G4endl;
     722      }
     723    }
     724  }
     725}
     726
     727//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     728
     729void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
     730{
     731  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     732  const G4Region* reg = r;
     733  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
     734
     735  // the region is in the list
     736  if (nSCoffRegions) {
     737    for (G4int i=0; i<nSCoffRegions; i++) {
     738      if (reg == scoffRegions[i]) {
     739        if(!val) deRegions[i] = 0;
     740        return;
     741      }
     742    }
     743  }
     744
     745  // new region
     746  if(val) {
     747    useSubCutoff = true;
     748    scoffRegions.push_back(reg);
     749    nSCoffRegions++;
     750  } else {
     751    useSubCutoff = false;
     752  }
     753}
     754
     755//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     756
     757void G4VEnergyLossProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
     758{
     759  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     760  const G4Region* reg = r;
     761  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
     762
     763  // the region is in the list
     764  if (nDERegions) {
     765    for (G4int i=0; i<nDERegions; i++) {
     766      if (reg == deRegions[i]) {
     767        if(!val) deRegions[i] = 0;
     768        return;
     769      }
     770    }
     771  }
     772
     773  // new region
     774  if(val) {
     775    useDeexcitation = true;
     776    deRegions.push_back(reg);
     777    nDERegions++;
     778  } else {
     779    useDeexcitation = false;
     780  }
    647781}
    648782
     
    674808  //  <<" stepLimit= "<<x<<G4endl;
    675809  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);
    688810}
    689811
     
    806928  if (length >= fRange) {
    807929    eloss = preStepKinEnergy;
    808     currentModel->CorrectionsAlongStep(currentCouple, dynParticle,
    809                                        eloss, esecdep, length);
     930    if (useDeexcitation) {
     931      if(idxDERegions[currentMaterialIndex]) {
     932        currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
     933        if(eloss < 0.0) eloss = 0.0;
     934      }
     935    }
    810936    fParticleChange.SetProposedKineticEnergy(0.0);
    811     fParticleChange.ProposeLocalEnergyDeposit(preStepKinEnergy);
     937    fParticleChange.ProposeLocalEnergyDeposit(eloss);
    812938    return &fParticleChange;
    813939  }
     
    9341060      G4double tmax =
    9351061        std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut);
     1062      G4double emean = eloss;
    9361063      eloss = fluc->SampleFluctuations(currentMaterial,dynParticle,
    937                                        tmax,length,eloss);
     1064                                       tmax,length,emean);
    9381065      /*                           
    9391066      if(-1 < verboseLevel)
     
    9501077  eloss += esecdep;
    9511078  if(eloss < 0.0) eloss = 0.0;
     1079
     1080  // deexcitation
     1081  else if (useDeexcitation) {
     1082    if(idxDERegions[currentMaterialIndex]) {
     1083      currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
     1084      if(eloss < 0.0) eloss = 0.0;
     1085    }
     1086  }
    9521087
    9531088  // Energy balanse
     
    11061241
    11071242  SelectModel(postStepScaledEnergy);
     1243  if(useDeexcitation) {
     1244    currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
     1245  }
     1246
    11081247  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
    11091248  G4double tcut = (*theCuts)[currentMaterialIndex];
     
    11401279//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    11411280
    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()
     1281G4bool G4VEnergyLossProcess::StorePhysicsTable(
     1282       const G4ParticleDefinition* part, const G4String& directory,
     1283       G4bool ascii)
     1284{
     1285  G4bool res = true;
     1286  if ( baseParticle || part != particle ) return res;
     1287
     1288  if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX"))
     1289    {res = false;}
     1290
     1291  if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr"))
     1292    {res = false;}
     1293
     1294  if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX"))
     1295    {res = false;}
     1296
     1297  if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation"))
     1298    {res = false;}
     1299
     1300  if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation"))
     1301    {res = false;}
     1302
     1303  if(isIonisation &&
     1304     !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange"))
     1305    {res = false;}
     1306
     1307  if(isIonisation &&
     1308     !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range"))
     1309    {res = false;}
     1310 
     1311  if(isIonisation &&
     1312     !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange"))
     1313    {res = false;}
     1314 
     1315  if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda"))
     1316    {res = false;}
     1317
     1318  if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda"))
     1319    {res = false;}
     1320
     1321  if ( res ) {
     1322    if(0 < verboseLevel) {
     1323      G4cout << "Physics tables are stored for " << particle->GetParticleName()
     1324             << " and process " << GetProcessName()
     1325             << " in the directory <" << directory
     1326             << "> " << G4endl;
     1327    }
     1328  } else {
     1329    G4cout << "Fail to store Physics Tables for "
     1330           << particle->GetParticleName()
     1331           << " and process " << GetProcessName()
     1332           << " in the directory <" << directory
     1333           << "> " << G4endl;
     1334  }
     1335  return res;
     1336}
     1337
     1338//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1339
     1340G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
     1341       const G4ParticleDefinition* part, const G4String& directory,
     1342       G4bool ascii)
     1343{
     1344  G4bool res = true;
     1345  const G4String particleName = part->GetParticleName();
     1346
     1347  if(1 < verboseLevel) {
     1348    G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
     1349           << particleName << " and process " << GetProcessName()
     1350           << "; tables_are_built= " << tablesAreBuilt
    11571351           << 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
     1352  }
     1353  if(particle == part) {
     1354
     1355    //    G4bool yes = true;
     1356    if ( !baseParticle ) {
     1357
     1358      G4bool fpi = true;
     1359      if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi))
     1360        {fpi = false;}
     1361
     1362      if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false))
     1363        {fpi = false;}
     1364
     1365      if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi))
     1366        {res = false;}
     1367
     1368      if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false))
     1369        {res = false;}
     1370
     1371      if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false))
     1372        {res = false;}
     1373
     1374      if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi))
     1375        {res = false;}
     1376
     1377      if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true))
     1378        {res = false;}
     1379
     1380      G4bool yes = false;
     1381      if(nSCoffRegions > 0) {yes = true;}
     1382
     1383      if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes))
     1384        {res = false;}
     1385
     1386      if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes))
     1387        {res = false;}
     1388
     1389      if(!fpi) yes = false;
     1390      if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes))
     1391        {res = false;}
     1392    }
     1393  }
     1394
     1395  return res;
     1396}
     1397
     1398//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1399
     1400G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part,
     1401                                        G4PhysicsTable* aTable, G4bool ascii,
     1402                                        const G4String& directory,
     1403                                        const G4String& tname)
     1404{
     1405  G4bool res = true;
     1406  if ( aTable ) {
     1407    const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
     1408    if( !aTable->StorePhysicsTable(name,ascii)) res = false;
     1409  }
     1410  return res;
     1411}
     1412
     1413//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1414
     1415G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part,
     1416                                           G4PhysicsTable* aTable, G4bool ascii,
     1417                                           const G4String& directory,
     1418                                           const G4String& tname,
     1419                                           G4bool mandatory)
     1420{
     1421  G4bool res = true;
     1422  G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
     1423  G4bool yes = aTable->ExistPhysicsTable(filename);
     1424  if(yes) {
     1425    yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii);
     1426    if((G4LossTableManager::Instance())->SplineFlag()) {
     1427      size_t n = aTable->length();
     1428      for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);}
     1429    }
     1430  }
     1431  if(yes) {
     1432    if (0 < verboseLevel) {
     1433      G4cout << tname << " table for " << part->GetParticleName()
     1434             << " is Retrieved from <" << filename << ">"
    11891435             << G4endl;
    1190       if(theCSDARangeTable && isIonisation) {
    1191         G4cout << (*theCSDARangeTable) << G4endl;
    1192       }
    1193       G4cout << "      RangeTableForLoss address= " << theRangeTableForLoss
     1436    }
     1437  } else {
     1438    if(mandatory) res = false;
     1439    if(mandatory || 1 < verboseLevel) {
     1440      G4cout << tname << " table for " << part->GetParticleName()
     1441             << " from file <"
     1442             << filename << "> is not Retrieved"
    11941443             << 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       }
     1444    }
     1445  }
     1446  return res;
     1447}
     1448
     1449//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1450
     1451G4double G4VEnergyLossProcess::GetDEDXDispersion(
     1452                                  const G4MaterialCutsCouple *couple,
     1453                                  const G4DynamicParticle* dp,
     1454                                        G4double length)
     1455{
     1456  DefineMaterial(couple);
     1457  G4double ekin = dp->GetKineticEnergy();
     1458  SelectModel(ekin*massRatio);
     1459  G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
     1460  tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]);
     1461  G4double d = 0.0;
     1462  G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
     1463  if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length);
     1464  return d;
     1465}
     1466
     1467//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1468
     1469G4double G4VEnergyLossProcess::CrossSectionPerVolume(
     1470         G4double kineticEnergy, const G4MaterialCutsCouple* couple)
     1471{
     1472  // Cross section per volume is calculated
     1473  DefineMaterial(couple);
     1474  G4double cross = 0.0;
     1475  G4bool b;
     1476  if(theLambdaTable) {
     1477    cross =
     1478      ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b);
     1479  } else {
     1480    SelectModel(kineticEnergy);
     1481    cross =
     1482      currentModel->CrossSectionPerVolume(currentMaterial,
     1483                                          particle, kineticEnergy,
     1484                                          (*theCuts)[currentMaterialIndex]);
     1485  }
     1486
     1487  return cross;
     1488}
     1489
     1490//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1491
     1492G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track)
     1493{
     1494  DefineMaterial(track.GetMaterialCutsCouple());
     1495  preStepLambda = GetLambdaForScaledEnergy(track.GetKineticEnergy()*massRatio);
     1496  G4double x = DBL_MAX;
     1497  if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
     1498  return x;
     1499}
     1500
     1501//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1502
     1503G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track,
     1504                                                   G4double x, G4double y,
     1505                                                   G4double& z)
     1506{
     1507  G4GPILSelection sel;
     1508  return AlongStepGetPhysicalInteractionLength(track, x, y, z, &sel);
     1509}
     1510
     1511//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1512
     1513G4double G4VEnergyLossProcess::GetMeanFreePath(
     1514                             const G4Track& track,
     1515                             G4double,
     1516                             G4ForceCondition* condition)
     1517
     1518{
     1519  *condition = NotForced;
     1520  return MeanFreePath(track);
     1521}
     1522
     1523//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1524
     1525G4double G4VEnergyLossProcess::GetContinuousStepLimit(
     1526                const G4Track&,
     1527                G4double, G4double, G4double&)
     1528{
     1529  return DBL_MAX;
     1530}
     1531
     1532//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1533
     1534G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
     1535                 const G4MaterialCutsCouple* couple, G4double cut)
     1536{
     1537  //  G4double cut  = (*theCuts)[couple->GetIndex()];
     1538  //  G4int nbins = nLambdaBins;
     1539  G4double tmin =
     1540    std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
     1541             minKinEnergy);
     1542  if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
     1543  G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
     1544  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
     1545
     1546  return v;
     1547}
     1548
     1549//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1550 
     1551void G4VEnergyLossProcess::AddCollaborativeProcess(
     1552            G4VEnergyLossProcess* p)
     1553{
     1554  G4bool add = true;
     1555  if(p->GetProcessName() != "eBrem") add = false;
     1556  if(add && nProcesses > 0) {
     1557    for(G4int i=0; i<nProcesses; i++) {
     1558      if(p == scProcesses[i]) {
     1559        add = false;
     1560        break;
     1561      }
     1562    }
     1563  }
     1564  if(add) {
     1565    scProcesses.push_back(p);
     1566    nProcesses++;
     1567    if (1 < verboseLevel) {
     1568      G4cout << "### The process " << p->GetProcessName()
     1569             << " is added to the list of collaborative processes of "
     1570             << GetProcessName() << G4endl;
    12111571    }
    12121572  }
     
    13771737//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    13781738
    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

    r991 r1005  
    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.10 2009/02/24 09:56:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    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),
     59  safetyHelper(0),
    5760  facrange(0.02),
    5861  facgeom(2.5),
     
    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
     79void G4VMscModel::InitialiseSafetyHelper()
     80{
     81  if(!safetyHelper) {
     82    safetyHelper = G4TransportationManager::GetTransportationManager()
     83      ->GetSafetyHelper();
     84    safetyHelper->InitialiseHelper();
     85  }
     86}
     87
     88//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     89
     90void G4VMscModel::ComputeDisplacement(G4ParticleChangeForMSC* fParticleChange, 
     91                                      const G4ThreeVector& dir,
     92                                      G4double displacement,
     93                                      G4double postsafety)
     94{
     95  const G4ThreeVector* pos = fParticleChange->GetProposedPosition();
     96  G4double r = displacement;
     97  if(r >  postsafety) {
     98    G4double newsafety = safetyHelper->ComputeSafety(*pos);
     99    if(r > newsafety) r = newsafety;
     100  }
     101  if(r > 0.) {
     102
     103    // compute new endpoint of the Step
     104    G4ThreeVector newPosition = *pos + r*dir;
     105
     106    // definitely not on boundary
     107    if(displacement == r) {
     108      safetyHelper->ReLocateWithinVolume(newPosition);
     109
     110    } else {
     111      // check safety after displacement
     112      G4double postsafety = safetyHelper->ComputeSafety(newPosition);
     113
     114      // displacement to boundary
     115      if(postsafety <= 0.0) {
     116        safetyHelper->Locate(newPosition,
     117                             *fParticleChange->GetProposedMomentumDirection());
     118
     119        // not on the boundary
     120      } else {
     121        safetyHelper->ReLocateWithinVolume(newPosition);
     122      }
     123    }
     124    fParticleChange->ProposePosition(newPosition);
     125  }
     126}
     127
     128//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     129
     130void G4VMscModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
     131                                    const G4MaterialCutsCouple*,
     132                                    const G4DynamicParticle*,
     133                                    G4double, G4double)
     134{}
     135
     136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracChangeset for help on using the changeset viewer.