Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (14 years ago)
Author:
garnier
Message:

update ti head

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

Legend:

Unmodified
Added
Removed

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

    r1315 r1340  
    1 $Id: History,v 1.424 2010/06/04 15:36:39 vnivanch Exp $
     1$Id: History,v 1.436 2010/11/04 12:55:09 vnivanch Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     2004 November 10: V.Ivant (emutils-V09-03-23)
     21- Fixed part of problems reported by the Coverity tool
     22  (mainly initialisation)
     23
     2425 October 10: V.Ivant (emutils-V09-03-22)
     25- Fixed problems reported by the Coverity tool (mainly initialisation)
     26
     2715 October 10: V.Ivant (emutils-V09-03-21)
     28- G4VEmAngularDistribution - new general interface
     29- G4VBremAngularDistribution moved from lowenergy
     30- G4VEmModel - G4VEmAngularDistribution added and Get/Set methods
     31- G4VEnergyLossProcess - fixed bug #1141 (L.Pandola)
     32
     3307 September 10: V.Ivant (emutils-V09-03-20)
     34- G4VMscModel - in computation of the displacement added
     35                a protection limited displacement to be incide current
     36                volume (address bug #1128)
     37- G4LossTableManager - a little cleanup of interfaces
     38- G4ElectronIonPair - added method SampleNumberOfIonsAlongStep
     39- G4EmCalculator - return back tag 18
     40
     4123 August 10: V.Ivant (emutils-V09-03-19)
     42G4VEnergyLossProcess - minor optimisation of PostStepDoIt (in some cases
     43  one call to the log of random number less)
     44
     4517 August 10: V.Ivant (emutils-V09-03-18)
     46G4VEmProcess, G4VEnergyLossProcess, G4VMultipleScattering, G4VEmModel,
     47G4EmMultiModel - substituted obsolete method GetDefinition() of the class
     48                 G4DynamicParticle by the new one GetParticleDefinition()
     49
     5029 July 10: V.Ivant (emutils-V09-03-17)
     51G4EmMultiModel - fixed and cleaned up
     52G4VMultipleScattering - added more detailed printout for kaon+
     53G4LossTableManager - added pointer and access method to G4ElectronIonPair
    1954
    20554 June 10: V.Ivant (emutils-V09-03-16)

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4ElectronIonPair.hh,v 1.2 2008/10/17 14:46:16 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4ElectronIonPair.hh,v 1.5 2010/10/25 17:23:01 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929//
    … …  
    6262#include "G4ParticleDefinition.hh"
    6363#include "G4ThreeVector.hh"
    64 #include "G4TrackVector.hh"
    6564#include "G4VProcess.hh"
    6665#include <vector>
    … …  
    8685  inline G4double MeanNumberOfIonsAlongStep(const G4Step*);
    8786
     87  inline G4int SampleNumberOfIonsAlongStep(const G4Step*);
     88
    8889  // returns pointer to the new vector of positions of
    8990  // ionisation points in the World coordinate system
    90   std::vector<G4ThreeVector>*
    91   SampleIonsAlongStep(const G4ThreeVector& prePosition,
    92                       const G4ThreeVector& postPosition,
    93                       G4double numberOfIonisations);
    94 
    9591  std::vector<G4ThreeVector>* SampleIonsAlongStep(const G4Step*);
    9692
    … …  
    115111private:
    116112
     113  void Initialise();
     114
     115  G4double FindMeanEnergyPerIonPair(const G4Material*);
     116
    117117  // hide assignment operator
    118118  G4ElectronIonPair & operator=(const G4ElectronIonPair &right);
    119119  G4ElectronIonPair(const G4ElectronIonPair&);
    120120
    121   G4double FindMeanEnergyPerIonPair(const G4Material*);
     121  // cash
     122  const G4Material*  curMaterial;
     123  G4double           curMeanEnergy;
    122124
    123   void Initialise();
    124 
    125   const G4ParticleDefinition* gamma;
    126 
    127   // cash
    128   const G4Material*           curMaterial;
    129   G4double                    curMeanEnergy;
    130 
     125  G4double FanoFactor;
     126 
    131127  G4int    verbose;             
    132128  G4int    nMaterials;
    133129
    134130  // list of G4 NIST materials with mean energy per ion defined
    135   std::vector<G4double>       g4MatData;
    136   std::vector<G4String>       g4MatNames;
    137 
     131  std::vector<G4double>  g4MatData;
     132  std::vector<G4String>  g4MatNames;
    138133};
    139134
    … …  
    141136G4ElectronIonPair::MeanNumberOfIonsAlongStep(const G4Step* step)
    142137{
    143   return MeanNumberOfIonsAlongStep(step->GetTrack()->GetDefinition(),
     138  return MeanNumberOfIonsAlongStep(step->GetTrack()->GetParticleDefinition(),
    144139                                   step->GetPreStepPoint()->GetMaterial(),
    145140                                   step->GetTotalEnergyDeposit(),
    … …  
    147142}
    148143
    149 inline std::vector<G4ThreeVector>*
    150 G4ElectronIonPair::SampleIonsAlongStep(const G4Step* step)
     144inline
     145G4int G4ElectronIonPair::SampleNumberOfIonsAlongStep(const G4Step* step)
    151146{
    152   return SampleIonsAlongStep(step->GetPreStepPoint()->GetPosition(),
    153                              step->GetPostStepPoint()->GetPosition(),
    154                              MeanNumberOfIonsAlongStep(step));
    155 }
     147  G4double meanion = MeanNumberOfIonsAlongStep(step);
     148  G4double sig = FanoFactor*std::sqrt(meanion);
     149  G4int nion = G4int(G4RandGauss::shoot(meanion,sig) + 0.5);
     150  return nion;
     151}
    156152
    157153inline
    … …  
    160156  G4int subtype = -1;
    161157  const G4VProcess* proc = step->GetPostStepPoint()->GetProcessDefinedStep();
    162   if(proc) subtype = proc->GetProcessSubType();
    163   return ResidualeChargePostStep(step->GetTrack()->GetDefinition(),
     158  if(proc) { subtype = proc->GetProcessSubType(); }
     159  return ResidualeChargePostStep(step->GetTrack()->GetParticleDefinition(),
    164160                                 step->GetSecondary(),
    165161                                 subtype);

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmMultiModel.hh,v 1.6 2007/05/22 17:31:57 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4EmMultiModel.hh,v 1.7 2010/07/04 17:51:09 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4040// Modifications:
    4141// 15-04-05 optimize internal interface (V.Ivanchenko)
     42// 04-07-10 updated interfaces according to g4 9.4 (V.Ivanchenko)
    4243//
    4344// Class Description:
    4445//
    45 // Energy loss model using several G4VEmModels
     46// EM model using several G4VEmModels for the same energy interval
    4647
    4748// -------------------------------------------------------------------
    … …  
    5556#include <vector>
    5657
    57 class G4Region;
    58 class G4PhysicsTable;
    5958class G4DynamicParticle;
    6059
    … …  
    6867  virtual ~G4EmMultiModel();
    6968
    70   void Initialise(const G4ParticleDefinition*, const G4DataVector&);
     69  void AddModel(G4VEmModel*);
    7170
    72   G4double MinEnergyCut(const G4ParticleDefinition*,
    73                         const G4MaterialCutsCouple*);
     71  virtual void Initialise(const G4ParticleDefinition*,
     72                          const G4DataVector&);
    7473
     74  virtual G4double ComputeDEDX(const G4MaterialCutsCouple*,
     75                               const G4ParticleDefinition*,
     76                               G4double kineticEnergy,
     77                               G4double cutEnergy);
    7578
    76   G4double ComputeDEDX(const G4MaterialCutsCouple*,
    77                        const G4ParticleDefinition*,
    78                        G4double kineticEnergy,
    79                        G4double cutEnergy);
     79  // main method to compute cross section per atom
     80  virtual
     81  G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
     82                                      G4double kinEnergy,
     83                                      G4double Z,
     84                                      G4double A = 0., /* amu */
     85                                      G4double cutEnergy = 0.0,
     86                                      G4double maxEnergy = DBL_MAX);
    8087
    81   G4double CrossSection(const G4MaterialCutsCouple*,
    82                         const G4ParticleDefinition*,
    83                         G4double kineticEnergy,
    84                         G4double cutEnergy,
    85                         G4double maxEnergy);
    86 
    87   void SampleSecondaries(std::vector<G4DynamicParticle*>*,
    88                          const G4MaterialCutsCouple*,
    89                          const G4DynamicParticle*,
    90                          G4double tmin,
    91                          G4double tmax);
    92 
    93   void DefineForRegion(const G4Region*);
    94 
    95   void AddModel(G4VEmModel*, G4double tmin, G4double tmax);
    96 
    97 protected:
    98 
    99   G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
    100                                     G4double kineticEnergy);
     88  virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
     89                                 const G4MaterialCutsCouple*,
     90                                 const G4DynamicParticle*,
     91                                 G4double tmin,
     92                                 G4double tmax);
    10193
    10294private:
    … …  
    108100  G4int                         nModels;
    109101  std::vector<G4VEmModel*>      model;
    110   G4DataVector                  tsecmin;
    111   G4DataVector                  cross_section;
     102  std::vector<G4double>         cross_section;
    112103
    113104};

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmSaturation.hh,v 1.7 2009/09/25 09:16:40 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4EmSaturation.hh,v 1.8 2010/08/17 17:36:58 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929//
    … …  
    140140{
    141141  G4Track* track = step->GetTrack();
    142   return VisibleEnergyDeposition(track->GetDefinition(),
     142  return VisibleEnergyDeposition(track->GetParticleDefinition(),
    143143                                 track->GetMaterialCutsCouple(),
    144144                                 step->GetStepLength(),

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LossTableManager.hh,v 1.58 2010/04/27 16:59:52 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4LossTableManager.hh,v 1.61 2010/09/03 10:09:45 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929//
    … …  
    9393class G4EmSaturation;
    9494class G4EmConfigurator;
     95class G4ElectronIonPair;
    9596class G4LossTableBuilder;
     97class G4VAtomDeexcitation;
    9698class G4Region;
    9799
    … …  
    259261  const std::vector<G4VMultipleScattering*>& GetMultipleScatteringVector();
    260262
    261   inline G4VEnergyLossProcess* GetEnergyLossProcess(const G4ParticleDefinition*);
     263  inline
     264  G4VEnergyLossProcess* GetEnergyLossProcess(const G4ParticleDefinition*);
    262265
    263266  G4EmCorrections* EmCorrections();
    … …  
    266269
    267270  G4EmConfigurator* EmConfigurator();
     271
     272  G4ElectronIonPair* ElectronIonPair();
     273
     274  G4VAtomDeexcitation* AtomDeexcitation();
     275
     276  void SetAtomDeexcitation(G4VAtomDeexcitation*);
    268277
    269278private:
    … …  
    345354  G4EmSaturation*             emSaturation;
    346355  G4EmConfigurator*           emConfigurator;
     356  G4ElectronIonPair*          emElectronIonPair;
     357  G4VAtomDeexcitation*        atomDeexcitation;
    347358
    348359  G4int nbinsLambda;
    … …  
    464475          G4double& length)
    465476{
    466   const G4ParticleDefinition* aParticle = dp->GetDefinition();
     477  const G4ParticleDefinition* aParticle = dp->GetParticleDefinition();
    467478  if(aParticle != currentParticle) {
    468479    std::map<PD,G4VEnergyLossProcess*,std::less<PD> >::const_iterator pos;

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VAtomDeexcitation.hh,v 1.3 2010/03/30 09:19:56 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VAtomDeexcitation.hh,v 1.5 2010/10/14 16:27:35 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5151
    5252#include "globals.hh"
     53#include "G4AtomicShell.hh"
     54#include "G4ProductionCutsTable.hh"
     55#include "G4Track.hh"
    5356#include <vector>
    5457
    55 class G4AtomicShell;
    5658class G4ParticleDefinition;
    5759class G4DynamicParticle;
     60class G4MaterialCutsCouple;
     61class G4VParticleChange;
     62
     63enum G4AtomicShellEnumerator
     64{
     65  fKShell = 0,
     66  fL1Shell,
     67  fL2Shell,
     68  fL3Shell,
     69  fM1Shell,
     70  fM2Shell,
     71  fM3Shell,
     72  fM4Shell,
     73  fM5Shell
     74};
    5875
    5976class G4VAtomDeexcitation {
    60 
    61   G4VAtomDeexcitation(const G4String& pname = "");
     77public:
     78
     79  G4VAtomDeexcitation(const G4String& modname = "Deexcitation",
     80                      const G4String& pixename = "");
    6281
    6382  virtual ~G4VAtomDeexcitation();
    … …  
    6584  //========== initialization ==========
    6685
    67   virtual void PreparePhysicsTable(const G4ParticleDefinition&);
    68   virtual void BuildPhysicsTable(const G4ParticleDefinition&);
     86  // Overall initialisation before new run
     87  void InitialiseAtomicDeexcitation();
     88
     89  // Initialisation of deexcitation at the beginning of run
     90  virtual void InitialiseForNewRun() = 0;
     91
     92  // Initialisation for a concrete atom
     93  // May be called at run time
     94  virtual void InitialiseForExtraAtom(G4int Z) = 0;
     95
     96  // Activation of deexcitation per detector region
     97  void SetDeexcitationActiveRegion(const G4String& rname = "");
     98
     99  // Activation of Auger electron production
     100  inline void SetAugerActive(G4bool);
     101  inline G4bool IsAugerActive() const;
     102
     103  // Activation of PIXE simulation
     104  inline void SetPIXEActive(G4bool);
     105  inline G4bool IsPIXEActive() const;
     106
     107  // Deexcitation model name
     108  inline const G4String& GetName() const;
    69109
    70110  // PIXE model name
    … …  
    72112  inline const G4String& PIXECrossSectionModel() const;
    73113
    74   // Activation of deexcitation per detector region
    75   void SetFluorescenceActiveRegion(const G4Region* region = 0);
    76   void SetAugerActiveRegion(const G4Region* region = 0);
    77   void SetPIXECrossSectionActiveRegion(const G4Region* region = 0);
    78 
    79   void SetFluorescenceActiveRegion(const G4String& rname = "");
    80   void SetAugerActiveRegion(const G4String& rname = "");
    81   void SetPIXECrossSectionActiveRegion(const G4String& rname = "");
     114  // Access to the list of atoms active for deexcitation
     115  inline const std::vector<G4bool>& GetListOfActiveAtoms() const;
     116
     117  // Verbosity level
     118  inline void SetVerboseLevel(G4int);
     119  inline G4int GetVerboseLevel() const;
    82120
    83121  //========== Run time methods ==========
    84122
    85123  // Check if deexcitation is active for a given geometry volume
    86   G4bool CheckFluorescenceActiveRegion(G4int coupleIndex);
    87 
    88   // Check if deexcitation is active for a given geometry volume
    89   G4bool CheckPIXEActiveRegion(G4int coupleIndex);
     124  inline G4bool CheckDeexcitationActiveRegion(G4int coupleIndex);
    90125
    91126  // Get atomic shell by shell index, used by discrete processes
    92127  // (for example, photoelectric), when shell vacancy sampled by the model
    93   const G4AtomicShell* GetAtomicShell(G4int Z, G4int ShellIndex);
    94 
    95   // selection of random shell for ionisation process
    96   virtual const G4AtomicShell* SelectRandomShell(const G4DynamicParticle*,
    97                                                  G4int Z);
     128  virtual
     129  const G4AtomicShell* GetAtomicShell(G4int Z,
     130                                      G4AtomicShellEnumerator shell) = 0;
    98131
    99132  // generation of deexcitation for given atom and shell vacancy
    100   virtual void GenerateParticles(std::vector<G4DynamicParticle*>*,
    101                                  const G4AtomicShell*, G4int Z) = 0;
     133  inline void GenerateParticles(std::vector<G4DynamicParticle*>* secVect, 
     134                                const G4AtomicShell*,
     135                                G4int Z,
     136                                G4int coupleIndex);
     137
     138  // generation of deexcitation for given atom and shell vacancy
     139  virtual void GenerateParticles(std::vector<G4DynamicParticle*>* secVect, 
     140                                 const G4AtomicShell*,
     141                                 G4int Z,
     142                                 G4double gammaCut,
     143                                 G4double eCut) = 0;
    102144
    103145  // access or compute PIXE cross section
    104   virtual G4double GetPIXECrossSection (const G4ParticleDefinition*,
    105                                         G4int Z, G4double kinE) = 0;
    106 
    107   // calculate PIXE cross section from the models
    108   virtual G4double CalculatePIXECrossSection(const G4ParticleDefinition*,
    109                                              G4int Z, G4double kinE) = 0;
     146  virtual G4double
     147  GetShellIonisationCrossSectionPerAtom(const G4ParticleDefinition*,
     148                                        G4int Z,
     149                                        G4AtomicShellEnumerator shell,
     150                                        G4double kinE) = 0;
     151
     152  // access or compute PIXE cross section
     153  virtual G4double
     154  ComputeShellIonisationCrossSectionPerAtom(const G4ParticleDefinition*,
     155                                            G4int Z,
     156                                            G4AtomicShellEnumerator shell,
     157                                            G4double kinE) = 0;
    110158
    111159  // Sampling of PIXE for ionisation processes
    112   virtual void
    113   AlongStepDeexcitation(std::vector<G4DynamicParticle*>* secVect,
    114                         const G4DynamicParticle* icidentParticle,
    115                         const G4MaterialCutsCouple*,
    116                         G4double trueStepLenght,
    117                         G4double eLoss) = 0;
    118 
    119 
     160  void AlongStepDeexcitation(G4VParticleChange* pParticleChange, 
     161                             const G4Step& step,
     162                             G4double& eLoss,
     163                             G4int coupleIndex);
    120164
    121165private:
    … …  
    125169  G4VAtomDeexcitation & operator=(const G4VAtomDeexcitation &right);
    126170
     171  G4ProductionCutsTable* theCoupleTable;
     172  G4int    verbose;
     173  G4String name;
    127174  G4String namePIXE;
    128 
     175  G4bool   flagAuger;
     176  G4bool   flagPIXE;
     177  std::vector<G4bool>   activeZ;
     178  std::vector<G4bool>   activeDeexcitationMedia;
     179  std::vector<G4String> activeRegions;
     180  std::vector<G4DynamicParticle*> vdyn;
     181  std::vector<G4Track*> secVect;
    129182};
     183
     184inline void G4VAtomDeexcitation::SetAugerActive(G4bool val)
     185{
     186  flagAuger = val;
     187}
     188
     189inline G4bool G4VAtomDeexcitation::IsAugerActive() const
     190{
     191  return flagAuger;
     192}
     193
     194inline void G4VAtomDeexcitation::SetPIXEActive(G4bool val)
     195{
     196  flagPIXE = val;
     197}
     198
     199inline G4bool G4VAtomDeexcitation::IsPIXEActive() const
     200{
     201  return flagPIXE;
     202}
     203
     204inline const G4String& G4VAtomDeexcitation::GetName() const
     205{
     206  return name;
     207}
    130208
    131209inline
    … …  
    141219}
    142220
     221inline const std::vector<G4bool>&
     222G4VAtomDeexcitation::GetListOfActiveAtoms() const
     223{
     224  return activeZ;
     225}
     226
     227inline void G4VAtomDeexcitation::SetVerboseLevel(G4int val)
     228{
     229  verbose = val;
     230}
     231
     232inline G4int G4VAtomDeexcitation::GetVerboseLevel() const
     233{
     234  return verbose;
     235}
     236
     237inline G4bool
     238G4VAtomDeexcitation::CheckDeexcitationActiveRegion(G4int coupleIndex)
     239{
     240  return activeDeexcitationMedia[coupleIndex];
     241}
     242
     243inline void
     244G4VAtomDeexcitation::GenerateParticles(std::vector<G4DynamicParticle*>* v, 
     245                                       const G4AtomicShell* as,
     246                                       G4int Z,
     247                                       G4int idx)
     248{
     249  G4double gCut = (*theCoupleTable->GetEnergyCutsVector(idx))[0];
     250  if(gCut < as->BindingEnergy()) {
     251    GenerateParticles(v, as, Z, gCut,
     252                      (*theCoupleTable->GetEnergyCutsVector(idx))[1]);
     253  }
     254}
     255
    143256#endif
    144257

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmModel.hh,v 1.75 2010/05/26 10:41:34 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VEmModel.hh,v 1.77 2010/10/14 16:27:35 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6767// 16-02-09 Moved implementations of virtual methods to source (VI)
    6868// 07-04-09 Moved msc methods from G4VEmModel to G4VMscModel (VI)
     69// 13-10-10 Added G4VEmAngularDistribution (VI)
    6970//
    7071// Class Description:
    … …  
    8788#include "G4DataVector.hh"
    8889#include "G4VEmFluctuationModel.hh"
     90#include "G4VEmAngularDistribution.hh"
    8991#include "G4EmElementSelector.hh"
    9092#include "Randomize.hh"
    … …  
    146148                                                                     
    147149  // min cut in kinetic energy allowed by the model
     150  // obsolete method will be removed
    148151  virtual G4double MinEnergyCut(const G4ParticleDefinition*,
    149152                                const G4MaterialCutsCouple*);
    … …  
    252255  //------------------------------------------------------------------------
    253256
     257  void SetParticleChange(G4VParticleChange*, G4VEmFluctuationModel* f=0);
     258
    254259  inline G4VEmFluctuationModel* GetModelOfFluctuations();
    255260
     261  inline G4VEmAngularDistribution* GetAngularDistribution();
     262
     263  inline void SetAngularDistribution(G4VEmAngularDistribution*);
     264
    256265  inline G4double HighEnergyLimit() const;
    257266
    … …  
    284293  inline void SetDeexcitationFlag(G4bool val);
    285294
    286   inline void ActivateNuclearStopping(G4bool);
    287 
    288295  inline G4double MaxSecondaryKinEnergy(const G4DynamicParticle* dynParticle);
    289296
    290297  inline const G4String& GetName() const;
    291 
    292   inline void SetParticleChange(G4VParticleChange*, G4VEmFluctuationModel*);
    293298
    294299  inline void SetCurrentCouple(const G4MaterialCutsCouple*);
    … …  
    310315  // ======== Parameters of the class fixed at construction =========
    311316
    312   G4VEmFluctuationModel* fluc;
     317  G4VEmFluctuationModel* flucModel;
     318  G4VEmAngularDistribution* anglModel;
    313319  const G4String   name;
    314320
    … …  
    329335
    330336  G4VParticleChange*  pParticleChange;
    331   G4bool              nuclearStopping;
     337  //  G4bool              nuclearStopping;
    332338
    333339  // ======== Cashed values - may be state dependent ================
    … …  
    377383G4double G4VEmModel::MaxSecondaryKinEnergy(const G4DynamicParticle* dynPart)
    378384{
    379   return MaxSecondaryEnergy(dynPart->GetDefinition(),
     385  return MaxSecondaryEnergy(dynPart->GetParticleDefinition(),
    380386                            dynPart->GetKineticEnergy());
    381387}
    … …  
    479485inline G4VEmFluctuationModel* G4VEmModel::GetModelOfFluctuations()
    480486{
    481   return fluc;
     487  return flucModel;
     488}
     489
     490//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     491
     492inline G4VEmAngularDistribution* G4VEmModel::GetAngularDistribution()
     493{
     494  return anglModel;
     495}
     496
     497//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     498
     499inline void G4VEmModel::SetAngularDistribution(G4VEmAngularDistribution* p)
     500{
     501  anglModel = p;
    482502}
    483503
    … …  
    587607}
    588608
    589 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    590 
    591 inline void G4VEmModel::ActivateNuclearStopping(G4bool val)
    592 {
    593   nuclearStopping = val;
    594 }
    595 
    596609//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    597610
    … …  
    603616//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    604617
    605 inline void G4VEmModel::SetParticleChange(G4VParticleChange* p, 
    606                                           G4VEmFluctuationModel* f = 0)
    607 {
    608   if(p && pParticleChange != p) { pParticleChange = p; }
    609   fluc = f;
    610 }
    611 
    612 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    613 
    614618#endif
    615619

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmProcess.hh,v 1.60 2010/04/28 14:43:13 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VEmProcess.hh,v 1.61 2010/08/17 17:36:59 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    413413inline void G4VEmProcess::InitialiseStep(const G4Track& track)
    414414{
    415   currentParticle = track.GetDefinition();
     415  currentParticle = track.GetParticleDefinition();
    416416  preStepKinEnergy = track.GetKineticEnergy();
    417417  DefineMaterial(track.GetMaterialCutsCouple());

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

    r1315 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEnergyLossProcess.hh,v 1.92 2010/04/28 14:43:13 vnivanch Exp $
     26// $Id: G4VEnergyLossProcess.hh,v 1.93 2010/10/14 16:27:35 vnivanch Exp $
    2727// GEANT4 tag $Name:
    2828//
    … …  
    114114class G4Region;
    115115class G4SafetyHelper;
     116class G4VAtomDeexcitation;
    116117
    117118//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    466467  std::vector<G4VEmModel*>              emModels;
    467468  G4VEmFluctuationModel*                fluctModel;
     469  G4VAtomDeexcitation*                  atomDeexcitation;
    468470  std::vector<const G4Region*>          scoffRegions;
    469471  std::vector<const G4Region*>          deRegions;

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMscModel.hh,v 1.9 2009/04/07 18:39:47 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VMscModel.hh,v 1.10 2010/09/07 16:05:33 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    140140  G4double dtrl;
    141141  G4double lambdalimit;
    142   G4double geommax;
     142  G4double geomMin;
     143  G4double geomMax;
    143144
    144145  G4MscStepLimitType steppingAlgorithm;
    … …  
    206207                                              G4double limit)
    207208{
    208   G4double res = geommax;
     209  G4double res = geomMax;
    209210  if(track.GetVolume() != safetyHelper->GetWorldVolume()) {
    210211    res = safetyHelper->CheckNextStep(

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4ElectronIonPair.cc,v 1.2 2008/10/17 14:46:16 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4ElectronIonPair.cc,v 1.5 2010/10/25 17:23:01 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4646
    4747#include "G4ElectronIonPair.hh"
    48 #include "G4Gamma.hh"
    4948#include "G4Material.hh"
    5049#include "G4MaterialTable.hh"
    … …  
    6362  curMeanEnergy = 0.0;
    6463  nMaterials = 0;
     64  FanoFactor = 0.2;
    6565  Initialise();
    6666}
    … …  
    9797        }
    9898      }
    99       if(curMeanEnergy > 0.0) nion = (edep - niel)/curMeanEnergy;
     99      if(curMeanEnergy > 0.0) { nion = (edep - niel)/curMeanEnergy; }
    100100    }
    101101  }
    … …  
    106106
    107107std::vector<G4ThreeVector>*
    108 G4ElectronIonPair::SampleIonsAlongStep(const G4ThreeVector& prePos,
    109                                        const G4ThreeVector& postPos,
    110                                        G4double meanion)
    111 {
    112   std::vector<G4ThreeVector>* v = new std::vector<G4ThreeVector>;
    113 
    114   G4double sig = 0.2*std::sqrt(meanion);
    115   G4int nion = G4int(G4RandGauss::shoot(meanion,sig) + 0.5);
     108G4ElectronIonPair::SampleIonsAlongStep(const G4Step* step)
     109{
     110  std::vector<G4ThreeVector>* v = 0;
     111
     112  G4int nion = SampleNumberOfIonsAlongStep(step);
    116113
    117114  // sample ionisation along step
    118115  if(nion > 0) {
    119116
    120     G4ThreeVector deltaPos = postPos - prePos; 
    121     for(G4int i=0; i<nion; i++) {
     117    v = new std::vector<G4ThreeVector>;
     118    G4ThreeVector prePos = step->GetPreStepPoint()->GetPosition();
     119    G4ThreeVector deltaPos = step->GetPostStepPoint()->GetPosition() - prePos; 
     120    for(G4int i=0; i<nion; ++i) {
    122121      v->push_back( prePos + deltaPos*G4UniformRand() );
    123122    }
    … …  
    138137  G4int nholes = 0;
    139138
    140   if(2 == subType || 12 == subType || 13 == subType) nholes = 1;
     139  if(2 == subType || 12 == subType || 13 == subType) { nholes = 1; }
    141140  return nholes;
    142141}
    … …  
    173172  if(nmat > 0) {
    174173    G4cout << "### G4ElectronIonPair: mean energy per ion pair avalable:" << G4endl;
    175     for(G4int i=0; i<nmat; i++) {
     174    for(G4int i=0; i<nmat; ++i) {
    176175      const G4Material* mat = (*mtable)[i];
    177176      G4double x = mat->GetIonisation()->GetMeanEnergyPerIonPair();
    … …  
    191190    G4cout << "### G4ElectronIonPair: mean energy per ion pair "
    192191           << " for Geant4 materials" << G4endl;
    193     for(G4int i=0; i<nMaterials; i++) {
     192    for(G4int i=0; i<nMaterials; ++i) {
    194193      G4cout << "   " << g4MatNames[i] << "    Epair= "
    195194             << g4MatData[i]/eV << " eV" << G4endl;

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmCalculator.cc,v 1.53 2010/04/13 10:58:03 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4EmCalculator.cc,v 1.57 2010/11/04 12:55:09 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9898  currentLambda      = 0;
    9999  currentModel       = 0;
     100  currentProcess     = 0;
    100101  loweModel          = 0;
    101102  chargeSquare       = 1.0;
    102103  massRatio          = 1.0;
     104  mass               = 0.0;
     105  currentCut         = 0.0;
    103106  currentParticleName= "";
    104107  currentMaterialName= "";

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmConfigurator.cc,v 1.8 2010/06/04 15:33:56 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4EmConfigurator.cc,v 1.9 2010/07/29 11:13:28 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    172172
    173173      G4VProcess* proc = 0;
    174       for(G4int i=0; i<np; i++) {
     174      for(G4int i=0; i<np; ++i) {
    175175        if(processName == (*plist)[i]->GetProcessName()) {
    176176          proc = (*plist)[i];
    … …  
    182182               << processName << "> for " << particleName << G4endl;
    183183        return;
    184 
    185184      }
    186185
    187186      if(mod) {
    188         if(!UpdateModelEnergyRange(mod, emin,emax)) { return; }
     187        if(!UpdateModelEnergyRange(mod, emin, emax)) { return; }
    189188      }
    190189      // classify process

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmModelManager.cc,v 1.60 2010/04/12 18:28:40 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4EmModelManager.cc,v 1.63 2010/10/15 10:22:13 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    463463      }
    464464    }
    465     /*
    466     G4int nm = setOfRegionModels[reg]->NumberOfModels();
    467     for(G4int j=0; j<nm; ++j) {
    468 
    469       G4VEmModel* model = models[setOfRegionModels[reg]->ModelIndex(j)];
    470       G4double tcutmin = model->MinEnergyCut(particle, couple);
    471 
    472       if(cut < tcutmin)    { cut = tcutmin; }
    473       if(subcut < tcutmin) { subcut = tcutmin; }
    474       if(1 < verboseLevel) {
    475             G4cout << "The model # " << j
    476                    << "; tcutmin(MeV)= " << tcutmin/MeV
    477                    << "; tcut(MeV)= " << cut/MeV
    478                    << "; tsubcut(MeV)= " << subcut/MeV
    479                    << " for " << particle->GetParticleName()
    480                    << G4endl;
    481       }
    482     }
    483     */
    484465    theCuts[i] = cut;
    485466    if(minSubRange < 1.0) { theSubCuts[i] = subcut; }
    … …  
    656637void G4EmModelManager::DumpModelList(G4int verb)
    657638{
    658   if(verb == 0) return;
     639  if(verb == 0) { return; }
    659640  for(G4int i=0; i<nRegions; ++i) {
    660641    G4RegionModels* r = setOfRegionModels[i];
    … …  
    665646             << " ======" << G4endl;;
    666647      for(G4int j=0; j<n; ++j) {
    667         const G4VEmModel* m = models[r->ModelIndex(j)];
     648        G4VEmModel* m = models[r->ModelIndex(j)];
    668649        G4cout << std::setw(20);
    669         G4cout << m->GetName() << " :     Emin= "
    670                << std::setw(10) << G4BestUnit(r->LowEdgeEnergy(j),"Energy")
    671                << "        Emax=   "
    672                << G4BestUnit(r->LowEdgeEnergy(j+1),"Energy")
    673                << G4endl;
     650        G4cout << m->GetName() << " :   Emin= "
     651               << std::setw(8) << G4BestUnit(r->LowEdgeEnergy(j),"Energy")
     652               << "      Emax=   "
     653               << std::setw(8) << G4BestUnit(r->LowEdgeEnergy(j+1),"Energy");
     654        G4VEmAngularDistribution* an = m->GetAngularDistribution();
     655        if(an) { G4cout << "  " << an->GetName(); }
     656        G4cout << G4endl;
    674657      } 
    675658    }

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmMultiModel.cc,v 1.6 2007/05/22 17:31:58 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4EmMultiModel.cc,v 1.8 2010/08/17 17:36:59 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4040// Modifications:
    4141// 15-04-05 optimize internal interface (V.Ivanchenko)
    42 //
    43 
    44 // Class Description:
    45 //
    46 // Energy loss model using several G4VEmModels
    47 
    48 // -------------------------------------------------------------------
     42// 04-07-10 updated interfaces according to g4 9.4 (V.Ivanchenko)
    4943//
    5044
    … …  
    5953
    6054G4EmMultiModel::G4EmMultiModel(const G4String& nam)
    61   : G4VEmModel(nam),
    62   nModels(0)
     55  : G4VEmModel(nam), nModels(0)
    6356{
    6457  model.clear();
    65   tsecmin.clear();
    6658  cross_section.clear();
    6759}
    … …  
    7062
    7163G4EmMultiModel::~G4EmMultiModel()
     64{}
     65
     66//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     67
     68void G4EmMultiModel::AddModel(G4VEmModel* p)
    7269{
    73   if(nModels) {
    74     for(G4int i=0; i<nModels; i++) {
    75       delete model[i];
    76     }
    77   }
     70  cross_section.push_back(0.0);
     71  model.push_back(p);
     72  ++nModels;
    7873}
    7974
    … …  
    8378                                const G4DataVector& cuts)
    8479{
    85   if(nModels) {
    86     for(G4int i=0; i<nModels; i++) {
     80  if(nModels > 0) {
     81    G4cout << "### Initialisation of EM MultiModel " << GetName()
     82           << " including following list of models:" << G4endl;
     83    for(G4int i=0; i<nModels; ++i) {
     84      G4cout << "    " << (model[i])->GetName();
     85      (model[i])->SetParticleChange(pParticleChange, GetModelOfFluctuations());
    8786      (model[i])->Initialise(p, cuts);
    8887    }
     88    G4cout << G4endl;
    8989  }
    90 }
    91 
    92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    93 
    94 G4double G4EmMultiModel::MinEnergyCut(const G4ParticleDefinition* p,
    95                                       const G4MaterialCutsCouple* couple)
    96 {
    97   G4double cut = DBL_MAX;
    98   if(nModels) {
    99     cut = (model[0])->MinEnergyCut(p, couple);
    100   }
    101   return cut;
    10290}
    10391
    … …  
    10694G4double G4EmMultiModel::ComputeDEDX(const G4MaterialCutsCouple* couple,
    10795                                     const G4ParticleDefinition* p,
    108                                            G4double kineticEnergy,
    109                                            G4double cutEnergy)
     96                                     G4double kineticEnergy,
     97                                     G4double cutEnergy)
    11098{
    111   G4double dedx  = 0.0;
     99  SetCurrentCouple(couple);
     100  G4double dedx = 0.0;
    112101
    113   if(nModels) {
    114     dedx =  (model[0])->ComputeDEDX(couple, p, cutEnergy, kineticEnergy);
     102  if(nModels > 0) {
     103    for(G4int i=0; i<nModels; i++) {
     104      dedx += (model[i])->ComputeDEDX(couple, p, cutEnergy, kineticEnergy);
     105    }
    115106  }
    116107
    … …  
    120111//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    121112
    122 G4double G4EmMultiModel::CrossSection(const G4MaterialCutsCouple* couple,
    123                                       const G4ParticleDefinition* p,
    124                                             G4double kineticEnergy,
    125                                             G4double cutEnergy,
    126                                             G4double maxKinEnergy)
     113G4double G4EmMultiModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
     114                                                    G4double kinEnergy,
     115                                                    G4double Z,
     116                                                    G4double A,
     117                                                    G4double cutEnergy,
     118                                                    G4double maxEnergy)
    127119{
    128   G4double cross     = 0.0;
    129   G4double t1        = cutEnergy;
    130   G4double t2        = cutEnergy;
    131   if(nModels) {
    132     for(G4int i=0; i<nModels; i++) {
    133       t1 = std::max(t2, tsecmin[i]);
    134       t2 = std::min(maxKinEnergy, tsecmin[i+1]);
    135       cross += (model[i])->CrossSection(couple, p, kineticEnergy, t1, t2);
     120  G4double cross = 0.0;
     121  if(nModels>0) {
     122    for(G4int i=0; i<nModels; ++i) {
     123      (model[i])->SetCurrentCouple(CurrentCouple());
     124      cross += (model[i])->ComputeCrossSectionPerAtom(p, kinEnergy, Z, A,
     125                                                      cutEnergy, maxEnergy);
    136126    }
    137127  }
    … …  
    144134                                       const G4MaterialCutsCouple* couple,
    145135                                       const G4DynamicParticle* dp,
    146                                        G4double tmin,
     136                                       G4double minEnergy,
    147137                                       G4double maxEnergy)
    148138{
     139  SetCurrentCouple(couple);
    149140  if(nModels > 0) {
    150141    G4int i;
    151142    G4double cross = 0.0;
    152     G4double t1    = tmin;
    153     G4double t2    = tmin;
    154     for(i=0; i<nModels; i++) {
    155       t1 = std::max(t2, tsecmin[i]);
    156       t2 = std::min(maxEnergy, tsecmin[i+1]);
    157       cross += (model[i])->CrossSection(couple, dp->GetDefinition(),
    158                                         dp->GetKineticEnergy(), t1, t2);
     143    for(i=0; i<nModels; ++i) {
     144      cross += (model[i])->CrossSection(couple, dp->GetParticleDefinition(),
     145                                        dp->GetKineticEnergy(), minEnergy, maxEnergy);
    159146      cross_section[i] = cross;
    160147    }
    161148
    162149    cross *= G4UniformRand();
    163     t2 = tmin;
    164150
    165     for(i=0; i<nModels; i++) {
    166       t1 = std::max(t2, tsecmin[i]);
    167       t2 = std::min(maxEnergy, tsecmin[i+1]);
     151    for(i=0; i<nModels; ++i) {
    168152      if(cross <= cross_section[i]) {
    169         (model[i])->SampleSecondaries(vdp, couple, dp, t1, t2);
    170         break;
     153        (model[i])->SampleSecondaries(vdp, couple, dp, minEnergy, maxEnergy);
     154        return;
    171155      }
    172156    }
    … …  
    176160//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    177161
    178 G4double G4EmMultiModel:: MaxSecondaryEnergy(const G4ParticleDefinition*,
    179                                                    G4double kinEnergy)
    180 {
    181   G4cout << "Warning! G4EmMultiModel::"
    182          << "MaxSecondaryEnergy(const G4ParticleDefinition*,G4double kinEnergy)"
    183          << " should not be used!" << G4endl;
    184   return kinEnergy;
    185 }
    186 
    187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    188 
    189 void G4EmMultiModel::DefineForRegion(const G4Region* r)
    190 {
    191   if(nModels) {
    192     for(G4int i=0; i<nModels; i++) {(model[i])->DefineForRegion(r);}
    193   }
    194 }
    195 
    196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    197 
    198 void G4EmMultiModel::AddModel(G4VEmModel* p, G4double tmin, G4double tmax)
    199 {
    200   if(tmin < tmax && 0.0 < tmin) {
    201 
    202     if(nModels == 0) {
    203       tsecmin.push_back(tmin);
    204       tsecmin.push_back(tmax);
    205       cross_section.push_back(0.0);
    206       model.push_back(p);
    207       nModels++;
    208 
    209     } else {
    210       G4int i, j;
    211       G4bool increment = false;
    212       for(i=0; i<nModels; i++) {
    213 
    214         if(tmin < tsecmin[i]) {
    215           G4double t2 = std::min(tsecmin[i+1],tmax);
    216           if(tmin < t2) {
    217             tsecmin.push_back(0.0);
    218             cross_section.push_back(0.0);
    219             model.push_back(0);
    220             for(j=nModels; j>i; j--) {
    221               model[j] = model[j-1];
    222               tsecmin[j+1] = tsecmin[j];
    223             }
    224             model[i] = p;
    225             tsecmin[i+1] = t2;
    226             tsecmin[i]   = tmin;
    227             increment = true;
    228           }
    229         } else if(i == nModels-1) {
    230           G4double t1 = std::min(tsecmin[i+1],tmin);
    231           G4double t2 = std::max(tsecmin[i+1],tmax);
    232           if(t1 < t2) {
    233             tsecmin.push_back(t2);
    234             cross_section.push_back(0.0);
    235             model.push_back(p);
    236             increment = true;
    237           }
    238         }
    239       }
    240       if(increment) nModels++;
    241     }
    242   }
    243 }
    244 
    245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    246 

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EmSaturation.cc,v 1.10 2009/09/25 09:16:40 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4EmSaturation.cc,v 1.11 2010/10/25 17:23:01 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5959  verbose = 1;
    6060  manager = 0;
     61
    6162  curMaterial = 0;
    62   curBirks = 0.0;
    63   curRatio = 1.0;
     63  curBirks    = 0.0;
     64  curRatio    = 1.0;
    6465  curChargeSq = 1.0;
    65   nMaterials = 0;
     66  nMaterials  = 0;
     67
    6668  electron = 0;
     69  proton   = 0;
     70  nist     = G4NistManager::Instance();
     71
    6772  Initialise();
    6873}
    … …  
    8287                                      G4double niel)
    8388{
    84   if(edep <= 0.0) return 0.0;
     89  if(edep <= 0.0) { return 0.0; }
    8590
    8691  G4double evis = edep;
    … …  
    109114
    110115      // continues energy loss
    111       if(eloss > 0.0) eloss /= (1.0 + bfactor*eloss/length);
     116      if(eloss > 0.0) { eloss /= (1.0 + bfactor*eloss/length); }
    112117 
    113118      // non-ionizing energy loss
    114119      if(nloss > 0.0) {
    115         if(!proton) {proton = G4Proton::Proton();}
     120        if(!proton) { proton = G4Proton::Proton(); }
    116121        G4double escaled = nloss*curRatio;
    117122        G4double s = manager->GetRange(proton,escaled,couple)/curChargeSq;
    … …  
    133138  // is this material in the vector?
    134139 
    135   for(G4int j=0; j<nG4Birks; j++) {
     140  for(G4int j=0; j<nG4Birks; ++j) {
    136141    if(name == g4MatNames[j]) {
    137142      if(verbose > 0)
    … …  
    152157  if(!manager) {
    153158    manager = G4LossTableManager::Instance();
    154     nist    = G4NistManager::Instance();
    155159    electron= G4Electron::Electron();
    156     proton  = 0;
    157   }
    158 
    159   if(mat == curMaterial) return curBirks;
     160  }
     161
     162  if(mat == curMaterial) { return curBirks; }
    160163
    161164  curMaterial = mat;
    … …  
    165168
    166169  // seach in the run-time list
    167   for(G4int i=0; i<nMaterials; i++) {
     170  for(G4int i=0; i<nMaterials; ++i) {
    168171    if(mat == matPointers[i]) {
    169172      curBirks = mat->GetIonisation()->GetBirksConstant();
    … …  
    180183  // seach in the Geant4 list
    181184  if(curBirks == 0.0) {
    182     for(G4int j=0; j<nG4Birks; j++) {
     185    for(G4int j=0; j<nG4Birks; ++j) {
    183186      if(name == g4MatNames[j]) {
    184187        mat->GetIonisation()->SetBirksConstant(g4MatData[j]);
    … …  
    201204  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
    202205  size_t nelm = mat->GetNumberOfElements();
    203   for (size_t i=0; i<nelm; i++) {
     206  for (size_t i=0; i<nelm; ++i) {
    204207    const G4Element* elm = (*theElementVector)[i];
    205208    G4double Z = elm->GetZ();
    … …  
    231234  if(nMaterials > 0) {
    232235    G4cout << "### Birks coeffitients used in run time" << G4endl;
    233     for(G4int i=0; i<nMaterials; i++) {
     236    for(G4int i=0; i<nMaterials; ++i) {
    234237      G4double br = matPointers[i]->GetIonisation()->GetBirksConstant();
    235238      G4cout << "   " << matNames[i] << "     "
    … …  
    248251  if(nG4Birks > 0) {
    249252    G4cout << "### Birks coeffitients for Geant4 materials" << G4endl;
    250     for(G4int i=0; i<nG4Birks; i++) {
     253    for(G4int i=0; i<nG4Birks; ++i) {
    251254      G4cout << "   " << g4MatNames[i] << "   "
    252255             << g4MatData[i]*MeV/mm << " mm/MeV" << G4endl;

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LossTableManager.cc,v 1.101 2010/06/04 15:33:56 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4LossTableManager.cc,v 1.105 2010/11/04 12:55:09 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9595#include "G4EmSaturation.hh"
    9696#include "G4EmConfigurator.hh"
     97#include "G4ElectronIonPair.hh"
    9798#include "G4EmTableType.hh"
    9899#include "G4LossTableBuilder.hh"
     100#include "G4VAtomDeexcitation.hh"
    99101#include "G4Region.hh"
    100102
    … …  
    140142  delete emCorrections;
    141143  delete emSaturation;
     144  delete emConfigurator;
     145  delete emElectronIonPair;
     146  delete atomDeexcitation;
    142147}
    143148
    … …  
    182187  emSaturation = new G4EmSaturation();
    183188  emConfigurator = new G4EmConfigurator(verbose);
     189  emElectronIonPair = new G4ElectronIonPair();
    184190  tableBuilder->SetSplineFlag(splineFlag);
     191  atomDeexcitation = 0;
    185192}
    186193
    … …  
    438445    emConfigurator->Clear();
    439446    firstParticle = aParticle;
     447  }
     448  if(startInitialisation && atomDeexcitation) {
     449    atomDeexcitation->InitialiseAtomicDeexcitation();
    440450  }
    441451  startInitialisation = false;
    … …  
    590600
    591601  G4int n_dedx = t_list.size();
    592   if (!n_dedx) {
     602  if (0 == n_dedx || !em) {
    593603    G4cout << "G4LossTableManager WARNING: no DEDX processes for "
    594604           << aParticle->GetParticleName() << G4endl;
    … …  
    896906  //emCorrections->SetVerbose(val);
    897907  emSaturation->SetVerbose(val);
     908  emElectronIonPair->SetVerbose(val);
     909  if(atomDeexcitation) { atomDeexcitation->SetVerboseLevel(val); }
    898910}
    899911
    … …  
    10441056}
    10451057
     1058//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     1059
     1060G4ElectronIonPair* G4LossTableManager::ElectronIonPair()
     1061{
     1062  return emElectronIonPair;
     1063}
     1064
    10461065//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1066
     1067G4VAtomDeexcitation* G4LossTableManager::AtomDeexcitation()
     1068{
     1069  return atomDeexcitation;
     1070}
     1071
     1072//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1073 
     1074void G4LossTableManager::SetAtomDeexcitation(G4VAtomDeexcitation* p)
     1075{
     1076  atomDeexcitation = p;
     1077}
     1078//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmModel.cc,v 1.35 2010/05/26 18:06:25 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VEmModel.cc,v 1.37 2010/10/14 16:27:35 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6161
    6262G4VEmModel::G4VEmModel(const G4String& nam):
    63   fluc(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV),
     63  flucModel(0),anglModel(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV),
    6464  eMinActive(0.0),eMaxActive(DBL_MAX),
    6565  polarAngleLimit(0.0),secondaryThreshold(DBL_MAX),theLPMflag(false),
    66   pParticleChange(0),nuclearStopping(false),
     66  pParticleChange(0),/*nuclearStopping(false),*/
    6767  currentCouple(0),currentElement(0),
    6868  nsec(5),flagDeexcitation(false)
    … …  
    8484    }
    8585  }
     86  delete anglModel;
    8687}
    8788
    … …  
    246247G4double G4VEmModel::ChargeSquareRatio(const G4Track& track)
    247248{
    248   return GetChargeSquareRatio(track.GetDefinition(),
     249  return GetChargeSquareRatio(track.GetParticleDefinition(),
    249250                              track.GetMaterial(), track.GetKineticEnergy());
    250251}
    … …  
    296297
    297298//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     299
     300void
     301G4VEmModel::SetParticleChange(G4VParticleChange* p, G4VEmFluctuationModel* f)
     302{
     303  if(p && pParticleChange != p) { pParticleChange = p; }
     304  flucModel = f;
     305}
     306
     307//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEmProcess.cc,v 1.87 2010/05/10 13:35:32 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VEmProcess.cc,v 1.88 2010/08/17 17:36:59 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    516516    for (G4int i=0; i<num; ++i) {
    517517      G4DynamicParticle* dp = secParticles[i];
    518       const G4ParticleDefinition* p = dp->GetDefinition();
     518      const G4ParticleDefinition* p = dp->GetParticleDefinition();
    519519      G4double e = dp->GetKineticEnergy();
    520520      G4bool good = true;

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEnergyLossProcess.cc,v 1.167 2010/05/26 10:41:34 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VEnergyLossProcess.cc,v 1.171 2010/10/15 10:22:13 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    108108// 27-10-07 Virtual functions moved to source (V.Ivanchenko)
    109109// 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
     110// 15-10-10 Fixed 4-momentum balance if deexcitation is active (L.Pandola)
    110111//
    111112// Class Description:
    … …  
    144145#include "G4TransportationManager.hh"
    145146#include "G4EmConfigurator.hh"
     147#include "G4VAtomDeexcitation.hh"
    146148
    147149//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    222224  (G4LossTableManager::Instance())->Register(this);
    223225  fluctModel = 0;
     226  atomDeexcitation = 0;
    224227
    225228  scTracks.reserve(5);
    … …  
    570573
    571574  // Added tracking cut to avoid tracking artifacts
    572   if(isIonisation) { fParticleChange.SetLowEnergyLimit(lowestKinEnergy); }
     575  if(isIonisation) {
     576    fParticleChange.SetLowEnergyLimit(lowestKinEnergy);
     577    atomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
     578    if(atomDeexcitation) { useDeexcitation = true; }
     579  }
    573580
    574581  if(1 < verboseLevel) {
    … …  
    852859    for (G4int i=0; i<nDERegions; ++i) {
    853860      if (reg == deRegions[i]) {
    854         if(!val) deRegions[i] = 0;
     861        if(!val) { deRegions[i] = 0; }
    855862        return;
    856863      }
    … …  
    911918  DefineMaterial(track.GetMaterialCutsCouple());
    912919
    913   const G4ParticleDefinition* currPart = track.GetDefinition();
     920  const G4ParticleDefinition* currPart = track.GetParticleDefinition();
    914921  if(theGenericIon == particle) {
    915922    massRatio = proton_mass_c2/currPart->GetPDGMass();
    … …  
    9991006  /* 
    10001007  if(-1 < verboseLevel) {
    1001     const G4ParticleDefinition* d = track.GetDefinition();
     1008    const G4ParticleDefinition* d = track.GetParticleDefinition();
    10021009    G4cout << "AlongStepDoIt for "
    10031010           << GetProcessName() << " and particle "
    … …  
    10191026    eloss = preStepKinEnergy;
    10201027    if (useDeexcitation) {
    1021       if(idxDERegions[currentMaterialIndex]) {
     1028      if(atomDeexcitation) {
     1029        atomDeexcitation->AlongStepDeexcitation(&fParticleChange, step,
     1030                                                eloss, currentMaterialIndex);
     1031      } else if(idxDERegions[currentMaterialIndex]) {
    10221032        currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
    1023         if(eloss < 0.0) eloss = 0.0;
    1024       }
     1033      }
     1034      if(eloss < 0.0) { eloss = 0.0; }
    10251035    }
    10261036    fParticleChange.SetProposedKineticEnergy(0.0);
    … …  
    11281138            G4Track* t = scTracks[i];
    11291139            G4double e = t->GetKineticEnergy();
    1130             if (t->GetDefinition() == thePositron) { e += 2.0*electron_mass_c2; }
     1140            if (t->GetParticleDefinition() == thePositron) {
     1141              e += 2.0*electron_mass_c2;
     1142            }
    11311143            esec += e;
    11321144            pParticleChange->AddSecondary(t);
    … …  
    11691181  // deexcitation
    11701182  if (useDeexcitation) {
    1171     if(idxDERegions[currentMaterialIndex] && eloss > 0.0) {
     1183    G4double eloss_before = eloss;
     1184    if(atomDeexcitation) {
     1185      atomDeexcitation->AlongStepDeexcitation(&fParticleChange, step,
     1186                                              eloss, currentMaterialIndex);
     1187    } else if(idxDERegions[currentMaterialIndex] && eloss > 0.0) {
    11721188      currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
    1173       if(eloss < 0.0) { eloss = 0.0; }
    1174     }
     1189    }
     1190    esec += eloss_before - eloss;
    11751191  }
    11761192
    … …  
    11781194  G4double finalT = preStepKinEnergy - eloss - esec;
    11791195  if (finalT <= lowestKinEnergy) {
    1180     eloss  = preStepKinEnergy - esec;
     1196    eloss += finalT;
    11811197    finalT = 0.0;
    11821198  } else if(isIon) {
    11831199    fParticleChange.SetProposedCharge(
    1184       currentModel->GetParticleCharge(track.GetDefinition(),currentMaterial,finalT));
    1185   }
    1186 
     1200      currentModel->GetParticleCharge(track.GetParticleDefinition(),
     1201                                      currentMaterial,finalT));
     1202  }
     1203
     1204  if(eloss < 0.0) { eloss = 0.0; }
    11871205  fParticleChange.SetProposedKineticEnergy(finalT);
    11881206  fParticleChange.ProposeLocalEnergyDeposit(eloss);
    … …  
    12601278      /*
    12611279      // do not track very low-energy delta-electrons
    1262       if(theSecondaryRangeTable && (*it)->GetDefinition() == theElectron) {
     1280      if(theSecondaryRangeTable && (*it)->GetParticleDefinition() == theElectron) {
    12631281        G4double ekin = (*it)->GetKineticEnergy();
    12641282        G4double rg = ((*theSecondaryRangeTable)[idx]->Value(ekin));
    … …  
    12781296        /*     
    12791297        if(-1 < verboseLevel)
    1280           G4cout << "New track " << t->GetDefinition()->GetParticleName()
     1298          G4cout << "New track " << t->GetParticleDefinition()->GetParticleName()
    12811299                 << " e(keV)= " << t->GetKineticEnergy()/keV
    12821300                 << " fragment= " << fragment
    … …  
    12931311                                                      const G4Step&)
    12941312{
     1313  // In all cases clear number of interaction lengths
     1314  theNumberOfInteractionLengthLeft = -1.0;
     1315
    12951316  fParticleChange.InitializeForPostStep(track);
    12961317  G4double finalT = track.GetKineticEnergy();
    1297   if(finalT <= lowestKinEnergy) return &fParticleChange;
     1318  if(finalT <= lowestKinEnergy) { return &fParticleChange; }
    12981319
    12991320  G4double postStepScaledEnergy = finalT*massRatio;
    … …  
    13211342    }
    13221343    */
    1323     if(preStepLambda*G4UniformRand() > lx) {
    1324       ClearNumberOfInteractionLengthLeft();
     1344    if(lx <= 0.0) {
     1345      return &fParticleChange;
     1346    } else if(preStepLambda*G4UniformRand() > lx) {
    13251347      return &fParticleChange;
    13261348    }
    … …  
    13601382  }
    13611383  */
    1362   ClearNumberOfInteractionLengthLeft();
    13631384  return &fParticleChange;
    13641385}

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMscModel.cc,v 1.16 2010/04/06 09:24:46 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VMscModel.cc,v 1.18 2010/09/07 16:05:33 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6464  dtrl(0.05),
    6565  lambdalimit(mm),
    66   geommax(1.e50*mm),
     66  geomMin(1.e-6*CLHEP::mm),
     67  geomMax(1.e50*CLHEP::mm),
    6768  steppingAlgorithm(fUseSafety),
    6869  samplez(false),
    … …  
    106107                                      G4double postsafety)
    107108{
     109  if(displacement <= geomMin) { return; }
    108110  const G4ThreeVector* pos = fParticleChange->GetProposedPosition();
    109   G4double r = displacement;
    110   if(r >  postsafety) {
    111     G4double newsafety = safetyHelper->ComputeSafety(*pos);
    112     if(r > newsafety) r = newsafety;
     111
     112  // displaced point is definitely within the volume
     113  if(displacement < postsafety) {
     114
     115    // compute new endpoint of the Step
     116    G4ThreeVector newPosition = *pos + displacement*dir;
     117    safetyHelper->ReLocateWithinVolume(newPosition);
     118    fParticleChange->ProposePosition(newPosition);
     119    return;
    113120  }
    114   if(r > 0.) {
     121
     122  // displaced point may be outside the volume
     123  G4double newsafety = safetyHelper->ComputeSafety(*pos);
     124
     125  // add a factor which ensure numerical stability
     126  G4double r = std::min(displacement, newsafety*0.99);
     127
     128  if(r > geomMin) {
    115129
    116130    // compute new endpoint of the Step
    117131    G4ThreeVector newPosition = *pos + r*dir;
    118 
    119     // definitely not on boundary
    120     if(displacement == r) {
    121       safetyHelper->ReLocateWithinVolume(newPosition);
    122 
    123     } else {
    124       // check safety after displacement
    125       G4double postsafety = safetyHelper->ComputeSafety(newPosition);
    126 
    127       // displacement to boundary
    128       if(postsafety <= 0.0) {
    129         safetyHelper->Locate(newPosition,
    130                              *fParticleChange->GetProposedMomentumDirection());
    131 
    132         // not on the boundary
    133       } else {
    134         safetyHelper->ReLocateWithinVolume(newPosition);
    135       }
    136     }
     132    safetyHelper->ReLocateWithinVolume(newPosition);
    137133    fParticleChange->ProposePosition(newPosition);
    138134  }

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VMultipleScattering.cc,v 1.82 2010/04/12 11:45:03 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4VMultipleScattering.cc,v 1.86 2010/10/26 11:30:46 vnivanch Exp $
     27// GEANT4 tag $Name: emutils-V09-03-23 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    294294  if(verboseLevel>0 && ( num == "e-" || num == "mu+" || 
    295295                         num == "proton" || num == "pi-" ||
    296                          num == "GenericIon")) {
     296                         num == "kaon+" || num == "GenericIon")) {
    297297    PrintInfoDefinition();
    298298    if(2 < verboseLevel && theLambdaTable) G4cout << *theLambdaTable << G4endl;
    … …  
    348348  DefineMaterial(track.GetMaterialCutsCouple());
    349349  G4double ekin = track.GetKineticEnergy();
    350   if(isIon) { ekin *= proton_mass_c2/track.GetDefinition()->GetPDGMass(); }
     350  if(isIon) { ekin *= proton_mass_c2/track.GetParticleDefinition()->GetPDGMass(); }
    351351  currentModel = static_cast<G4VMscModel*>(SelectModel(ekin));
    352352  if(x > 0.0 && ekin > 0.0 && currentModel->IsActive(ekin)) {
    353353    G4double tPathLength =
    354354      currentModel->ComputeTruePathLengthLimit(track, theLambdaTable, x);
    355     if (tPathLength < x) *selection = CandidateForSelection;
     355    if (tPathLength < x) { *selection = CandidateForSelection; }
    356356    x = currentModel->ComputeGeomPathLength(tPathLength);
    357357    //  G4cout << "tPathLength= " << tPathLength
    … …  
    404404                                       G4double& currentSafety)
    405405{
    406   G4GPILSelection* selection = 0;
    407   return AlongStepGetPhysicalInteractionLength(track,previousStepSize,currentMinimalStep,
    408                                                currentSafety, selection);
     406  G4GPILSelection selection = NotCandidateForSelection;
     407  G4double x = AlongStepGetPhysicalInteractionLength(track,previousStepSize,
     408                                                     currentMinimalStep,
     409                                                     currentSafety, &selection);
     410  return x;
    409411}
    410412
Note: See TracChangeset for help on using the changeset viewer.