Ignore:
Timestamp:
Apr 6, 2009, 12:21:12 PM (15 years ago)
Author:
garnier
Message:

update processes

Location:
trunk/source/processes/electromagnetic/muons
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/electromagnetic/muons/History

    r819 r961  
    1 $Id: History,v 1.106 2007/11/12 10:34:23 vnivanch Exp $
     1$Id: History,v 1.126 2009/02/26 11:04:20 vnivanch Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     2026 February 09: V.Ivant (emmuons-V09-02-01)
     21G4MuIonisation - fixed initialisation alowing to configure external model
     22                 of fluctuations
     23
     2420 February 09: V.Ivant (emmuons-V09-02-00)
     25- Cleanup: improved comments, move virtual methods from .hh to .cc
     26
     2712 November 08: V.Ivant (emmuons-V09-01-15)
     28G4EnergyLossForExtrapolator - added method TrueStepLength; fixed initialisation
     29                              before a step
     30
     3127 October 08: V.Ivant (emmuons-V09-01-14)
     32G4EnergyLossForExtrapolator - make method ComputeTrueStep public and cleanup
     33
     3416 October 08: V.Ivant (emmuons-V09-01-13)
     35G4MuMscModel - remove obsolete
     36G4EnergyLossForExtrapolator - added spline option
     37G4MuIonisation, G4MuBremsstrahlung, G4MuPairProduction,
     38G4MuMultipleScattering - change SubType and improved cout
     39
     404 August 08: V.Ivant (emmuons-V09-01-12)
     41G4MuMscModel - added protection for ions
     42
     4331 July 08: V.Ivant (emmuons-V09-01-11)
     44G4MuMscModel - do not define min and max energy in constructor but use Set
     45               methods
     46G4MuMultipleScattering - added cout of model names
     47
     4821 April 08:  V.Ivanchenko (emmuons-V09-01-10)
     49G4MuBremsstrahlungModel, G4MuPairProductionModel - use CrossSectionPerVolume
     50                    from the base class, do not use A in CrossSEctionPerAtom
     51G4MuMscModel - do not use A in SetupTarget
     52
     5304 April 08:  V.Ivanchenko (emmuons-V09-01-09)
     54G4MuMultipleScattering - use G4WentzelVIModel model
     55                         build table for particles with mass < GeV
     56
     5704 April 08:  V.Ivanchenko (emmuons-V09-01-08)
     58- G4MuBremsstrahlungModel - instead of static const use members of a class,
     59                            this allows to reuse the model for different
     60                            particle type
     61
     6227 March 08:  V.Ivanchenko (emmuons-V09-01-07)
     63- G4MuPairProductionModel - fixed nan value at initialisation
     64  of the sampling table
     65
     6626 March 08:  V.Ivanchenko (emmuons-V09-01-06)
     67- G4MuMscModel - fixed outstanding bug in sampling of scattering
     68
     6925 March 08:  V.Ivanchenko (emmuons-V09-01-05)
     70- G4MuMscModel - added shift along particle direction for displacement
     71- G4MuBetheBlochModel - update computation of correction
     72
     7317 March 08:  V.Ivanchenko (emmuons-V09-01-04)
     74- G4MuMscModel - fixed sampling
     75
     7614 March 08:  V.Ivanchenko (emmuons-V09-01-03)
     77- G4MuMscModel - use G4VMscModel interface
     78
     7906 March 08:  V.Ivanchenko (emmuons-V09-01-02)
     80- G4MuBremsstrahlungModel - remove ignoreCut flag, remove obsolete methods
     81                            and members, set some members protected to
     82                            be used by G4hBremsstrahlungModel
     83- G4MuPairProductionModel - remove ignoreCut flag,  set some members protected
     84                            to be used by G4hBremsstrahlungModel
     85- SubType for all processes is initialized
     86
     8722 February 08:  V.Ivanchenko (emmuons-V09-01-01)
     88G4MuMscModel - added sampling of tail distribution
     89
     9014 January 08:  V.Ivanchenko (emmuons-V09-01-00)
     91G4MuMscModel - added computation of the second moment of the distribution;
     92               fixed sampling
     93G4MuMultipleScattering - modified default RangeFactor
    1994
    209512 November 07:  V.Ivanchenko (emmuons-V09-00-04)

  • trunk/source/processes/electromagnetic/muons/include/G4EnergyLossForExtrapolator.hh

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EnergyLossForExtrapolator.hh,v 1.9 2007/07/28 13:44:25 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: G4EnergyLossForExtrapolator.hh,v 1.12 2008/11/13 14:14:07 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929//---------------------------------------------------------------------------
    … …  
    8181
    8282  G4double EnergyBeforeStep(G4double kinEnergy, G4double step,
    83                             const G4Material*, const G4ParticleDefinition*);
     83                            const G4Material*, const G4ParticleDefinition*);
     84
     85  G4double TrueStepLength(G4double kinEnergy, G4double step,
     86                          const G4Material*, const G4ParticleDefinition* part);
    8487
    8588  inline G4double EnergyAfterStep(G4double kinEnergy, G4double step,
    86                            const G4Material*, const G4String& particleName);
     89                                  const G4Material*, const G4String& particleName);
    8790
    8891  inline G4double EnergyBeforeStep(G4double kinEnergy, G4double step,
    89                             const G4Material*, const G4String& particleName);
     92                                   const G4Material*, const G4String& particleName);
    9093
    9194  inline G4double AverageScatteringAngle(G4double kinEnergy, G4double step,
    92                                   const G4Material*, const G4ParticleDefinition* part);
     95                                         const G4Material*,
     96                                         const G4ParticleDefinition* part);
    9397
    9498  inline G4double AverageScatteringAngle(G4double kinEnergy, G4double step,
    95                                          const G4Material*, const G4String& particleName);
     99                                         const G4Material*,
     100                                         const G4String& particleName);
     101
     102  inline G4double ComputeTrueStep(const G4Material*, const G4ParticleDefinition* part,
     103                                  G4double kinEnergy, G4double stepLength);
    96104
    97105  inline G4double EnergyDispersion(G4double kinEnergy, G4double step,
    98                             const G4Material*, const G4ParticleDefinition*);
     106                                   const G4Material*, const G4ParticleDefinition*);
    99107
    100108  inline G4double EnergyDispersion(G4double kinEnergy, G4double step,
    … …  
    113121  void Initialisation();
    114122
     123  G4bool SetupKinematics(const G4ParticleDefinition*, const G4Material*,
     124                         G4double kinEnergy);
     125
    115126  G4PhysicsTable* PrepareTable();
    116127
    … …  
    123134  void ComputeProtonDEDX(const G4ParticleDefinition* part, G4PhysicsTable* table);
    124135
    125   G4double ComputeTrueStep(const G4Material*, const G4ParticleDefinition* part,
    126                            G4double kinEnergy, G4double stepLength);
     136  void ComputeTrasportXS(const G4ParticleDefinition* part, G4PhysicsTable* table);
    127137
    128138  inline G4double ComputeValue(G4double x, const G4PhysicsTable* table);
    129 
    130   inline G4double ComputeScatteringAngle(G4double x);
    131139
    132140  // hide assignment operator
    … …  
    158166  G4PhysicsTable*          invRangeMuon;
    159167  G4PhysicsTable*          invRangeProton;
     168  G4PhysicsTable*          mscElectron;
    160169
    161170  const G4Material* currentMaterial;
    … …  
    214223//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    215224
    216 inline G4double G4EnergyLossForExtrapolator::EnergyDispersion(G4double kinEnergy,
    217                                                               G4double step,
    218                                                               const G4Material* mat,
    219                                                               const G4String& name)
     225inline G4double
     226G4EnergyLossForExtrapolator::EnergyDispersion(G4double kinEnergy,
     227                                              G4double step,
     228                                              const G4Material* mat,
     229                                              const G4String& name)
    220230{
    221231  return EnergyDispersion(kinEnergy,step,mat,FindParticle(name));
    … …  
    224234//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    225235
    226 inline G4double G4EnergyLossForExtrapolator::AverageScatteringAngle(G4double kinEnergy,
    227                                                              G4double stepLength,
    228                                                              const G4Material* mat,
    229                                                              const G4ParticleDefinition* part)
    230 {
    231   if(!isInitialised) Initialisation();
     236inline G4double
     237G4EnergyLossForExtrapolator::AverageScatteringAngle(G4double kinEnergy,
     238                                                    G4double stepLength,
     239                                                    const G4Material* mat,
     240                                                    const G4ParticleDefinition* part)
     241{
    232242  G4double theta = 0.0;
    233   if(mat && part && kinEnergy > 0.0) {
    234     G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
    235     if(step > 0.001*radLength) theta = ComputeScatteringAngle(stepLength);
     243  if(SetupKinematics(part, mat, kinEnergy)) {
     244    G4double t = stepLength/radLength;
     245    G4double y = std::max(0.001, t);
     246    theta = 19.23*MeV*std::sqrt(charge2*t)*(1.0 + 0.038*std::log(y))/(beta2*gam*mass);
    236247  }
    237248  return theta;
    … …  
    240251//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    241252
    242 inline G4double G4EnergyLossForExtrapolator::ComputeScatteringAngle(G4double x)
    243 {
    244   G4double t = x/radLength;
    245   return 19.23*MeV*std::sqrt(charge2*t)*(1.0 + 0.038*std::log(t))/(beta2*gam*mass);
    246 }
     253inline G4double
     254G4EnergyLossForExtrapolator::ComputeTrueStep(const G4Material* mat,
     255                                             const G4ParticleDefinition* part,
     256                                             G4double kinEnergy,
     257                                             G4double stepLength)
     258{
     259  G4double theta = AverageScatteringAngle(kinEnergy,stepLength,mat,part);
     260  return stepLength*std::sqrt(1.0 + 0.625*theta*theta);
     261}
    247262
    248263//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    249264 
    250 inline G4double G4EnergyLossForExtrapolator::EnergyDispersion(
    251                                                        G4double kinEnergy,
     265inline
     266G4double G4EnergyLossForExtrapolator::EnergyDispersion(G4double kinEnergy,
    252267                                                       G4double stepLength,
    253268                                                       const G4Material* mat,
    254269                                                       const G4ParticleDefinition* part)
    255270{
    256   if(!isInitialised) Initialisation();
    257271  G4double sig2 = 0.0;
    258   if(mat && part ) {
     272  if(SetupKinematics(part, mat, kinEnergy)) {
    259273    G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
    260     sig2 = (1.0/beta2 - 0.5)* twopi_mc2_rcl2 * tmax*step * electronDensity * charge2;
     274    sig2 = (1.0/beta2 - 0.5)*twopi_mc2_rcl2*tmax*step*electronDensity*charge2;
    261275  }
    262276  return sig2;
    … …  
    269283{
    270284  G4double res = 0.0;
    271   bool b;
    272   res = ((*table)[index])->GetValue(x, b);
     285  G4bool b;
     286  if(table) res = ((*table)[index])->GetValue(x, b);
    273287  return res;
    274288}

  • trunk/source/processes/electromagnetic/muons/include/G4MuBetheBlochModel.hh

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuBetheBlochModel.hh,v 1.17 2007/05/22 17:35:58 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: G4MuBetheBlochModel.hh,v 1.18 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    7676  virtual ~G4MuBetheBlochModel();
    7777
    78   void Initialise(const G4ParticleDefinition*, const G4DataVector&);
     78  virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
    7979
    80   G4double MinEnergyCut(const G4ParticleDefinition*,
    81                         const G4MaterialCutsCouple*);
     80  virtual G4double MinEnergyCut(const G4ParticleDefinition*,
     81                                const G4MaterialCutsCouple*);
    8282                       
    8383  virtual G4double ComputeCrossSectionPerElectron(
    … …  
    113113protected:
    114114
    115   G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
    116                                     G4double kinEnergy);
     115  virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
     116                                      G4double kinEnergy);
    117117
    118118private:
    119119
    120   void SetParticle(const G4ParticleDefinition* p);
     120  inline void SetParticle(const G4ParticleDefinition* p);
    121121
    122122  // hide assignment operator
    … …  
    141141};
    142142
    143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     143//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    144144
    145 inline G4double G4MuBetheBlochModel::MaxSecondaryEnergy(
    146           const G4ParticleDefinition*,
    147                 G4double kinEnergy)
     145inline void G4MuBetheBlochModel::SetParticle(const G4ParticleDefinition* p)
    148146{
    149   G4double tau  = kinEnergy/mass;
    150   G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
    151                   (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
    152   return tmax;
     147  if(!particle) {
     148    particle = p;
     149    mass = particle->GetPDGMass();
     150    massSquare = mass*mass;
     151    ratio = electron_mass_c2/mass;
     152  }
    153153}
    154154

  • trunk/source/processes/electromagnetic/muons/include/G4MuBremsstrahlung.hh

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuBremsstrahlung.hh,v 1.29 2007/05/23 08:49:32 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuBremsstrahlung.hh,v 1.31 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8686  virtual ~G4MuBremsstrahlung();
    8787
    88   G4bool IsApplicable(const G4ParticleDefinition& p);
     88  virtual G4bool IsApplicable(const G4ParticleDefinition& p);
    8989
    90   G4double MinPrimaryEnergy(const G4ParticleDefinition* p,
    91                             const G4Material*,
    92                             G4double cut);
     90  virtual G4double MinPrimaryEnergy(const G4ParticleDefinition* p,
     91                                    const G4Material*,
     92                                    G4double cut);
    9393
    9494  // Print out of the class parameters
    95   void PrintInfo();
     95  virtual void PrintInfo();
    9696
    9797protected:
    9898
    99   void InitialiseEnergyLossProcess(const G4ParticleDefinition*,
    100                                    const G4ParticleDefinition*);
     99  virtual void InitialiseEnergyLossProcess(const G4ParticleDefinition*,
     100                                           const G4ParticleDefinition*);
    101101
    102102private:
    … …  
    114114
    115115//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    117 
    118 inline G4bool G4MuBremsstrahlung::IsApplicable(const G4ParticleDefinition& p)
    119 {
    120   return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV);
    121 }
    122 
    123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    124 
    125 inline G4double G4MuBremsstrahlung::MinPrimaryEnergy(const G4ParticleDefinition*,
    126                                                      const G4Material*,
    127                                                            G4double)
    128 {
    129   return lowestKinEnergy;
    130 }
    131 
    132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    133116
    134117#endif

  • trunk/source/processes/electromagnetic/muons/include/G4MuBremsstrahlungModel.hh

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuBremsstrahlungModel.hh,v 1.17 2007/10/11 09:25:31 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuBremsstrahlungModel.hh,v 1.22 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4545// 10-02-04 Add lowestKinEnergy (V.Ivanchenko)
    4646// 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
    47 // 13-02-06 add ComputeCrossSectionPerAtom (mma)
     47// 13-02-06 Add ComputeCrossSectionPerAtom (mma)
    4848// 11-10-07 Add ignoreCut flag (V.Ivanchenko)
     49// 28-02-08 Reorganized protected methods and members (V.Ivanchenko)
     50// 06-03-08 Remove obsolete methods and members (V.Ivanchenko)
    4951//
    5052
    … …  
    5254// Class Description:
    5355//
    54 // Implementation of energy loss for gamma emission by muons
     56// Implementation of bremssrahlung by muons
    5557
    5658// -------------------------------------------------------------------
    … …  
    6163
    6264#include "G4VEmModel.hh"
     65#include "G4NistManager.hh"
    6366
    6467class G4Element;
    … …  
    7578  virtual ~G4MuBremsstrahlungModel();
    7679
    77   void SetParticle(const G4ParticleDefinition*);
     80  virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
    7881
    79   void Initialise(const G4ParticleDefinition*, const G4DataVector&);
    80 
    81   void SetLowestKineticEnergy(G4double e) {lowestKinEnergy = e;};
    82 
    83   G4double MinEnergyCut(const G4ParticleDefinition*,
    84                         const G4MaterialCutsCouple*);
     82  virtual G4double MinEnergyCut(const G4ParticleDefinition*,
     83                                const G4MaterialCutsCouple*);
    8584                             
    8685  virtual G4double ComputeCrossSectionPerAtom(
    … …  
    9089                                 G4double cutEnergy,
    9190                                 G4double maxEnergy);
    92                                  
    93   virtual G4double CrossSectionPerVolume(const G4Material*,
    94                          const G4ParticleDefinition*,
    95                                G4double kineticEnergy,
    96                                G4double cutEnergy,
    97                                G4double maxEnergy);
    98                                
     91                                                               
    9992  virtual G4double ComputeDEDXPerVolume(const G4Material*,
    10093                                const G4ParticleDefinition*,
    … …  
    10295                                G4double cutEnergy);
    10396                             
    104   void SampleSecondaries(std::vector<G4DynamicParticle*>*,
    105                          const G4MaterialCutsCouple*,
    106                          const G4DynamicParticle*,
    107                          G4double tmin,
    108                          G4double maxEnergy);
     97  virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
     98                                 const G4MaterialCutsCouple*,
     99                                 const G4DynamicParticle*,
     100                                 G4double tmin,
     101                                 G4double maxEnergy);
     102
     103  inline void SetLowestKineticEnergy(G4double e);
    109104
    110105protected:
    111106
    112   G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
    113                               G4double kineticEnergy);
     107  G4double ComputMuBremLoss(G4double Z, G4double tkin, G4double cut);
     108 
     109  G4double ComputeMicroscopicCrossSection(G4double tkin,
     110                                          G4double Z,
     111                                          G4double cut);
    114112
    115 public:
     113  virtual G4double ComputeDMicroscopicCrossSection(G4double tkin,
     114                                                   G4double Z,
     115                                                   G4double gammaEnergy);
    116116
    117  G4double ComputMuBremLoss(G4double Z, G4double A, G4double tkin, G4double cut);
    118 
    119  G4double ComputeMicroscopicCrossSection(G4double tkin,
    120                                            G4double Z,
    121                                            G4double A,
    122                                            G4double cut);
    123 
    124  G4double ComputeDMicroscopicCrossSection(G4double tkin,
    125                                           G4double Z,
    126                                           G4double A,
    127                                           G4double gammaEnergy);
    128 
    129   inline void SetIgnoreCutFlag(G4bool);
    130 
    131   inline G4bool IgnoreCutFlag() const;
     117  inline void SetParticle(const G4ParticleDefinition*);
    132118
    133119private:
    134120
    135  G4DataVector* ComputePartialSumSigma(const G4Material* material,
    136                                       G4double tkin, G4double cut);
     121  G4DataVector* ComputePartialSumSigma(const G4Material* material,
     122                                       G4double tkin, G4double cut);
    137123
    138  const G4Element* SelectRandomAtom(const G4MaterialCutsCouple* couple) const;
    139 
    140  void MakeSamplingTables();
    141 
     124  const G4Element* SelectRandomAtom(const G4MaterialCutsCouple* couple) const;
    142125
    143126  // hide assignment operator
    … …  
    145128  G4MuBremsstrahlungModel(const  G4MuBremsstrahlungModel&);
    146129
     130protected:
     131
     132  const G4ParticleDefinition* particle;
     133  G4NistManager* nist;
     134  G4double mass;
     135  G4double rmass;
     136  G4double cc;
     137  G4double coeff;
     138  G4double sqrte;
     139  G4double bh;
     140  G4double bh1;
     141  G4double btf;
     142  G4double btf1;
     143
     144private:
     145
    147146  G4ParticleDefinition*       theGamma;
    148   const G4ParticleDefinition* particle;
    149147  G4ParticleChangeForLoss*    fParticleChange;
    150148
    … …  
    153151  G4double lowestKinEnergy;
    154152  G4double minThreshold;
    155   G4double mass;
    156 
    157   // tables for sampling
    158   G4int nzdat,ntdat,NBIN;
    159   static G4double zdat[5],adat[5],tdat[8];
    160   G4double ya[1001], proba[5][8][1001];
    161   G4double cutFixed;
    162 
    163   G4bool  ignoreCut;
    164153
    165154  std::vector<G4DataVector*> partialSumSigma;
    166   G4bool  samplingTablesAreFilled;
    167 
    168155};
    169156
    170157//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    171158
    172 inline G4double G4MuBremsstrahlungModel::MaxSecondaryEnergy(
    173                                  const G4ParticleDefinition*,
    174                                  G4double kineticEnergy)
     159inline void G4MuBremsstrahlungModel::SetLowestKineticEnergy(G4double e)
    175160{
    176   return kineticEnergy;
     161  lowestKinEnergy = e;
    177162}
    178163
    179164//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    180165
    181 inline void G4MuBremsstrahlungModel::SetIgnoreCutFlag(G4bool val)
     166inline
     167void G4MuBremsstrahlungModel::SetParticle(const G4ParticleDefinition* p)
    182168{
    183   ignoreCut = val;
    184 }
    185 
    186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    187 
    188 inline G4bool G4MuBremsstrahlungModel::IgnoreCutFlag() const
    189 {
    190   return ignoreCut;
     169  if(!particle) {
     170    particle = p;
     171    mass = particle->GetPDGMass();
     172    rmass=mass/electron_mass_c2 ;
     173    cc=classic_electr_radius/rmass ;
     174    coeff= 16.*fine_structure_const*cc*cc/3. ;
     175  }
    191176}
    192177

  • trunk/source/processes/electromagnetic/muons/include/G4MuIonisation.hh

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuIonisation.hh,v 1.30 2007/05/23 08:49:32 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: G4MuIonisation.hh,v 1.31 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9595  virtual ~G4MuIonisation();
    9696
    97   G4bool IsApplicable(const G4ParticleDefinition& p);
     97  virtual G4bool IsApplicable(const G4ParticleDefinition& p);
    9898
    99   G4double MinPrimaryEnergy(const G4ParticleDefinition* p,
    100                             const G4Material*, G4double cut);
     99  virtual G4double MinPrimaryEnergy(const G4ParticleDefinition* p,
     100                                    const G4Material*, G4double cut);
    101101
    102102  // Print out of the class parameters
    103   void PrintInfo();
     103  virtual void PrintInfo();
    104104
    105105protected:
    … …  
    127127
    128128//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    130 
    131 inline G4bool G4MuIonisation::IsApplicable(const G4ParticleDefinition& p)
    132 {
    133   return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV);
    134 }
    135 
    136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    137 
    138 inline G4double G4MuIonisation::MinPrimaryEnergy(const G4ParticleDefinition*,
    139                                                  const G4Material*,
    140                                                  G4double cut)
    141 {
    142   G4double x = 0.5*cut/electron_mass_c2;
    143   G4double g = x*ratio + std::sqrt((1. + x)*(1. + x*ratio*ratio));
    144   return mass*(g - 1.0);
    145 }
    146 
    147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    148129
    149130#endif

  • trunk/source/processes/electromagnetic/muons/include/G4MuMultipleScattering.hh

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuMultipleScattering.hh,v 1.1 2007/10/26 09:52:37 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuMultipleScattering.hh,v 1.2 2008/04/13 17:19:13 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -----------------------------------------------------------------------------
    … …  
    6060//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    6161
    62 class G4MuMscModel;
     62class G4VMscModel;
    6363
    6464class G4MuMultipleScattering : public G4VMultipleScattering
    … …  
    9494private:        // data members
    9595
    96   G4MuMscModel* mscModel;
     96  G4VMscModel* mscModel;
    9797
    9898  G4double thetaLimit;

  • trunk/source/processes/electromagnetic/muons/include/G4MuPairProduction.hh

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuPairProduction.hh,v 1.29 2007/05/23 08:49:32 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuPairProduction.hh,v 1.31 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8484  virtual ~G4MuPairProduction();
    8585
    86   G4bool IsApplicable(const G4ParticleDefinition& p);
     86  virtual G4bool IsApplicable(const G4ParticleDefinition& p);
    8787
    88   G4double MinPrimaryEnergy(const G4ParticleDefinition* p,
    89                             const G4Material*, G4double cut);
     88  virtual G4double MinPrimaryEnergy(const G4ParticleDefinition* p,
     89                                    const G4Material*, G4double cut);
    9090
    9191  // Print out of the class parameters
    92   void PrintInfo();
     92  virtual void PrintInfo();
    9393
    9494protected:
    9595
    96   void InitialiseEnergyLossProcess(const G4ParticleDefinition*,
    97                                    const G4ParticleDefinition*);
     96  virtual void InitialiseEnergyLossProcess(const G4ParticleDefinition*,
     97                                           const G4ParticleDefinition*);
    9898
    9999private:
    … …  
    113113
    114114//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    116 
    117 inline G4bool G4MuPairProduction::IsApplicable(const G4ParticleDefinition& p)
    118 {
    119   return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV);
    120 }
    121 
    122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    123 
    124 inline G4double G4MuPairProduction::MinPrimaryEnergy(const G4ParticleDefinition*,
    125                                                      const G4Material*,
    126                                                            G4double)
    127 {
    128   return lowestKinEnergy;
    129 }
    130 
    131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    132115
    133116#endif

  • trunk/source/processes/electromagnetic/muons/include/G4MuPairProductionModel.hh

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuPairProductionModel.hh,v 1.24 2007/10/11 13:52:03 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuPairProductionModel.hh,v 1.28 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4848// 12-05-06 Add parameter to SelectRandomAtom (A.Bogdanov)
    4949// 11-10-07 Add ignoreCut flag (V.Ivanchenko)
     50// 28-02-08 Reorganized protected methods and members (V.Ivanchenko)
    5051
    5152//
    … …  
    6263
    6364#include "G4VEmModel.hh"
     65#include "G4NistManager.hh"
    6466#include <vector>
    6567
    … …  
    7072class G4MuPairProductionModel : public G4VEmModel
    7173{
    72 
    7374public:
    7475
    … …  
    7879  virtual ~G4MuPairProductionModel();
    7980
    80   void SetParticle(const G4ParticleDefinition*);
    81 
    82   void Initialise(const G4ParticleDefinition*, const G4DataVector&);
    83 
    84   void SetLowestKineticEnergy(G4double e) {lowestKinEnergy = e;};
    85 
    86   G4double MinEnergyCut(const G4ParticleDefinition*,
    87                         const G4MaterialCutsCouple*);
     81  virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
     82
    8883                       
    8984  virtual G4double ComputeCrossSectionPerAtom(
    … …  
    9489                                 G4double maxEnergy);
    9590                                 
    96   virtual G4double CrossSectionPerVolume(const G4Material*,
    97                          const G4ParticleDefinition*,
    98                                G4double kineticEnergy,
    99                                G4double cutEnergy,
    100                                G4double maxEnergy);
    101                                
    10291  virtual G4double ComputeDEDXPerVolume(const G4Material*,
    10392                                const G4ParticleDefinition*,
    … …  
    10594                                G4double cutEnergy);
    10695
    107   void SampleSecondaries(std::vector<G4DynamicParticle*>*,
    108                          const G4MaterialCutsCouple*,
    109                          const G4DynamicParticle*,
    110                          G4double tmin,
    111                          G4double maxEnergy);
     96  virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
     97                                 const G4MaterialCutsCouple*,
     98                                 const G4DynamicParticle*,
     99                                 G4double tmin,
     100                                 G4double maxEnergy);
     101
     102  virtual G4double MinEnergyCut(const G4ParticleDefinition*,
     103                               const G4MaterialCutsCouple*);
     104
     105  inline void SetLowestKineticEnergy(G4double e);
     106
     107  inline void SetParticle(const G4ParticleDefinition*);
    112108
    113109protected:
    114 
    115   inline G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
    116                                      G4double kineticEnergy);
    117 
    118 
    119 public:
    120110
    121111  G4double ComputMuPairLoss(G4double Z, G4double tkin, G4double cut,
    … …  
    126116                                          G4double cut);
    127117
    128   G4double ComputeDMicroscopicCrossSection(G4double tkin,
    129                                            G4double Z,
    130                                            G4double pairEnergy);
    131 
    132   inline void SetIgnoreCutFlag(G4bool);
    133 
    134   inline G4bool IgnoreCutFlag() const;
     118  virtual G4double ComputeDMicroscopicCrossSection(G4double tkin,
     119                                                   G4double Z,
     120                                                   G4double pairEnergy);
     121
     122  virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
     123                                      G4double kineticEnergy);
     124
     125  inline void SetCurrentElement(G4double Z);
    135126
    136127private:
    … …  
    144135  void MakeSamplingTables();
    145136
    146   void SetCurrentElement(G4double Z);
    147 
    148137  inline G4double InterpolatedIntegralCrossSection(
    149138                     G4double dt, G4double dz, G4int iz,
    … …  
    154143  G4MuPairProductionModel(const  G4MuPairProductionModel&);
    155144
    156   G4ParticleDefinition*       theElectron;
    157   G4ParticleDefinition*       thePositron;
    158   G4ParticleChangeForLoss*    fParticleChange;
    159   G4ParticleChangeForGamma*   gParticleChange;
    160 
    161   G4double minPairEnergy;
    162   G4double lowestKinEnergy;
     145protected:
     146
     147  const G4ParticleDefinition* particle;
     148  G4NistManager*              nist;
    163149
    164150  G4double factorForCross;
    … …  
    170156  G4double lnZ;
    171157
    172   const G4ParticleDefinition* particle;
     158  static G4double xgi[8],wgi[8];
     159
     160private:
     161
     162  G4ParticleDefinition*       theElectron;
     163  G4ParticleDefinition*       thePositron;
     164  G4ParticleChangeForLoss*    fParticleChange;
     165
     166  G4double minPairEnergy;
     167  G4double lowestKinEnergy;
    173168
    174169  // tables for sampling
    … …  
    177172  G4int nbiny;
    178173  size_t nmaxElements;
    179   static G4double zdat[5],adat[5],tdat[8],xgi[8],wgi[8];
     174  static G4double zdat[5],adat[5],tdat[8];
    180175  G4double ya[1001],proba[5][8][1001];
    181176
    … …  
    184179  G4double dy;
    185180
    186   G4bool  ignoreCut;
    187 
    188181  G4bool  samplingTablesAreFilled;
    189182  std::vector<G4double>   partialSum;
    … …  
    192185//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    193186
    194 inline G4double G4MuPairProductionModel::MaxSecondaryEnergy(
    195                                  const G4ParticleDefinition*,
    196                                        G4double kineticEnergy)
    197 {
    198   G4double maxPairEnergy = kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
    199   return maxPairEnergy;
     187inline void G4MuPairProductionModel::SetLowestKineticEnergy(G4double e)
     188{
     189  lowestKinEnergy = e;
     190}
     191
     192//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     193
     194inline
     195void G4MuPairProductionModel::SetParticle(const G4ParticleDefinition* p)
     196{
     197  if(!particle) {
     198    particle = p;
     199    particleMass = particle->GetPDGMass();
     200  }
    200201}
    201202
    … …  
    206207  if(Z != currentZ) {
    207208    currentZ = Z;
    208     z13 = std::pow(Z,0.333333333);
     209    G4int iz = G4int(Z);
     210    z13 = nist->GetZ13(iz);
    209211    z23 = z13*z13;
    210     lnZ = std::log(Z);
     212    lnZ = nist->GetLOGZ(iz);
    211213  }
    212214}
    … …  
    229231//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    230232
    231 inline void G4MuPairProductionModel::SetIgnoreCutFlag(G4bool val)
    232 {
    233   ignoreCut = val;
    234 }
    235 
    236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    237 
    238 inline G4bool G4MuPairProductionModel::IgnoreCutFlag() const
    239 {
    240   return ignoreCut;
    241 }
    242 
    243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    244 
    245233#endif

  • trunk/source/processes/electromagnetic/muons/src/G4EnergyLossForExtrapolator.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EnergyLossForExtrapolator.cc,v 1.13 2007/07/28 13:44:25 vnivanch Exp $
    27 // GEANT4 tag $Name: $
     26// $Id: G4EnergyLossForExtrapolator.cc,v 1.18 2008/11/13 14:14:07 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929//---------------------------------------------------------------------------
    … …  
    6767#include "G4MuBremsstrahlungModel.hh"
    6868#include "G4ProductionCuts.hh"
     69#include "G4LossTableManager.hh"
     70#include "G4WentzelVIModel.hh"
    6971
    7072//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    8890  delete invRangePositron;
    8991  delete invRangeProton;
     92  delete mscElectron;
    9093  delete cuts;
    9194}
    … …  
    100103  if(!isInitialised) Initialisation();
    101104  G4double kinEnergyFinal = kinEnergy;
    102   if(mat && part) {
    103     G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
     105  if(SetupKinematics(part, mat, kinEnergy)) {
     106    G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
    104107    G4double r  = ComputeRange(kinEnergy,part);
    105108    if(r <= step) {
    … …  
    125128  G4double kinEnergyFinal = kinEnergy;
    126129
    127   if(mat && part) {
    128     G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
     130  if(SetupKinematics(part, mat, kinEnergy)) {
     131    G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
    129132    G4double r  = ComputeRange(kinEnergy,part);
    130133
    … …  
    141144//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    142145
    143 G4double G4EnergyLossForExtrapolator::ComputeTrueStep(const G4Material* mat,
    144                                                       const G4ParticleDefinition* part,
    145                                                       G4double kinEnergy, G4double stepLength)
    146 {
     146G4double G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy,
     147                                                     G4double stepLength,
     148                                                     const G4Material* mat,
     149                                                     const G4ParticleDefinition* part)
     150{
     151  G4double res = stepLength;
     152  if(!isInitialised) Initialisation();
     153  if(SetupKinematics(part, mat, kinEnergy)) {
     154    if(part == electron || part == positron) {
     155      G4double x = stepLength*ComputeValue(kinEnergy, mscElectron);
     156      if(x < 0.2) res *= (1.0 + 0.5*x + x*x/3.0);
     157      else if(x < 0.9999) res = -std::log(1.0 - x)*stepLength/x;
     158      else res = ComputeRange(kinEnergy,part);
     159   
     160    } else {
     161      res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
     162    }
     163  }
     164  return res;
     165}
     166
     167//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     168
     169G4bool G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part,
     170                                                    const G4Material* mat,
     171                                                    G4double kinEnergy)
     172{
     173  if(!part || !mat || kinEnergy < keV) return false;
    147174  if(!isInitialised) Initialisation();
    148175  G4bool flag = false;
    … …  
    182209    if(tmax > maxEnergyTransfer) tmax = maxEnergyTransfer;
    183210  }
    184   G4double theta = ComputeScatteringAngle(stepLength);
    185   return stepLength*std::sqrt(1.0 + 0.625*theta*theta);
     211  return true;
    186212}
    187213
    … …  
    231257  invRangeMuon     = PrepareTable();
    232258  invRangeProton   = PrepareTable();
     259  mscElectron      = PrepareTable();
    233260
    234261  G4LossTableBuilder builder;
    … …  
    262289  builder.BuildInverseRangeTable(rangeProton, invRangeProton); 
    263290
     291  ComputeTrasportXS(electron, mscElectron);
    264292}
    265293
    … …  
    273301
    274302    G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins);
     303    v->SetSpline(G4LossTableManager::Instance()->SplineFlag());
    275304    table->push_back(v);
    276305  }
    … …  
    488517//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    489518
     519void G4EnergyLossForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part,
     520                                                    G4PhysicsTable* table)
     521{
     522  G4DataVector v;
     523  G4WentzelVIModel* msc = new G4WentzelVIModel();
     524  msc->SetPolarAngleLimit(CLHEP::pi);
     525  msc->Initialise(part, v);
     526
     527  mass    = part->GetPDGMass();
     528  charge2 = 1.0;
     529  currentParticle = part;
     530
     531  const G4MaterialTable* mtable = G4Material::GetMaterialTable();
     532
     533  if(0<verbose) {
     534    G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName()
     535           << G4endl;
     536  }
     537 
     538  for(G4int i=0; i<nmat; i++) { 
     539
     540    const G4Material* mat = (*mtable)[i];
     541    if(1<verbose)
     542      G4cout << "i= " << i << "  mat= " << mat->GetName() << G4endl;
     543    G4PhysicsVector* aVector = (*table)[i];
     544    for(G4int j=0; j<nbins; j++) {
     545       
     546       G4double e = aVector->GetLowEdgeEnergy(j);
     547       G4double xs = msc->CrossSectionPerVolume(mat,part,e);
     548       aVector->PutValue(j,xs);
     549       if(1<verbose) {
     550         G4cout << "j= " << j << "  e(MeV)= " << e/MeV 
     551                << " xs(1/mm)= " << xs*mm << G4endl;
     552       }
     553    }
     554  }
     555  delete msc;
     556}
     557
     558//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     559

  • trunk/source/processes/electromagnetic/muons/src/G4MuBetheBlochModel.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuBetheBlochModel.cc,v 1.23 2007/05/22 17:35:58 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuBetheBlochModel.cc,v 1.25 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8888  theElectron = G4Electron::Electron();
    8989  corr = G4LossTableManager::Instance()->EmCorrections();
     90  fParticleChange = 0;
    9091
    9192  if(p) SetParticle(p);
    … …  
    99100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    100101
    101 void G4MuBetheBlochModel::SetParticle(const G4ParticleDefinition* p)
    102 {
    103   if(!particle) {
    104     particle = p;
    105     mass = particle->GetPDGMass();
    106     massSquare = mass*mass;
    107     ratio = electron_mass_c2/mass;
    108   }
    109 }
    110 
    111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    112 
    113102G4double G4MuBetheBlochModel::MinEnergyCut(const G4ParticleDefinition*,
    114103                                           const G4MaterialCutsCouple* couple)
    … …  
    117106}
    118107
     108//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     109
     110G4double G4MuBetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
     111                                                 G4double kinEnergy)
     112{
     113  G4double tau  = kinEnergy/mass;
     114  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
     115                  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
     116  return tmax;
     117}
     118
    119119//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    120120
    … …  
    124124  if(p) SetParticle(p);
    125125
    126   if(pParticleChange)
    127     fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>
    128                                                              (pParticleChange);
    129   else
    130     fParticleChange = new G4ParticleChangeForLoss();
     126  if(!fParticleChange) {
     127    if(pParticleChange) {
     128      fParticleChange =
     129        reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
     130    } else {
     131      fParticleChange = new G4ParticleChangeForLoss();
     132    }
     133  }
    131134}
    132135
    … …  
    275278
    276279  //High order corrections
    277   dedx += corr->HighOrderCorrections(p,material,kineticEnergy);
     280  dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
    278281
    279282  return dedx;

  • trunk/source/processes/electromagnetic/muons/src/G4MuBremsstrahlung.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuBremsstrahlung.cc,v 1.38 2007/05/22 17:35:58 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuBremsstrahlung.cc,v 1.42 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8282    lowestKinEnergy(1.*GeV),
    8383    isInitialised(false)
    84 {}
     84{
     85  SetProcessSubType(fBremsstrahlung);
     86}
    8587
    8688//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    9193//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    9294
    93 void G4MuBremsstrahlung::InitialiseEnergyLossProcess(const G4ParticleDefinition* part,
    94                                                      const G4ParticleDefinition*)
     95G4bool G4MuBremsstrahlung::IsApplicable(const G4ParticleDefinition& p)
     96{
     97  return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV);
     98}
     99
     100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     101
     102G4double G4MuBremsstrahlung::MinPrimaryEnergy(const G4ParticleDefinition*,
     103                                              const G4Material*,
     104                                              G4double)
     105{
     106  return lowestKinEnergy;
     107}
     108
     109//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     110
     111void G4MuBremsstrahlung::InitialiseEnergyLossProcess(
     112                                 const G4ParticleDefinition* part,
     113                                 const G4ParticleDefinition*)
    95114{
    96115  if(!isInitialised) {
    … …  
    105124    em->SetLowestKineticEnergy(lowestKinEnergy);
    106125
    107     G4VEmFluctuationModel* fm = new G4UniversalFluctuation();
    108     em->SetLowEnergyLimit(0.1*keV);
    109     em->SetHighEnergyLimit(100.0*TeV);
     126    G4VEmFluctuationModel* fm = 0;
     127    em->SetLowEnergyLimit(MinKinEnergy());
     128    em->SetHighEnergyLimit(MaxKinEnergy());
    110129    AddEmModel(1, em, fm);
    111130  }
    … …  
    115134
    116135void G4MuBremsstrahlung::PrintInfo()
    117 {
    118   G4cout << "      Parametrised model "
    119          << G4endl;
    120 }
     136{}
    121137
    122138//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

  • trunk/source/processes/electromagnetic/muons/src/G4MuBremsstrahlungModel.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuBremsstrahlungModel.cc,v 1.24 2007/11/08 11:48:28 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuBremsstrahlungModel.cc,v 1.33 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4646// 13-02-03 Add name (V.Ivanchenko)
    4747// 10-02-04 Add lowestKinEnergy (V.Ivanchenko)
    48 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
    49 // 03-08-05 Angular correlations according to PRM (V.Ivantchenko)
     48// 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
     49// 03-08-05 Angular correlations according to PRM (V.Ivanchenko)
    5050// 13-02-06 add ComputeCrossSectionPerAtom (mma)
    5151// 21-03-06 Fix problem of initialisation in case when cuts are not defined (VI)
    5252// 07-11-07 Improve sampling of final state (A.Bogdanov)
     53// 28-02-08 Use precomputed Z^1/3 and Log(A) (V.Ivanchenko)
    5354//
    5455
    … …  
    7475
    7576//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    76 
    77 // static members
    78 //
    79 G4double G4MuBremsstrahlungModel::zdat[]={1., 4., 13., 29., 92.};
    80 G4double G4MuBremsstrahlungModel::adat[]={1.01, 9.01, 26.98, 63.55, 238.03};
    81 G4double G4MuBremsstrahlungModel::tdat[]={1.e3, 1.e4, 1.e5, 1.e6, 1.e7, 1.e8,
    82                                           1.e9, 1.e10};
    83 
    8477//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    8578
    … …  
    9083  : G4VEmModel(nam),
    9184    particle(0),
     85    sqrte(sqrt(exp(1.))),
     86    bh(202.4),
     87    bh1(446.),
     88    btf(183.),
     89    btf1(1429.),
     90    fParticleChange(0),
    9291    lowestKinEnergy(1.0*GeV),
    93     minThreshold(1.0*keV),
    94     nzdat(5),
    95     ntdat(8),
    96     NBIN(1000),
    97     cutFixed(0.98*keV),
    98     ignoreCut(false),
    99     samplingTablesAreFilled(false)
     92    minThreshold(1.0*keV)
    10093{
    10194  theGamma = G4Gamma::Gamma();
     95  nist = G4NistManager::Instance();
    10296  if(p) SetParticle(p);
    10397}
    … …  
    125119//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    126120
    127 void G4MuBremsstrahlungModel::SetParticle(const G4ParticleDefinition* p)
    128 {
    129   if(!particle) {
    130     particle = p;
    131     mass = particle->GetPDGMass();
    132   }
    133 }
    134 
    135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    136 
    137121void G4MuBremsstrahlungModel::Initialise(const G4ParticleDefinition* p,
    138122                                         const G4DataVector& cuts)
    … …  
    142126  highKinEnergy = HighEnergyLimit();
    143127
     128  // partial cross section is computed for fixed energy
    144129  G4double fixedEnergy = 0.5*highKinEnergy;
    145130
    … …  
    148133  if(theCoupleTable) {
    149134    G4int numOfCouples = theCoupleTable->GetTableSize();
    150    
     135
     136    // clear old data   
    151137    G4int nn = partialSumSigma.size();
    152138    G4int nc = cuts.size();
    153139    if(nn > 0) {
    154140      for (G4int ii=0; ii<nn; ii++){
    155         G4DataVector* a=partialSumSigma[ii];
     141        G4DataVector* a = partialSumSigma[ii];
    156142        if ( a )  delete a;   
    157143      }
    158144      partialSumSigma.clear();
    159145    }
     146    // fill new data
    160147    if (numOfCouples>0) {
    161148      for (G4int i=0; i<numOfCouples; i++) {
    162149        G4double cute = DBL_MAX;
     150
     151        // protection for usage with extrapolator
    163152        if(i < nc) cute = cuts[i];
    164         if(cute < cutFixed || ignoreCut) cute = cutFixed;
     153
    165154        const G4MaterialCutsCouple* couple =
    166155          theCoupleTable->GetMaterialCutsCouple(i);
    … …  
    171160    }
    172161  }
    173   if(!samplingTablesAreFilled) MakeSamplingTables();
    174   if(pParticleChange)
    175     fParticleChange =
    176       reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
    177   else
    178     fParticleChange = new G4ParticleChangeForLoss();
     162
     163  // define pointer to G4ParticleChange
     164  if(!fParticleChange) {
     165    if(pParticleChange)
     166      fParticleChange =
     167        reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
     168    else
     169      fParticleChange = new G4ParticleChangeForLoss();
     170  }
    179171}
    180172
    … …  
    188180{
    189181  G4double dedx = 0.0;
    190   if (kineticEnergy <= lowestKinEnergy || ignoreCut) return dedx;
     182  if (kineticEnergy <= lowestKinEnergy) return dedx;
    191183
    192184  G4double tmax = kineticEnergy;
    193   G4double cut  = min(cutEnergy,tmax);
    194   if(cut < cutFixed) cut = cutFixed;
     185  G4double cut  = std::min(cutEnergy,tmax);
     186  if(cut < minThreshold) cut = minThreshold;
    195187
    196188  const G4ElementVector* theElementVector = material->GetElementVector();
    197189  const G4double* theAtomicNumDensityVector =
    198                                           material->GetAtomicNumDensityVector();
     190    material->GetAtomicNumDensityVector();
    199191
    200192  //  loop for elements in the material
    201193  for (size_t i=0; i<material->GetNumberOfElements(); i++) {
    202194
    203     G4double Z = (*theElementVector)[i]->GetZ();
    204     G4double A = (*theElementVector)[i]->GetA()/(g/mole) ;
    205 
    206     G4double loss = ComputMuBremLoss(Z, A, kineticEnergy, cut);
     195    G4double loss =
     196      ComputMuBremLoss((*theElementVector)[i]->GetZ(), kineticEnergy, cut);
    207197
    208198    dedx += loss*theAtomicNumDensityVector[i];
    209199  }
     200  //  G4cout << "BR e= " << kineticEnergy << "  dedx= " << dedx << G4endl;
    210201  if(dedx < 0.) dedx = 0.;
    211202  return dedx;
    … …  
    214205//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    215206
    216 G4double G4MuBremsstrahlungModel::ComputMuBremLoss(G4double Z, G4double A,
     207G4double G4MuBremsstrahlungModel::ComputMuBremLoss(G4double Z,
    217208                                                   G4double tkin, G4double cut)
    218209{
    … …  
    239230    {
    240231      G4double ep = (aa + xgi[i]*hhh)*totalEnergy;
    241       loss += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, A, ep);
     232      loss += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
    242233    }
    243234    aa += hhh;
    … …  
    254245                                           G4double tkin,
    255246                                           G4double Z,
    256                                            G4double A,
    257247                                           G4double cut)
    258248{
    … …  
    272262  G4double bbb = log(vmax);
    273263  G4int    kkk = (G4int)((bbb-aaa)/ak1)+k2 ;
    274   G4double hhh = (bbb-aaa)/float(kkk);
     264  G4double hhh = (bbb-aaa)/G4double(kkk);
    275265
    276266  G4double aa = aaa;
    … …  
    281271    {
    282272      G4double ep = exp(aa + xgi[i]*hhh)*totalEnergy;
    283       cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, A, ep);
     273      cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
    284274    }
    285275    aa += hhh;
    … …  
    287277
    288278  cross *=hhh;
     279
     280  //G4cout << "BR e= " << tkin<< "  cross= " << cross/barn << G4endl;
    289281
    290282  return cross;
    … …  
    296288                                           G4double tkin,
    297289                                           G4double Z,
    298                                            G4double A,
    299290                                           G4double gammaEnergy)
    300291//  differential cross section
    301292{
    302   static const G4double sqrte=sqrt(exp(1.)) ;
    303   static const G4double bh=202.4,bh1=446.,btf=183.,btf1=1429. ;
    304   static const G4double rmass=mass/electron_mass_c2 ;
    305   static const G4double cc=classic_electr_radius/rmass ;
    306   static const G4double coeff= 16.*fine_structure_const*cc*cc/3. ;
    307 
    308293  G4double dxsection = 0.;
    309294
    … …  
    315300  G4double rab0=delta*sqrte ;
    316301
    317   G4double z13 = exp(-log(Z)/3.) ;
    318   G4double dn  = 1.54*exp(0.27*log(A)) ;
     302  G4int iz = G4int(Z);
     303  if(iz < 1) iz = 1;
     304
     305  G4double z13 = 1.0/nist->GetZ13(iz);
     306  G4double dn  = 1.54*nist->GetA27(iz);
    319307
    320308  G4double b,b1,dnstar ;
    321309
    322   if(Z<1.5)
     310  if(1 == iz)
    323311  {
    324     b=bh;
    325     b1=bh1;
    326     dnstar=dn ;
     312    b  = bh;
     313    b1 = bh1;
     314    dnstar = dn;
    327315  }
    328316  else
    329317  {
    330     b=btf;
    331     b1=btf1;
    332     dnstar = exp((1.-1./Z)*log(dn)) ;
     318    b  = btf;
     319    b1 = btf1;
     320    dnstar = dn/std::pow(dn, 1./Z);
    333321  }
    334322
    … …  
    359347                                           const G4ParticleDefinition*,
    360348                                                 G4double kineticEnergy,
    361                                                  G4double Z, G4double A,
     349                                                 G4double Z, G4double,
    362350                                                 G4double cutEnergy,
    363                                                  G4double)
    364 {
    365   G4double cut  = min(cutEnergy, kineticEnergy);
    366   if(cut < cutFixed || ignoreCut) cut = cutFixed;
    367   G4double cross =
    368     ComputeMicroscopicCrossSection (kineticEnergy, Z, A/(g/mole), cut);
    369   return cross;
    370 }
    371 
    372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    373 
    374 G4double G4MuBremsstrahlungModel::CrossSectionPerVolume(
    375                                                const G4Material* material,
    376                                                const G4ParticleDefinition*,
    377                                                      G4double kineticEnergy,
    378                                                      G4double cutEnergy,
    379                                                      G4double maxEnergy)
     351                                                 G4double maxEnergy)
    380352{
    381353  G4double cross = 0.0;
    382   if (cutEnergy >= maxEnergy || kineticEnergy <= lowestKinEnergy) return cross;
    383  
    384   G4double tmax = min(maxEnergy, kineticEnergy);
    385   G4double cut  = min(cutEnergy, tmax);
    386   if(cut < cutFixed || ignoreCut) cut = cutFixed;
    387 
    388   const G4ElementVector* theElementVector = material->GetElementVector();
    389   const G4double* theAtomNumDensityVector =
    390                                           material->GetAtomicNumDensityVector();
    391 
    392   for (size_t i=0; i<material->GetNumberOfElements(); i++) {
    393 
    394     G4double Z = (*theElementVector)[i]->GetZ();
    395     G4double A = (*theElementVector)[i]->GetA()/(g/mole);
    396 
    397     G4double cr = ComputeMicroscopicCrossSection(kineticEnergy, Z, A, cut);
    398 
    399     if(tmax < kineticEnergy) {
    400       cr -= ComputeMicroscopicCrossSection(kineticEnergy, Z, A, tmax);
    401     }
    402     cross += theAtomNumDensityVector[i] * cr;
    403   }
    404 
     354  if (kineticEnergy <= lowestKinEnergy) return cross;
     355  G4double tmax = std::min(maxEnergy, kineticEnergy);
     356  G4double cut  = std::min(cutEnergy, kineticEnergy);
     357  if(cut < minThreshold) cut = minThreshold;
     358  if (cut >= tmax) return cross;
     359
     360  cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
     361  if(tmax < kineticEnergy) {
     362    cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
     363  }
    405364  return cross;
    406365}
    … …  
    410369G4DataVector* G4MuBremsstrahlungModel::ComputePartialSumSigma(
    411370                                       const G4Material* material,
    412                                              G4double kineticEnergy,
    413                                              G4double cut)
    414 
    415 // Build the table of cross section per element. The table is built for MATERIAL
    416 // This table is used by DoIt to select randomly an element in the material.
     371                                       G4double kineticEnergy,
     372                                       G4double cut)
     373
     374// Build the table of cross section per element.
     375// The table is built for material
     376// This table is used to select randomly an element in the material.
    417377{
    418378  G4int nElements = material->GetNumberOfElements();
    419379  const G4ElementVector* theElementVector = material->GetElementVector();
    420380  const G4double* theAtomNumDensityVector =
    421                                           material->GetAtomicNumDensityVector();
     381    material->GetAtomicNumDensityVector();
    422382
    423383  G4DataVector* dv = new G4DataVector();
    … …  
    426386
    427387  for (G4int i=0; i<nElements; i++ ) {
    428 
    429     G4double Z = (*theElementVector)[i]->GetZ();
    430     G4double A = (*theElementVector)[i]->GetA()/(g/mole) ;
    431388    cross += theAtomNumDensityVector[i]
    432              * ComputeMicroscopicCrossSection(kineticEnergy, Z, A, cut);
     389      * ComputeMicroscopicCrossSection(kineticEnergy,
     390                                       (*theElementVector)[i]->GetZ(), cut);
    433391    dv->push_back(cross);
    434392  }
    … …  
    438396//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    439397
    440 void G4MuBremsstrahlungModel::MakeSamplingTables()
    441 {
    442 
    443   G4double AtomicNumber,AtomicWeight,KineticEnergy,
    444            TotalEnergy,Maxep;
    445 
    446   for (G4int iz=0; iz<nzdat; iz++)
    447    {
    448      AtomicNumber = zdat[iz];
    449      AtomicWeight = adat[iz]*g/mole ;
    450 
    451      for (G4int it=0; it<ntdat; it++)
    452      {
    453        KineticEnergy = tdat[it];
    454        TotalEnergy = KineticEnergy + mass;
    455        Maxep = KineticEnergy ;
    456 
    457        G4double CrossSection = 0.0 ;
    458 
    459        // calculate the differential cross section
    460        // numerical integration in
    461        //  log ...............
    462        G4double c = log(Maxep/cutFixed) ;
    463        G4double ymin = -5. ;
    464        G4double ymax = 0. ;
    465        G4double dy = (ymax-ymin)/NBIN ;
    466 
    467        G4double y = ymin - 0.5*dy ;
    468        G4double yy = ymin - dy ;
    469        G4double x = exp(y);
    470        G4double fac = exp(dy);
    471        G4double dx = exp(yy)*(fac - 1.0);
    472 
    473        for (G4int i=0 ; i<NBIN; i++)
    474        {
    475          y += dy ;
    476          x *= fac;
    477          dx*= fac;
    478          G4double ep = cutFixed*exp(c*x) ;
    479 
    480          CrossSection += ep*dx*ComputeDMicroscopicCrossSection(
    481                                                  KineticEnergy,AtomicNumber,
    482                                                  AtomicWeight,ep) ;
    483          ya[i]=y ;
    484          proba[iz][it][i] = CrossSection ;
    485 
    486        }
    487 
    488        proba[iz][it][NBIN] = CrossSection ;
    489        ya[NBIN] = 0. ;   //   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    490 
    491        if(CrossSection > 0.)
    492        {
    493          for(G4int ib=0; ib<=NBIN; ib++)
    494          {
    495            proba[iz][it][ib] /= CrossSection ;
    496          }
    497        }
    498      }
    499    }
    500   samplingTablesAreFilled = true;
    501 }
    502 
    503 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    504 
    505 void G4MuBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
    506                                                 const G4MaterialCutsCouple* couple,
    507                                                 const G4DynamicParticle* dp,
    508                                                 G4double minEnergy,
    509                                                 G4double maxEnergy)
     398void G4MuBremsstrahlungModel::SampleSecondaries(
     399                              std::vector<G4DynamicParticle*>* vdp,
     400                              const G4MaterialCutsCouple* couple,
     401                              const G4DynamicParticle* dp,
     402                              G4double minEnergy,
     403                              G4double maxEnergy)
    510404{
    511405  G4double kineticEnergy = dp->GetKineticEnergy();
    512406  // check against insufficient energy
    513   G4double tmax = min(kineticEnergy, maxEnergy);
    514   G4double tmin = min(kineticEnergy, minEnergy);
    515   if(tmin < cutFixed || ignoreCut) tmin = cutFixed;
     407  G4double tmax = std::min(kineticEnergy, maxEnergy);
     408  G4double tmin = std::min(kineticEnergy, minEnergy);
     409  if(tmin < minThreshold) tmin = minThreshold;
    516410  if(tmin >= tmax) return;
    517411
    518   // ===== the begining of a new code  ======
    519412  // ===== sampling of energy transfer ======
    520413
    … …  
    523416  // select randomly one element constituing the material
    524417  const G4Element* anElement = SelectRandomAtom(couple);
     418  G4double Z = anElement->GetZ();
    525419
    526420  G4double totalEnergy   = kineticEnergy + mass;
    527421  G4double totalMomentum = sqrt(kineticEnergy*(kineticEnergy + 2.0*mass));
    528422
    529   G4double AtomicNumber = anElement->GetZ();
    530   G4double AtomicWeight = anElement->GetA()/(g/mole);
    531 
    532   G4double func1 = tmin*ComputeDMicroscopicCrossSection(
    533                                     kineticEnergy,AtomicNumber,
    534                                     AtomicWeight,tmin);
     423  G4double func1 = tmin*
     424    ComputeDMicroscopicCrossSection(kineticEnergy,Z,tmin);
    535425
    536426  G4double lnepksi, epksi;
    537427  G4double func2;
    538   G4double ksi2;
    539428
    540429  do {
    541430    lnepksi = log(tmin) + G4UniformRand()*log(kineticEnergy/tmin);
    542431    epksi   = exp(lnepksi);
    543     func2   = epksi*ComputeDMicroscopicCrossSection(
    544                                 kineticEnergy,AtomicNumber,
    545                                 AtomicWeight,epksi);
    546     ksi2 = G4UniformRand();
    547 
    548   } while(func2/func1 < ksi2);
    549 
    550   // ===== the end of a new code =====
    551 
    552   // create G4DynamicParticle object for the Gamma
     432    func2   = epksi*ComputeDMicroscopicCrossSection(kineticEnergy,Z,epksi);
     433
     434  } while(func2 < func1*G4UniformRand());
     435
    553436  G4double gEnergy = epksi;
    554437
    555   // sample angle
     438  // ===== sample angle =====
     439
    556440  G4double gam  = totalEnergy/mass;
    557   G4double rmax = gam*min(1.0, totalEnergy/gEnergy - 1.0);
    558   rmax *= rmax;
    559   G4double x = G4UniformRand()*rmax/(1.0 + rmax);
     441  G4double rmax = gam*std::min(1.0, totalEnergy/gEnergy - 1.0);
     442  G4double rmax2= rmax*rmax;
     443  G4double x = G4UniformRand()*rmax2/(1.0 + rmax2);
    560444
    561445  G4double theta = sqrt(x/(1.0 - x))/gam;
    … …  
    577461
    578462  // save secondary
    579   G4DynamicParticle* aGamma = new G4DynamicParticle(theGamma,gDirection,gEnergy);
     463  G4DynamicParticle* aGamma =
     464    new G4DynamicParticle(theGamma,gDirection,gEnergy);
    580465  vdp->push_back(aGamma);
    581466}

  • trunk/source/processes/electromagnetic/muons/src/G4MuIonisation.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuIonisation.cc,v 1.54 2007/05/22 17:35:58 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuIonisation.cc,v 1.59 2009/02/26 11:04:20 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8686#include "G4MuBetheBlochModel.hh"
    8787#include "G4UniversalFluctuation.hh"
     88#include "G4IonFluctuations.hh"
    8889#include "G4BohrFluctuations.hh"
    8990#include "G4UnitsTable.hh"
    … …  
    99100    isInitialised(false)
    100101{
    101   SetStepFunction(0.2, 1*mm);
    102   SetIntegral(true);
    103   SetVerboseLevel(1);
     102  //  SetStepFunction(0.2, 1*mm);
     103  //SetIntegral(true);
     104  //SetVerboseLevel(1);
     105  SetProcessSubType(fIonisation);
    104106}
    105107
    … …  
    108110G4MuIonisation::~G4MuIonisation()
    109111{}
     112
     113//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     114
     115G4bool G4MuIonisation::IsApplicable(const G4ParticleDefinition& p)
     116{
     117  return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV);
     118}
     119
     120//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     121
     122G4double G4MuIonisation::MinPrimaryEnergy(const G4ParticleDefinition*,
     123                                          const G4Material*,
     124                                          G4double cut)
     125{
     126  G4double x = 0.5*cut/electron_mass_c2;
     127  G4double g = x*ratio + std::sqrt((1. + x)*(1. + x*ratio*ratio));
     128  return mass*(g - 1.0);
     129}
    110130
    111131//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    122142    SetSecondaryParticle(G4Electron::Electron());
    123143
    124     flucModel = new G4UniversalFluctuation();
     144    // Bragg peak model
     145    if (!EmModel(1)) SetEmModel(new G4BraggModel(),1);
     146    EmModel(1)->SetLowEnergyLimit(MinKinEnergy());
     147    EmModel(1)->SetHighEnergyLimit(0.2*MeV);
     148    AddEmModel(1, EmModel(1), new G4IonFluctuations());
    125149
    126     G4VEmModel* em = new G4BraggModel();
    127     em->SetLowEnergyLimit(0.1*keV);
    128     em->SetHighEnergyLimit(0.2*MeV);
    129     AddEmModel(1, em, flucModel);
    130     G4VEmModel* em1 = new G4BetheBlochModel();
    131     em1->SetLowEnergyLimit(0.2*MeV);
    132     em1->SetHighEnergyLimit(1.0*GeV);
    133     AddEmModel(2, em1, flucModel);
    134     G4VEmModel* em2 = new G4MuBetheBlochModel();
    135     em2->SetLowEnergyLimit(1.0*GeV);
    136     em2->SetHighEnergyLimit(100.0*TeV);
    137     AddEmModel(3, em2, flucModel);
     150    // high energy fluctuation model
     151    if (!FluctModel()) SetFluctModel(new G4UniversalFluctuation());
     152
     153    // moderate energy model
     154    if (!EmModel(2)) SetEmModel(new G4BetheBlochModel(),2);
     155    EmModel(2)->SetLowEnergyLimit(0.2*MeV);
     156    EmModel(2)->SetHighEnergyLimit(1.0*GeV);
     157    AddEmModel(2, EmModel(2), FluctModel());
     158
     159    // high energy model
     160    if (!EmModel(3)) SetEmModel(new G4MuBetheBlochModel(),3);
     161    EmModel(3)->SetLowEnergyLimit(1.0*GeV);
     162    EmModel(3)->SetHighEnergyLimit(MaxKinEnergy());
     163    AddEmModel(3, EmModel(3), FluctModel());
    138164
    139165    ratio = electron_mass_c2/mass;
    … …  
    145171
    146172void G4MuIonisation::PrintInfo()
    147 {
    148   G4cout << "      Bether-Bloch model for E > 0.2 MeV, "
    149          << "parametrisation of Bragg peak below, "
    150          << G4endl;
    151   G4cout << "      radiative corrections for E > 1 GeV" << G4endl;
    152 }
     173{}
    153174
    154175//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

  • trunk/source/processes/electromagnetic/muons/src/G4MuMultipleScattering.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuMultipleScattering.cc,v 1.3 2007/11/09 19:48:10 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuMultipleScattering.cc,v 1.12 2008/10/16 13:37:04 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -----------------------------------------------------------------------------
    … …  
    4747
    4848#include "G4MuMultipleScattering.hh"
    49 #include "G4MuMscModel.hh"
     49#include "G4WentzelVIModel.hh"
    5050#include "G4MscStepLimitType.hh"
    5151
    … …  
    6161  samplez           = false ;
    6262  isInitialized     = false; 
    63   SetRangeFactor(0.04);
     63  SetRangeFactor(0.2);
    6464  SetLateralDisplasmentFlag(true);
    6565}
    … …  
    8484  if(isInitialized) {
    8585
    86     if (p->GetParticleType() != "nucleus") {
     86    if (p->GetParticleType() != "nucleus" && p->GetPDGMass() < GeV) {
    8787      mscModel->SetStepLimitType(StepLimitType());
    8888      mscModel->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
    89       //mscModel->SetThetaLimit(thetaLimit);
    9089      mscModel->SetRangeFactor(RangeFactor());
    9190    }
     91    mscModel->SetPolarAngleLimit(PolarAngleLimit());
    9292    return;
    9393  }
    9494
    95   if (p->GetParticleType() == "nucleus") {
     95  if (p->GetParticleType() == "nucleus" || p->GetPDGMass() > GeV) {
    9696    SetLateralDisplasmentFlag(false);
    9797    SetBuildLambdaTable(false);
    98     //    SetRangeFactor(0.2);
    9998  }
    10099
    101   // initialisation of parameters
    102   //  G4String part_name = p->GetParticleName();
    103   mscModel = new G4MuMscModel(RangeFactor(),thetaLimit);
     100  // initialisation of the model
     101
     102  mscModel = new G4WentzelVIModel();
     103  mscModel->SetStepLimitType(StepLimitType());
    104104  mscModel->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
     105  mscModel->SetRangeFactor(RangeFactor());
     106  mscModel->SetPolarAngleLimit(PolarAngleLimit());
     107  mscModel->SetLowEnergyLimit(MinKinEnergy());
     108  mscModel->SetHighEnergyLimit(MaxKinEnergy());
    105109
    106110  AddEmModel(1,mscModel);
    … …  
    112116void G4MuMultipleScattering::PrintInfo()
    113117{
    114   G4cout << "      Boundary/stepping algorithm is active with RangeFactor= "
    115          << RangeFactor()
    116          << "  Step limit type " << StepLimitType()
    117         << G4endl;
     118  G4cout << "      RangeFactor= " << RangeFactor()
     119         << ", step limit type: " << StepLimitType()
     120         << ", lateralDisplacement: " << LateralDisplasmentFlag()
     121        << G4endl;
    118122}
    119123

  • trunk/source/processes/electromagnetic/muons/src/G4MuPairProduction.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuPairProduction.cc,v 1.48 2007/05/22 17:35:58 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuPairProduction.cc,v 1.52 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8989    lowestKinEnergy(1.*GeV),
    9090    isInitialised(false)
    91 {}
     91{
     92  SetProcessSubType(fPairProdByCharged);
     93}
    9294
    9395//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    98100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    99101
    100 void G4MuPairProduction::InitialiseEnergyLossProcess(const G4ParticleDefinition* part,
    101                                                      const G4ParticleDefinition*)
     102G4bool G4MuPairProduction::IsApplicable(const G4ParticleDefinition& p)
     103{
     104  return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV);
     105}
     106
     107//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     108
     109G4double G4MuPairProduction::MinPrimaryEnergy(const G4ParticleDefinition*,
     110                                              const G4Material*,
     111                                              G4double)
     112{
     113  return lowestKinEnergy;
     114}
     115
     116//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     117
     118void G4MuPairProduction::InitialiseEnergyLossProcess(
     119                         const G4ParticleDefinition* part,
     120                         const G4ParticleDefinition*)
    102121{
    103122  if (!isInitialised) {
    … …  
    110129    G4MuPairProductionModel* em = new G4MuPairProductionModel();
    111130    em->SetLowestKineticEnergy(lowestKinEnergy);
    112     G4VEmFluctuationModel* fm = new G4UniversalFluctuation();
    113     em->SetLowEnergyLimit(0.1*keV);
    114     em->SetHighEnergyLimit(100.0*TeV);
     131    G4VEmFluctuationModel* fm = 0;
     132    em->SetLowEnergyLimit(MinKinEnergy());
     133    em->SetHighEnergyLimit(MaxKinEnergy());
    115134    AddEmModel(1, em, fm);
    116135  }
    … …  
    120139
    121140void G4MuPairProduction::PrintInfo()
    122 {
    123   G4cout << "      Parametrised model "
    124          << G4endl;
    125 }
     141{}
    126142
    127143//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

  • trunk/source/processes/electromagnetic/muons/src/G4MuPairProductionModel.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4MuPairProductionModel.cc,v 1.35 2007/10/11 13:52:04 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4MuPairProductionModel.cc,v 1.40 2009/02/20 14:48:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    104104                                                 const G4String& nam)
    105105  : G4VEmModel(nam),
    106   minPairEnergy(4.*electron_mass_c2),
    107   lowestKinEnergy(1.*GeV),
    108   factorForCross(4.*fine_structure_const*fine_structure_const
     106    particle(0),
     107    factorForCross(4.*fine_structure_const*fine_structure_const
    109108                   *classic_electr_radius*classic_electr_radius/(3.*pi)),
    110109    sqrte(sqrt(exp(1.))),
    111110    currentZ(0),
    112     particle(0),
     111    fParticleChange(0),
     112    minPairEnergy(4.*electron_mass_c2),
     113    lowestKinEnergy(1.*GeV),
    113114    nzdat(5),
    114115    ntdat(8),
    … …  
    118119    ymax(0.),
    119120    dy((ymax-ymin)/nbiny),
    120     ignoreCut(false),
    121121    samplingTablesAreFilled(false)
    122122{
    123123  SetLowEnergyLimit(minPairEnergy);
     124  nist = G4NistManager::Instance();
    124125
    125126  theElectron = G4Electron::Electron();
    … …  
    144145//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    145146
    146 void G4MuPairProductionModel::SetParticle(const G4ParticleDefinition* p)
    147 {
    148   if(!particle) {
    149     particle = p;
    150     particleMass = particle->GetPDGMass();
    151   }
     147G4double G4MuPairProductionModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
     148                                                     G4double kineticEnergy)
     149{
     150  G4double maxPairEnergy = kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
     151  return maxPairEnergy;
    152152}
    153153
    … …  
    161161    MakeSamplingTables();
    162162  }
    163   if(pParticleChange) {
    164     if(ignoreCut) {
    165       gParticleChange =
    166         reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
    167       fParticleChange = 0;
    168     } else {
     163  if(!fParticleChange) {
     164    if(pParticleChange)
    169165      fParticleChange =
    170166        reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
    171       gParticleChange = 0;
    172     }
    173   } else {
    174     fParticleChange = new G4ParticleChangeForLoss();
    175     gParticleChange = 0;
     167    else
     168      fParticleChange = new G4ParticleChangeForLoss();
    176169  }
    177170}
    … …  
    186179{
    187180  G4double dedx = 0.0;
    188   if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy
    189       || ignoreCut)
     181  if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
    190182    return dedx;
    191183
    … …  
    216208  G4double loss = 0.0;
    217209
    218   G4double cut  = min(cutEnergy,tmax);
     210  G4double cut = std::min(cutEnergy,tmax);
    219211  if(cut <= minPairEnergy) return loss;
    220212
    … …  
    251243                                           G4double Z,
    252244                                           G4double cut)
    253 
    254 {
    255   G4double cross = 0. ;
    256 
     245{
     246  G4double cross = 0.;
    257247  SetCurrentElement(Z);
    258248  G4double tmax = MaxSecondaryEnergy(particle, tkin);
    259 
    260249  if (tmax <= cut) return cross;
    261250
    … …  
    400389                                                 G4double Z, G4double,
    401390                                                 G4double cutEnergy,
    402                                                  G4double)
    403 {
    404   G4double cut  = max(minPairEnergy,cutEnergy);
    405   if(ignoreCut) cut = minPairEnergy;
    406   G4double cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
    407   return cross;
    408 }
    409 
    410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    411 
    412 G4double G4MuPairProductionModel::CrossSectionPerVolume(
    413                                                const G4Material* material,
    414                                                const G4ParticleDefinition*,
    415                                                      G4double kineticEnergy,
    416                                                      G4double cutEnergy,
    417                                                      G4double maxEnergy)
     391                                                 G4double maxEnergy)
    418392{
    419393  G4double cross = 0.0;
    420394  if (kineticEnergy <= lowestKinEnergy) return cross;
    421395
    422   maxEnergy += particleMass;
    423 
    424   const G4ElementVector* theElementVector = material->GetElementVector();
    425   const G4double* theAtomNumDensityVector = material->
    426                                                     GetAtomicNumDensityVector();
    427 
    428   for (size_t i=0; i<material->GetNumberOfElements(); i++) {
    429     G4double Z = (*theElementVector)[i]->GetZ();
    430     SetCurrentElement(Z);
    431     G4double tmax = min(maxEnergy,MaxSecondaryEnergy(particle, kineticEnergy));
    432     G4double cut  = max(minPairEnergy,cutEnergy);
    433     if(ignoreCut) cut = minPairEnergy;
    434     if(cut < tmax) {
    435       G4double cr = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut)
    436                   - ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
    437 
    438       cross += theAtomNumDensityVector[i] * cr;
    439     }
     396  SetCurrentElement(Z);
     397  G4double tmax = std::min(maxEnergy, kineticEnergy);
     398  G4double cut  = std::min(cutEnergy, kineticEnergy);
     399  if(cut < minPairEnergy) cut = minPairEnergy;
     400  if (cut >= tmax) return cross;
     401
     402  cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
     403  if(tmax < kineticEnergy) {
     404    cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
    440405  }
    441406  return cross;
    … …  
    451416    SetCurrentElement(Z);
    452417
    453     for (G4int it=0; it<ntdat; it++)
    454     {
     418    for (G4int it=0; it<ntdat; it++) {
     419
    455420      G4double kineticEnergy = tdat[it];
    456421      G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
    457 
     422      // G4cout << "Z= " << currentZ << " z13= " << z13
     423      //<< " mE= " << maxPairEnergy << G4endl;
    458424      G4double CrossSection = 0.0 ;
    459425
    460       G4double y = ymin - 0.5*dy ;
    461       G4double yy = ymin - dy ;
    462       G4double x = exp(y);
    463       G4double fac = exp(dy);
    464       G4double dx = exp(yy)*(fac - 1.0);
    465 
    466       G4double c = log(maxPairEnergy/minPairEnergy);
    467 
    468       for (G4int i=0 ; i<nbiny; i++)
    469       {
    470         y += dy ;
    471         if(c > 0.0) {
    472           x *= fac;
    473           dx*= fac;
    474           G4double ep = minPairEnergy*exp(c*x) ;
    475           CrossSection += ep*dx*ComputeDMicroscopicCrossSection(
    476                                       kineticEnergy, Z, ep);
    477         }
    478         ya[i] = y;
    479         proba[iz][it][i] = CrossSection;
     426      if(maxPairEnergy > minPairEnergy) {
     427
     428        G4double y = ymin - 0.5*dy ;
     429        G4double yy = ymin - dy ;
     430        G4double x = exp(y);
     431        G4double fac = exp(dy);
     432        G4double dx = exp(yy)*(fac - 1.0);
     433
     434        G4double c = log(maxPairEnergy/minPairEnergy);
     435
     436        for (G4int i=0 ; i<nbiny; i++) {
     437          y += dy ;
     438          if(c > 0.0) {
     439            x *= fac;
     440            dx*= fac;
     441            G4double ep = minPairEnergy*exp(c*x) ;
     442            CrossSection +=
     443              ep*dx*ComputeDMicroscopicCrossSection(kineticEnergy, Z, ep);
     444          }
     445          ya[i] = y;
     446          proba[iz][it][i] = CrossSection;
     447        }
     448       
     449      } else {
     450        for (G4int i=0 ; i<nbiny; i++) {
     451          proba[iz][it][i] = CrossSection;
     452        }
    480453      }
    481454
    482455      ya[nbiny]=ymax;
    483 
    484456      proba[iz][it][nbiny] = CrossSection;
    485457
    … …  
    515487  G4double maxEnergy     = std::min(tmax, maxPairEnergy);
    516488  G4double minEnergy     = std::max(tmin, minPairEnergy);
    517   if(ignoreCut)minEnergy = minPairEnergy;
     489
    518490  if(minEnergy >= maxEnergy) return;
    519491  //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
    … …  
    618590  // primary change
    619591  kineticEnergy -= (ElectronEnergy + PositronEnergy);
    620   if(fParticleChange)
    621     fParticleChange->SetProposedKineticEnergy(kineticEnergy);
    622   else
    623     gParticleChange->SetProposedKineticEnergy(kineticEnergy);
     592  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
    624593
    625594  vdp->push_back(aParticle1);
    … …  
    655624    G4double maxPairEnergy = MaxSecondaryEnergy(particle,kinEnergy);
    656625    G4double minEnergy     = std::max(tmin, minPairEnergy);
    657     if(ignoreCut)minEnergy = minPairEnergy;
    658626
    659627    G4int iz;
Note: See TracChangeset for help on using the changeset viewer.