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

update ti head

Location:
trunk/source/processes/electromagnetic
Files:
98 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/electromagnetic/highenergy/History

    r1315 r1340  
    1 $Id: History,v 1.36 2010/06/04 10:23:31 vnivanch Exp $
     1$Id: History,v 1.37 2010/10/26 15:40:03 vnivanch Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     2026 October 10: V.Ivanchenko (emhighenergy-V09-03-02)
     21- Fixed problem reported by the Coverity tools (mainly pedantic
     22  initialisation)
     23- Added G4mplIonisationWithDeltaModel which is substituted
     24  G4mplIonisationModel by default - delat-ray production is required
     25  both by ATLAS and CMS
    1926
    202704 March 10: V.Ivanchenko (emhighenergy-V09-03-01)

  • trunk/source/processes/electromagnetic/highenergy/include/G4GammaConversionToMuons.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4GammaConversionToMuons.hh,v 1.2 2006/06/29 19:32:18 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4GammaConversionToMuons.hh,v 1.3 2010/10/26 14:15:40 vnivanch Exp $
     28// GEANT4 tag $Name: emhighenergy-V09-03-02 $
    2929//
    3030//         ------------ G4GammaConversionToMuons physics process ------
    … …  
    132132     G4double HighestEnergyLimit ;    // high energy limit of the tables
    133133
    134      G4double fminimalEnergy;         // minimalEnergy of produced particles
    135 
    136134     G4double MeanFreePath;           // actual MeanFreePath (current medium)
    137135     G4double CrossSecFactor;         // factor to artificially increase

  • trunk/source/processes/electromagnetic/highenergy/include/G4hhIonisation.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4hhIonisation.hh,v 1.6 2009/02/20 16:38:33 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4hhIonisation.hh,v 1.7 2010/10/26 14:15:40 vnivanch Exp $
     27// GEANT4 tag $Name: emhighenergy-V09-03-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9090  G4double   mass;
    9191  G4double   ratio;
    92   G4double   minKinEnergy;
    93   G4double   maxKinEnergy;
    9492
    9593  const G4ParticleDefinition* theParticle;
    … …  
    9896
    9997  G4bool                      isInitialised;
    100 
    101   G4double                    eth;
    102 
    10398};
    10499

  • trunk/source/processes/electromagnetic/highenergy/include/G4mplIonisationModel.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4mplIonisationModel.hh,v 1.7 2009/02/20 16:38:33 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4mplIonisationModel.hh,v 1.8 2010/10/26 15:40:03 vnivanch Exp $
     27// GEANT4 tag $Name: emhighenergy-V09-03-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    109109  G4double beta2lim;
    110110  G4double bg2lim;
    111   G4double factlow;
    112111  G4double chargeSquare;
    113112  G4double dedxlim;
    114113  G4int    nmpl;
    115114  G4double pi_hbarc2_over_mc2;
    116   G4double approxConst;
    117115};
    118116

  • trunk/source/processes/electromagnetic/highenergy/src/G4GammaConversionToMuons.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4GammaConversionToMuons.cc,v 1.7 2008/10/16 14:29:48 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4GammaConversionToMuons.cc,v 1.8 2010/10/26 14:15:40 vnivanch Exp $
     28// GEANT4 tag $Name: emhighenergy-V09-03-02 $
    2929//
    3030//         ------------ G4GammaConversionToMuons physics process ------
    … …  
    5353{
    5454  SetProcessSubType(15);
     55  MeanFreePath = DBL_MAX;
    5556}
    5657

  • trunk/source/processes/electromagnetic/highenergy/src/G4eeToHadronsModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4eeToHadronsModel.cc,v 1.9 2008/07/10 18:06:39 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4eeToHadronsModel.cc,v 1.10 2010/10/26 14:15:40 vnivanch Exp $
     27// GEANT4 tag $Name: emhighenergy-V09-03-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    7676{
    7777  theGamma = G4Gamma::Gamma();
     78  highKinEnergy = HighEnergyLimit();
     79  lowKinEnergy  = LowEnergyLimit();
     80  emin = lowKinEnergy;
     81  emax = highKinEnergy;
     82  peakKinEnergy = highKinEnergy;
     83  epeak = emax;
    7884}
    7985
    … …  
    9298                                    const G4DataVector&)
    9399{
    94   if(isInitialised) return;
     100  if(isInitialised) { return; }
    95101  isInitialised  = true;
    96102

  • trunk/source/processes/electromagnetic/highenergy/src/G4eeToHadronsMultiModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4eeToHadronsMultiModel.cc,v 1.8 2009/04/12 17:48:21 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4eeToHadronsMultiModel.cc,v 1.9 2010/10/26 14:15:40 vnivanch Exp $
     27// GEANT4 tag $Name: emhighenergy-V09-03-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6868    verbose(ver),
    6969    isInitialised(false)
    70 {}
     70{
     71  thKineticEnergy  = DBL_MAX;
     72  maxKineticEnergy = 1.2*GeV;
     73}
    7174
    7275//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    9093  if(!isInitialised) {
    9194    isInitialised = true;
    92 
    93     thKineticEnergy  = DBL_MAX;
    94     maxKineticEnergy = 1.2*GeV;
    9595
    9696    cross = new G4eeCrossSections();

  • trunk/source/processes/electromagnetic/highenergy/src/G4hhIonisation.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4hhIonisation.cc,v 1.10 2010/06/04 10:23:31 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4hhIonisation.cc,v 1.11 2010/10/26 14:15:40 vnivanch Exp $
     27// GEANT4 tag $Name: emhighenergy-V09-03-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    102102                                                 const G4ParticleDefinition* bpart)
    103103{
    104   if(isInitialised) return;
     104  if(isInitialised) { return; }
    105105
    106106  theParticle = part;
    107   if(bpart) G4cout << "G4hhIonisation::InitialiseEnergyLossProcess WARNING: no "
    108                    << "base particle should be defined for the process "
    109                    << GetProcessName() << G4endl;
    110 
     107  if(bpart) {
     108    G4cout << "G4hhIonisation::InitialiseEnergyLossProcess WARNING: no "
     109           << "base particle should be defined for the process "
     110           << GetProcessName() << G4endl;
     111  }
    111112  SetBaseParticle(0);
    112113  SetSecondaryParticle(G4Electron::Electron());
    … …  
    114115  mass  = theParticle->GetPDGMass();
    115116  ratio = electron_mass_c2/mass;
    116   eth = 2.0*MeV*mass/proton_mass_c2;
     117  G4double eth = 2*MeV*mass/proton_mass_c2;
    117118  flucModel = new G4BohrFluctuations();
    118119
    119120  G4int nm = 1;
    120121
    121   minKinEnergy = MinKinEnergy();
     122  G4double minKinEnergy = MinKinEnergy();
    122123
    123124  if(eth > minKinEnergy) {

  • trunk/source/processes/electromagnetic/highenergy/src/G4mplIonisation.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4mplIonisation.cc,v 1.10 2010/03/28 16:45:38 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4mplIonisation.cc,v 1.11 2010/10/26 15:40:03 vnivanch Exp $
     27// GEANT4 tag $Name: emhighenergy-V09-03-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4949#include "G4Electron.hh"
    5050#include "G4mplIonisationModel.hh"
     51#include "G4mplIonisationWithDeltaModel.hh"
    5152#include "G4BohrFluctuations.hh"
    5253
    … …  
    8586                                                  const G4ParticleDefinition*)
    8687{
    87   if(isInitialised) return;
     88  if(isInitialised) { return; }
    8889
    8990  SetBaseParticle(0);
    9091  SetSecondaryParticle(G4Electron::Electron());
    9192
    92   G4mplIonisationModel* ion  = new G4mplIonisationModel(magneticCharge,"PAI");
     93  G4mplIonisationWithDeltaModel* ion =
     94    new G4mplIonisationWithDeltaModel(magneticCharge,"PAI");
    9395  ion->SetLowEnergyLimit(MinKinEnergy());
    9496  ion->SetHighEnergyLimit(MaxKinEnergy());

  • trunk/source/processes/electromagnetic/highenergy/src/G4mplIonisationModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4mplIonisationModel.cc,v 1.7 2009/04/12 17:35:41 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4mplIonisationModel.cc,v 1.8 2010/10/26 15:40:03 vnivanch Exp $
     27// GEANT4 tag $Name: emhighenergy-V09-03-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    7474{
    7575  nmpl         = G4int(abs(magCharge) * 2 * fine_structure_const + 0.5);
    76   if(nmpl > 6)      nmpl = 6;
    77   else if(nmpl < 1) nmpl = 1;
     76  if(nmpl > 6)      { nmpl = 6; }
     77  else if(nmpl < 1) { nmpl = 1; }
    7878  pi_hbarc2_over_mc2 = pi * hbarc * hbarc / electron_mass_c2;
    7979  chargeSquare = magCharge * magCharge;
    8080  dedxlim = 45.*nmpl*nmpl*GeV*cm2/g;
    8181  fParticleChange = 0;
     82  mass = 0.0;
    8283}
    8384
    … …  
    9495  monopole = p;
    9596  mass     = monopole->GetPDGMass();
    96   if(!fParticleChange) fParticleChange = GetParticleChangeForLoss();
     97  if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
    9798}
    9899

  • trunk/source/processes/electromagnetic/lowenergy/History

    r1337 r1340  
    1 $Id: History,v 1.442 2010/06/15 08:04:11 gcosmo Exp $
     1$Id: History,v 1.462 2010/11/04 14:52:17 sincerti Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     2004.11.2010, S. Incerti, tag emlowen-V09-03-54
     21            - removed warnings in Rudd ionization classes (SI)
     22            - G4IonParametrisedLossModel - (VI) moved few virtual methods from
     23            inline to source, minor cleanup of initialisation 
     24
     2503.11.2010, S. Incerti, tag emlowen-V09-03-53
     26            restricted momentum conservation to electrons
     27            in G4DNA ionisation
     28
     2903.11.2010, S. Incerti, tag emlowen-V09-03-52
     30            new preliminary Geant4-DNA ionisation class for ions by Z. Francis
     31            to be used with G4LEDATA 6.18
     32
     3317.10.2010, S. Incerti, tag emlowen-V09-03-51
     34            extended low energy coverage of G4DNA electron models
     35            to be used with G4LEDATA 6.17
     36
     3714.10.2010, V. Ivanchenko, tag emlowen-V09-03-50
     38            G4GeneratorBS - optimise computations to speedup, fixed comments
     39            G4GeneratorBN - fixed comments
     40            G4VBremAngularDistribution - moved to utils
     41            G4ModifiedTsai - moved to standard
     42
     4313.10.2010, L. Pandola, tag emlowen-V09-03-49
     44            Update G4LivermoreIonisationModel to produce fluorescence AlongStep
     45            only if above the production cuts. Stricter check for energy
     46            conservation
     47
     4808.10.2010, S. Incerti, tag emlowen-V09-03-48
     49            Added new excitation model for H
     50
     5115.09.2010, S. Incerti, tag emlowen-V09-03-47
     52            Added protection in xs file opening for G4DNA Sanche excitation
     53
     5415.09.2010, S. Incerti, tag emlowen-V09-03-46
     55            Corrected data file names in G4DNA Melton and Sanche       
     56
     5708.09.2010, S. Incerti, tag emlowen-V09-03-45
     58            Updated high energy limits of G4DNAExcitation and G4DNAIonisation
     59
     6008.09.2010, S. Incerti, tag emlowen-V09-03-44
     61            Set high energy limit of G4DNAScreenedRutherfordModel to 1 MeV
     62
     6308.09.2010, S. Incerti, tag emlowen-V09-03-43
     64            Added new G4DNA processes and models for vib. exc & attachment
     65            Provided by Z. Francis et al. - Appl. Rad. Isot. (2010)
     66            http://dx.doi.org/10.1016/j.apradiso.2010.08.011
     67            to be used with G4LEDATA 6.16
     68
     6908.09.2010, S. Incerti, tag emlowen-V09-03-42
     70            Decreased low energy limit of G4DNAScreenedRutherfordModel
     71
     7205.09.2010, S. Incerti, tag emlowen-V09-03-41
     73            Bugzilla 1120
     74            Modified G4PhotoElectricAngularGeneratorSauterGavrila.cc
     75            as proposed by J. Goldberg
     76
     7725.08.2010, S. Incerti, tag emlowen-V09-03-40
     78            - updated & extended Rudd and Miller & Green models
     79            - to be used with G4LEDATA 6.15
     80
     8125.08.2010, S. Incerti, tag emlowen-V09-03-39
     82            -Adapted all high energy limits of G4DNA electron models
     83
     8424.08.2010, S. Incerti, tag emlowen-V09-03-38
     85            -Changed low energy limit of G4DNA elastic scattering models for e-
     86            -Switched default excitation model for e- to Born
     87            -to be used with G4LEDATA 6.14
     88
     8928.07.2010, L. Pandola, tag emlowen-V09-03-37
     90            First full version of G4Penelope08IonisationModel, model for e+/e-
     91            ionisation according to Penelope v2008. Still beta version.
     92
     9326.07.2010, L. Pandola, tag emlowen-V09-03-36
     94            Added class G4PenelopeCrossSection to store/handle cross sections
     95            (and higher momenta, like stopping powers) for the updated
     96            Penelope08 e+/e- models (ionisation and bremsstrahlung).
    1997
    209815.06.2010, G. Cosmo, tag emlowen-V09-03-35

  • trunk/source/processes/electromagnetic/lowenergy/include/G4DNAScreenedRutherfordElasticModel.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DNAScreenedRutherfordElasticModel.hh,v 1.2 2010/01/07 18:10:19 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4DNAScreenedRutherfordElasticModel.hh,v 1.3 2010/09/08 13:39:11 sincerti Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929
    … …  
    6767  G4double killBelowEnergy; 
    6868  G4double lowEnergyLimit; 
    69   G4double lowEnergyLimitOfModel; 
    70   G4double intermediateEnergyLimit;
     69  G4double intermediateEnergyLimit;
    7170  G4double highEnergyLimit;
    7271  G4bool isInitialised;

  • trunk/source/processes/electromagnetic/lowenergy/include/G4Generator2BN.hh

    r819 r1340  
    2424// ********************************************************************
    2525//
     26// $Id: G4Generator2BN.hh,v 1.4 2010/10/14 14:00:29 vnivanch Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2628//
    2729// -------------------------------------------------------------------
    … …  
    6365public:
    6466
    65   G4Generator2BN(const G4String& name);
     67  G4Generator2BN(const G4String& name = "");
    6668
    67   ~G4Generator2BN();
     69  virtual ~G4Generator2BN();
    6870
    6971  G4double PolarAngle(const G4double initial_energy,

  • trunk/source/processes/electromagnetic/lowenergy/include/G4Generator2BS.hh

    r819 r1340  
    2424// ********************************************************************
    2525//
     26// $Id: G4Generator2BS.hh,v 1.5 2010/10/14 14:00:29 vnivanch Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2628//
    2729// -------------------------------------------------------------------
    … …  
    3941//
    4042// Modifications:
    41 // 02 Jun 2003                           First implementation acording with new design
     43// 02 Jun 2003  First implementation acording with new design
     44// 12 Oct 2010  V.Ivanchenko moved RejectionFunction inline
    4245//               
    4346//
    4447// Class Description:
    4548//
    46 // Concrete class for Bremsstrahlung Angular Distribution Generation - 2BS Distribution
    47 // Further documentation available from http://www.ge.infn.it/geant4/lowE
     49// Concrete class for Bremsstrahlung Angular Distribution Generation
     50// 2BS Distribution
     51//
    4852
    4953// -------------------------------------------------------------------
    … …  
    5761#include "G4VBremAngularDistribution.hh"
    5862
     63class G4Pow;
     64
    5965class G4Generator2BS : public G4VBremAngularDistribution
    6066{
    … …  
    6268public:
    6369
    64   G4Generator2BS(const G4String& name);
     70  G4Generator2BS(const G4String& name="");
    6571
    66   ~G4Generator2BS();
     72  virtual ~G4Generator2BS();
    6773
    6874  G4double PolarAngle(const G4double initial_energy,
    … …  
    7480protected:
    7581
    76   G4double RejectionFunction(G4double value) const;
     82  inline G4double RejectionFunction(G4double value) const;
    7783
    7884private:
     85
    7986  G4double z;
    8087  G4double rejection_argument1, rejection_argument2, rejection_argument3;
    8188  G4double EnergyRatio;
    8289
     90  G4Pow* g4pow;
     91
    8392  // hide assignment operator
    84      G4Generator2BS & operator=(const  G4Generator2BS &right);
    85      G4Generator2BS(const  G4Generator2BS&);
     93  G4Generator2BS & operator=(const  G4Generator2BS &right);
     94  G4Generator2BS(const  G4Generator2BS&);
    8695
    8796};
    8897
     98inline G4double G4Generator2BS::RejectionFunction(G4double value) const
     99{
     100  G4double argument = (1+value)*(1+value);
     101  return (4+std::log(rejection_argument3+(z/argument)))*
     102    ((4*EnergyRatio*value/argument)-rejection_argument1)+rejection_argument2;
     103}
     104
    89105#endif
    90106

  • trunk/source/processes/electromagnetic/lowenergy/include/G4IonParametrisedLossModel.hh

    r1228 r1340  
    2424// ********************************************************************
    2525//
    26 //
     26// $Id: G4IonParametrisedLossModel.hh,v 1.8 2010/11/04 12:21:47 vnivanch Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2728//
    2829// ===========================================================================
    … …  
    136137
    137138   // Function, which computes the mean energy transfer rate to delta rays
    138    G4double DeltaRayMeanEnergyTransferRate(
     139   inline G4double DeltaRayMeanEnergyTransferRate(
    139140                                 const G4Material*, // Target Material
    140141                                 const G4ParticleDefinition*, // Projectile
    … …  
    183184   // Function checking the applicability of physics tables to ion-material
    184185   // combinations (Note: the energy range of tables is not checked)
    185    LossTableList::iterator IsApplicable(
     186   inline LossTableList::iterator IsApplicable(
    186187                      const G4ParticleDefinition*,  // Projectile (ion)
    187188                      const G4Material*);           // Target material
    … …  
    208209   
    209210   // Function for setting energy loss limit for stopping power integration
    210    void SetEnergyLossLimit(G4double ionEnergyLossLimit);
     211   inline void SetEnergyLossLimit(G4double ionEnergyLossLimit);
    211212
    212213 protected:
     214
     215   virtual
    213216   G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
    214217                               G4double);   // Kinetic energy of projectile
    … …  
    231234
    232235   // Function, which updates parameters concering particle properties
    233    void UpdateCache(
     236   inline void UpdateCache(
    234237                  const G4ParticleDefinition*);  // Projectile (ion)
    235238 
    … …  
    287290
    288291   // Pointer to particle change object, which is used to set e.g. the
    289    // energy loss due to nuclear stopping
     292   // energy loss and secondary delta-electron
     293   // used indicating if model is initialized 
    290294   G4ParticleChangeForLoss* particleChangeLoss;
    291 
    292    // Flag indicating if model is initialized (i.e. if
    293    // G4ParticleChangeForLoss was created)
    294    G4bool modelIsInitialised;
    295295
    296296   // ######################################################################

  • trunk/source/processes/electromagnetic/lowenergy/include/G4IonParametrisedLossModel.icc

    r1196 r1340  
    2424// ********************************************************************
    2525//
    26 //
     26// $Id: G4IonParametrisedLossModel.icc,v 1.7 2010/11/04 12:21:47 vnivanch Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2728//
    2829// ===========================================================================
    … …  
    4344//                               Minor bug fix in ComputeDEDXPerVolume (AL)
    4445//                20. 11. 2009 - Added set-method for energy loss limit (AL)
     46//                04. 11. 2010 - Moved virtual methods to the source (VI)
    4547//
    4648// Class description:
    … …  
    5456//
    5557// ===========================================================================
    56 
    5758
    5859inline G4double G4IonParametrisedLossModel::DeltaRayMeanEnergyTransferRate(
    … …  
    116117}
    117118
    118 
    119 inline
    120 G4double G4IonParametrisedLossModel::MaxSecondaryEnergy(
    121                              const G4ParticleDefinition* particle,
    122                              G4double kineticEnergy) {
    123 
    124   // ############## Maximum energy of secondaries ##########################
    125   // Function computes maximum energy of secondary electrons which are
    126   // released by an ion
    127   //
    128   // See Geant4 physics reference manual (version 9.1), section 9.1.1
    129   //
    130   // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
    131   //       C.Caso et al. (Part. Data Group), Europ. Phys. Jour. C 3 1 (1998).
    132   //       B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
    133   //
    134   // (Implementation adapted from G4BraggIonModel)
    135 
    136   if(particle != cacheParticle) UpdateCache(particle);
    137 
    138   G4double tau  = kineticEnergy/cacheMass;
    139   G4double tmax = 2.0 * electron_mass_c2 * tau * (tau + 2.) /
    140                   (1. + 2.0 * (tau + 1.) * cacheElecMassRatio +
    141                   cacheElecMassRatio * cacheElecMassRatio);
    142 
    143   return tmax;
    144 }
    145 
    146 
    147119inline
    148120void G4IonParametrisedLossModel::UpdateCache(
    … …  
    155127  cacheChargeSquare = q * q;
    156128}
    157 
    158 
    159 inline
    160 G4double G4IonParametrisedLossModel::GetChargeSquareRatio(
    161                              const G4ParticleDefinition* particle,
    162                              const G4Material* material,
    163                              G4double kineticEnergy) {    // Kinetic energy
    164 
    165   G4double chargeSquareRatio = corrections ->
    166                                      EffectiveChargeSquareRatio(particle,
    167                                                                 material,
    168                                                                 kineticEnergy);
    169   corrFactor = chargeSquareRatio *
    170                        corrections -> EffectiveChargeCorrection(particle,
    171                                                                 material,
    172                                                                 kineticEnergy);
    173   return corrFactor;
    174 }
    175 
    176 
    177 inline
    178 G4double G4IonParametrisedLossModel::GetParticleCharge(
    179                              const G4ParticleDefinition* particle,
    180                              const G4Material* material,
    181                              G4double kineticEnergy) {   // Kinetic energy
    182 
    183   return corrections -> GetParticleCharge(particle, material, kineticEnergy);
    184 }
    185 
    186129
    187130inline

  • trunk/source/processes/electromagnetic/lowenergy/src/G4DNABornExcitationModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DNABornExcitationModel.cc,v 1.9 2010/03/27 11:32:41 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4DNABornExcitationModel.cc,v 1.10 2010/08/24 13:51:06 sincerti Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929
    … …  
    7878  if (verboseLevel > 3)
    7979    G4cout << "Calling G4DNABornExcitationModel::Initialise()" << G4endl;
    80 
    81   // Energy limits
    82   // Energy limits
    8380
    8481  G4String fileElectron("dna/sigma_excitation_e_born");

  • trunk/source/processes/electromagnetic/lowenergy/src/G4DNABornIonisationModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DNABornIonisationModel.cc,v 1.16 2010/03/26 18:10:43 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4DNABornIonisationModel.cc,v 1.18 2010/11/03 12:22:36 sincerti Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929
    … …  
    376376    deltaDirection.rotateUz(primaryDirection);
    377377
    378     G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
    379 
    380     G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
    381     G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
    382     G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
    383     G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
    384     finalPx /= finalMomentum;
    385     finalPy /= finalMomentum;
    386     finalPz /= finalMomentum;
    387    
    388     G4ThreeVector direction;
    389     direction.set(finalPx,finalPy,finalPz);
    390    
    391     fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
     378    if (particle->GetDefinition() == G4Electron::ElectronDefinition())
     379    {
     380      G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
     381
     382      G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
     383      G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
     384      G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
     385      G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
     386      finalPx /= finalMomentum;
     387      finalPy /= finalMomentum;
     388      finalPz /= finalMomentum;
     389   
     390      G4ThreeVector direction;
     391      direction.set(finalPx,finalPy,finalPz);
     392   
     393      fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
     394    }
     395   
     396    else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
     397   
    392398    fParticleChangeForGamma->SetProposedKineticEnergy(k-bindingEnergy-secondaryKinetic);
    393399    fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
    … …  
    507513    G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
    508514    phi = twopi * G4UniformRand();
    509     cosTheta = std::sqrt(secKinetic / maxSecKinetic);
     515   
     516    // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
     517   
     518    // Restriction below 100 eV from Emfietzoglou (2000)
     519   
     520    if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
     521    else cosTheta = (2.*G4UniformRand())-1.;
     522       
    510523  }                     
    511524}

  • trunk/source/processes/electromagnetic/lowenergy/src/G4DNAChampionElasticModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DNAChampionElasticModel.cc,v 1.12 2010/04/08 17:29:08 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4DNAChampionElasticModel.cc,v 1.15 2010/10/17 11:28:51 sincerti Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929
    … …  
    4141{
    4242
    43   killBelowEnergy = 8.23*eV; // Minimum e- energy for energy loss by excitation
     43  killBelowEnergy = 0.025*eV; // Minimum e- energy for energy loss by excitation
    4444  lowEnergyLimit = 0 * eV;
    45   lowEnergyLimitOfModel = 7.4 * eV; // The model lower energy is 7.4 eV
    46   highEnergyLimit = 10 * MeV;
     45  lowEnergyLimitOfModel = 0.025 * eV;
     46  highEnergyLimit = 1. * MeV;
    4747  SetLowEnergyLimit(lowEnergyLimit);
    4848  SetHighEnergyLimit(highEnergyLimit);
    … …  
    139139
    140140    std::ostringstream eFullFileName;
    141     eFullFileName << path << "/dna/sigmadiff_cumulated_elastic_e_champion.dat";
     141    eFullFileName << path << "/dna/sigmadiff_cumulatedshort_elastic_e_champion.dat";
    142142    std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
    143143     

  • trunk/source/processes/electromagnetic/lowenergy/src/G4DNAElastic.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DNAElastic.cc,v 1.3 2009/03/04 13:28:49 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4DNAElastic.cc,v 1.4 2010/09/08 14:07:16 sincerti Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828
    2929#include "G4DNAElastic.hh"
    … …  
    6262    if(!Model()) SetModel(new G4DNAScreenedRutherfordElasticModel);
    6363    Model()->SetLowEnergyLimit(0*eV);
    64     Model()->SetHighEnergyLimit(10*MeV);
     64    Model()->SetHighEnergyLimit(1.*MeV);
    6565    AddEmModel(1, Model());
    6666  }

  • trunk/source/processes/electromagnetic/lowenergy/src/G4DNAExcitation.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DNAExcitation.cc,v 1.4 2010/03/27 11:32:41 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4DNAExcitation.cc,v 1.7 2010/10/08 08:53:17 sincerti Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828
    2929#include "G4DNAExcitation.hh"
    … …  
    5757       &p == G4Electron::Electron()
    5858    || &p == G4Proton::ProtonDefinition()
     59    || &p == instance->GetIon("hydrogen")
    5960    || &p == instance->GetIon("alpha++")
    6061    || &p == instance->GetIon("alpha+")
    … …  
    7778    {
    7879
    79       // First model
    80 
     80      // Emfietzoglou model
     81/*
    8182      if(!Model()) SetModel(new G4DNAEmfietzoglouExcitationModel);
    8283      Model()->SetLowEnergyLimit(8.23*eV);
    8384      Model()->SetHighEnergyLimit(10*MeV);
     85*/
     86      // Born model
    8487
    85       // Alternative model
    86 /*
    8788      if(!Model()) SetModel(new G4DNABornExcitationModel);
    8889      Model()->SetLowEnergyLimit(9*eV);
    8990      Model()->SetHighEnergyLimit(1*MeV);
    90 */
     91
    9192      AddEmModel(1, Model());   
    9293    }
    … …  
    106107    }
    107108
     109    if(name == "hydrogen")
     110    {
     111      if(!Model()) SetModel(new G4DNAMillerGreenExcitationModel);
     112      Model()->SetLowEnergyLimit(10*eV);
     113      Model()->SetHighEnergyLimit(500*keV);
     114   
     115      AddEmModel(1, Model());   
     116    }
     117
     118
    108119    if( name == "alpha" || name == "alpha+" || name == "helium" )
    109120    {
    110121      if(!Model()) SetModel(new G4DNAMillerGreenExcitationModel);
    111122      Model()->SetLowEnergyLimit(1*keV);
    112       Model()->SetHighEnergyLimit(10*MeV);
     123      Model()->SetHighEnergyLimit(400*MeV);
    113124
    114125      AddEmModel(1, Model());   

  • trunk/source/processes/electromagnetic/lowenergy/src/G4DNAGenericIonsManager.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DNAGenericIonsManager.cc,v 1.6 2009/06/10 13:32:36 mantero Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4DNAGenericIonsManager.cc,v 1.7 2010/11/03 10:44:26 sincerti Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828
    2929#include "G4DNAGenericIonsManager.hh"
    … …  
    7171 G4Ions *positronium2s;
    7272 
     73 G4Ions *carbon;
     74 G4Ions *nitrogen;
     75 G4Ions *oxygen;
     76 G4Ions *iron;
     77
     78 iron=     new G4Ions(
     79                        "iron",    52.5672*GeV,       0.0*MeV,  +26.0*eplus,
     80                        0,              +1,             0,
     81                        0,               0,             0,
     82                        "nucleus",              +26,            +56,           0,
     83                        true,                -1.0,             0,       
     84                        false,                "",               0,             
     85                        0.0);
     86
     87 oxygen=   new G4Ions(
     88                        "oxygen",    15.0074*GeV,       0.0*MeV,  +8.0*eplus,
     89                        0,              +1,             0,
     90                        0,               0,             0,
     91                        "nucleus",              +8,            +16,           0,
     92                        true,                -1.0,             0,       
     93                        false,                "",               0,             
     94                        0.0);
     95
     96
     97 nitrogen= new G4Ions(
     98                        "nitrogen",    13.132*GeV,       0.0*MeV,  +7.0*eplus,
     99                        0,              +1,             0,
     100                        0,               0,             0,
     101                        "nucleus",              +7,            +14,           0,
     102                        true,                -1.0,             0,       
     103                        false,                "",               0,             
     104                        0.0);
     105
     106 carbon=   new G4Ions(
     107                        "carbon",    11.267025440*GeV,       0.0*MeV,  +6.0*eplus,
     108                        0,              +1,             0,
     109                        0,               0,             0,
     110                        "nucleus",              +6,            +12,           0,
     111                        true,                -1.0,             0,       
     112                        false,                "",               0,             
     113                        0.0);
    73114 
    74  helium=     new G4Ions(
     115 helium=   new G4Ions(
    75116                        "helium",    3.727417*GeV,       0.0*MeV,  +0.0*eplus,
    76117                        0,              +1,             0,
    … …  
    81122                        0.0);
    82123
    83  alphaPlus=  new G4Ions("alpha+",    3.727417*GeV,       0.0*MeV,  +1.0*eplus,
     124 alphaPlus= new G4Ions("alpha+",    3.727417*GeV,       0.0*MeV,  +1.0*eplus,
    84125                               1,              +1,             0,
    85126                               0,               0,             0,
    … …  
    110151
    111152
    112  /*
    113  // molechules construction
    114 
    115  G4Ions* oxonium; // H3O -- it will become H3O+
    116  G4Ions* hydroxyl; // OH -- it will produce OH- too
    117  G4Ions* molHydrogen; // H2
    118  //G4Ions* hydroxide; // OH-
    119  G4Ions* hydroPeroxide; // H2O2
    120  G4Ions* water; // H2O -- it will become also H2O+
    121 
    122 
    123  G4double mass = 19.02*g/Avogadro - 11*electron_mass_c2;
    124 
    125  oxonium = new G4Ions("H3O",        mass,             0,  +11.0*eplus,
    126                       0,               0,             0,
    127                       0,               0,             0,
    128                       "molecule",      0,             0,           0,
    129                       true,         -1.0,             0,
    130                       false,          "",             0,             
    131                       0.0);
    132  
    133  mass = 17.00734*g/Avogadro - 9*electron_mass_c2;
    134  
    135  hydroxyl = new G4Ions("OH",        mass,             0,  +9.0*eplus,
    136                       0,               0,             0,
    137                       0,               0,             0,
    138                       "molecule",      0,             0,           0,
    139                       true,         -1.0,             0,
    140                       false,          "",             0,             
    141                       0.0);
    142 
    143  mass = 2.01588*g/Avogadro - 2*electron_mass_c2;
    144 
    145  molHydrogen = new G4Ions("H2",         mass,             0,  +2.0*eplus,
    146                           0,               0,             0,
    147                           0,               0,             0,
    148                           "molecule",      0,             0,           0,
    149                           true,         -1.0,             0,
    150                           false,          "",             0,             
    151                           0.0);
    152 
    153  mass = 34.01468*g/Avogadro - 18*electron_mass_c2;
    154 
    155  hydroPeroxide = new G4Ions("H2O2",      mass,             0, +18.0*eplus,
    156                            0,               0,             0,
    157                            0,               0,             0,
    158                            "molecule",      0,             0,           0,
    159                            true,         -1.0,             0,
    160                            false,          "",             0,             
    161                            0.0);
    162  
    163  mass = 18.015*g/Avogadro - 10*electron_mass_c2;
    164 
    165  water = new G4Ions("H2O",        mass,             0,  +10.0*eplus,
    166                     0,               0,             0,
    167                     0,               0,             0,
    168                     "molecule",      0,             0,           0,
    169                     true,         -1.0,             0,
    170                     false,          "",             0,             
    171                     0.0);
    172 
    173  map["H3O" ]  =oxonium;
    174  map["OH"  ]  =hydroxyl;
    175  map["H2"  ]  =molHydrogen;
    176  map["H2O2"]  =hydroPeroxide;
    177  map["H2O" ]  =water;
    178  */
    179 
    180 
    181153 map["helium"  ]=helium;
    182154 map["hydrogen"]=hydrogen;
    … …  
    185157 map["Ps-1s"   ]=positronium1s;
    186158 map["Ps-2s"   ]=positronium2s;
     159 map["carbon"  ]=carbon;
     160 map["nitrogen"]=nitrogen;
     161 map["oxygen"  ]=oxygen;
     162 map["iron"    ]=iron;
     163
    187164
    188165}

  • trunk/source/processes/electromagnetic/lowenergy/src/G4DNAIonisation.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DNAIonisation.cc,v 1.4 2009/11/02 17:00:11 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4DNAIonisation.cc,v 1.5 2010/09/08 14:30:45 sincerti Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828
    2929#include "G4DNAIonisation.hh"
    … …  
    9191      if(!Model(2)) SetModel(new G4DNABornIonisationModel,2);
    9292      Model(2)->SetLowEnergyLimit(500*keV);
    93       Model(2)->SetHighEnergyLimit(10*MeV);
     93      Model(2)->SetHighEnergyLimit(100*MeV);
    9494   
    9595      AddEmModel(1, Model(1));   
    … …  
    110110      if(!Model()) SetModel(new G4DNARuddIonisationModel);
    111111      Model()->SetLowEnergyLimit(0*keV);
    112       Model()->SetHighEnergyLimit(10*MeV);
     112      Model()->SetHighEnergyLimit(400*MeV);
    113113
    114114      AddEmModel(1, Model());   

  • trunk/source/processes/electromagnetic/lowenergy/src/G4DNAMillerGreenExcitationModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DNAMillerGreenExcitationModel.cc,v 1.9 2010/06/08 21:50:00 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4DNAMillerGreenExcitationModel.cc,v 1.11 2010/10/08 08:53:17 sincerti Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929
    … …  
    7474  instance = G4DNAGenericIonsManager::Instance();
    7575  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
     76  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
    7677  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
    7778  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
    … …  
    7980
    8081  G4String proton;
     82  G4String hydrogen;
    8183  G4String alphaPlusPlus;
    8284  G4String alphaPlus;
    … …  
    102104  }
    103105
     106  if (hydrogenDef != 0)
     107  {
     108    hydrogen = hydrogenDef->GetParticleName();
     109    lowEnergyLimit[hydrogen] = 10. * eV;
     110    highEnergyLimit[hydrogen] = 500. * keV;
     111   
     112    kineticEnergyCorrection[0] = 1.;
     113    slaterEffectiveCharge[0][0] = 0.;
     114    slaterEffectiveCharge[1][0] = 0.;
     115    slaterEffectiveCharge[2][0] = 0.;
     116    sCoefficient[0][0] = 0.;
     117    sCoefficient[1][0] = 0.;
     118    sCoefficient[2][0] = 0.;
     119  }
     120  else
     121  {
     122    G4Exception("G4DNAMillerGreenExcitationModel::Initialise: hydrogen is not defined");
     123   
     124  }
    104125  if (alphaPlusPlusDef != 0)
    105126  {
    106127    alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
    107128    lowEnergyLimit[alphaPlusPlus] = 1. * keV;
    108     highEnergyLimit[alphaPlusPlus] = 10. * MeV;
     129    highEnergyLimit[alphaPlusPlus] = 400. * MeV;
    109130
    110131    kineticEnergyCorrection[1] = 0.9382723/3.727417;
    … …  
    125146    alphaPlus = alphaPlusDef->GetParticleName();
    126147    lowEnergyLimit[alphaPlus] = 1. * keV;
    127     highEnergyLimit[alphaPlus] = 10. * MeV;
     148    highEnergyLimit[alphaPlus] = 400. * MeV;
    128149
    129150    kineticEnergyCorrection[2] = 0.9382723/3.727417;
    130151    slaterEffectiveCharge[0][2]=2.0;
    131     slaterEffectiveCharge[1][2]=1.15;
    132     slaterEffectiveCharge[2][2]=1.15;
     152
     153// Following values provided by M. Dingfelder
     154    slaterEffectiveCharge[1][2]=2.00;
     155    slaterEffectiveCharge[2][2]=2.00;
     156//
    133157    sCoefficient[0][2]=0.7;
    134158    sCoefficient[1][2]=0.15;
    … …  
    144168    helium = heliumDef->GetParticleName();
    145169    lowEnergyLimit[helium] = 1. * keV;
    146     highEnergyLimit[helium] = 10. * MeV;
     170    highEnergyLimit[helium] = 400. * MeV;
    147171   
    148172    kineticEnergyCorrection[3] = 0.9382723/3.727417;
    … …  
    153177    sCoefficient[1][3]=0.25;
    154178    sCoefficient[2][3]=0.25;
     179
    155180  }
    156181  else
    … …  
    163188    SetLowEnergyLimit(lowEnergyLimit[proton]);
    164189    SetHighEnergyLimit(highEnergyLimit[proton]);
     190  }
     191
     192  if (particle==hydrogenDef)
     193  {
     194    SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
     195    SetHighEnergyLimit(highEnergyLimit[hydrogen]);
    165196  }
    166197
    … …  
    231262      particleDefinition != G4Proton::ProtonDefinition()
    232263      &&
     264      particleDefinition != instance->GetIon("hydrogen")
     265      &&
    233266      particleDefinition != instance->GetIon("alpha++")
    234267      &&
    … …  
    272305
    273306      // add ONE or TWO electron-water excitation for alpha+ and helium
    274  
     307/* 
    275308      if ( particleDefinition == instance->GetIon("alpha+")
    276309           ||
    … …  
    278311         )
    279312      {
     313
    280314          G4DNAEmfietzoglouExcitationModel * excitationXS = new G4DNAEmfietzoglouExcitationModel();
     315          excitationXS->Initialise(G4Electron::ElectronDefinition());
    281316
    282317          G4double sigmaExcitation=0;
    … …  
    294329         
    295330          delete excitationXS;
     331
     332          // Alternative excitation model
     333
     334          G4DNABornExcitationModel * excitationXS = new G4DNABornExcitationModel();
     335          excitationXS->Initialise(G4Electron::ElectronDefinition());
     336
     337          G4double sigmaExcitation=0;
     338          G4double tmp=0;
     339
     340          if (k*0.511/3728 > 9*eV && k*0.511/3728 < 1*MeV ) sigmaExcitation =
     341            excitationXS->CrossSectionPerVolume(material,G4Electron::ElectronDefinition(),k*0.511/3728,tmp,tmp)
     342            /material->GetAtomicNumDensityVector()[1];
     343       
     344          if ( particleDefinition == instance->GetIon("alpha+") )
     345            crossSection = crossSection +  sigmaExcitation ;
     346         
     347          if ( particleDefinition == instance->GetIon("helium") )
     348            crossSection = crossSection + 2*sigmaExcitation ;
     349         
     350          delete excitationXS;
     351         
    296352      }     
     353*/     
    297354
    298355    }
    … …  
    327384  G4int level = RandomSelect(particleEnergy0,aDynamicParticle->GetDefinition());
    328385
    329   G4double excitationEnergy = waterExcitation.ExcitationEnergy(level);
     386  //  G4double excitationEnergy = waterExcitation.ExcitationEnergy(level);
     387
     388  // Dingfelder's excitation levels
     389  const G4double excitation[]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV};
     390  G4double excitationEnergy = excitation[level];
     391
    330392  G4double newEnergy = particleEnergy0 - excitationEnergy;
    331393 
    … …  
    364426  const G4double omegaj[]={0.85, 0.88, 0.88, 0.78, 0.78};
    365427 
     428  // Dingfelder's excitation levels
     429  const G4double Eliq[5]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV};
     430
    366431  G4int particleTypeIndex = 0;
    367432  G4DNAGenericIonsManager* instance;
    … …  
    369434
    370435  if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
     436  if (particleDefinition == instance->GetIon("hydrogen")) particleTypeIndex=0;
    371437  if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
    372438  if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
    … …  
    377443
    378444  // SI - added protection
    379   if (tCorrected < waterExcitation.ExcitationEnergy(excitationLevel)) return 0;
     445  if (tCorrected < Eliq[excitationLevel]) return 0;
    380446  //
    381447 
    … …  
    384450  G4double numerator;
    385451  numerator = std::pow(z * aj[excitationLevel], omegaj[excitationLevel]) *
    386     std::pow(tCorrected - waterExcitation.ExcitationEnergy(excitationLevel), nu);
     452    std::pow(tCorrected - Eliq[excitationLevel], nu);
     453
     454  // H case : see S. Uehara et al. IJRB 77, 2, 139-154 (2001) - section 3.3
     455 
     456  if (particleDefinition == instance->GetIon("hydrogen"))
     457     numerator = std::pow(z * 0.75*aj[excitationLevel], omegaj[excitationLevel]) *
     458     std::pow(tCorrected - Eliq[excitationLevel], nu);
     459
    387460
    388461  G4double power;
    … …  
    394467  G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
    395468
    396   zEff -= ( sCoefficient[0][particleTypeIndex] * S_1s(k, waterExcitation.ExcitationEnergy(excitationLevel), slaterEffectiveCharge[0][particleTypeIndex], 1.) +
    397             sCoefficient[1][particleTypeIndex] * S_2s(k, waterExcitation.ExcitationEnergy(excitationLevel), slaterEffectiveCharge[1][particleTypeIndex], 2.) +
    398             sCoefficient[2][particleTypeIndex] * S_2p(k, waterExcitation.ExcitationEnergy(excitationLevel), slaterEffectiveCharge[2][particleTypeIndex], 2.) );
     469  zEff -= ( sCoefficient[0][particleTypeIndex] * S_1s(k, Eliq[excitationLevel], slaterEffectiveCharge[0][particleTypeIndex], 1.) +
     470            sCoefficient[1][particleTypeIndex] * S_2s(k, Eliq[excitationLevel], slaterEffectiveCharge[1][particleTypeIndex], 2.) +
     471            sCoefficient[2][particleTypeIndex] * S_2p(k, Eliq[excitationLevel], slaterEffectiveCharge[2][particleTypeIndex], 2.) );
     472
     473  if (particleDefinition == instance->GetIon("hydrogen")) zEff = 1.;
    399474
    400475  G4double cross = sigma0 * zEff * zEff * numerator / denominator;
     476
    401477
    402478  return cross;
    … …  
    415491
    416492  if ( particle == instance->GetIon("alpha++") ||
    417        particle == G4Proton::ProtonDefinition()  )
     493       particle == G4Proton::ProtonDefinition()|| 
     494       particle == instance->GetIon("hydrogen")  ||
     495       particle == instance->GetIon("alpha+")  ||
     496       particle == instance->GetIon("helium")
     497     )
    418498  { 
    419499     while (i > 0)
    … …  
    437517  }
    438518
     519/*
    439520  // add ONE or TWO electron-water excitation for alpha+ and helium
    440521   
    … …  
    449530         
    450531          G4DNAEmfietzoglouExcitationModel * excitationXS = new G4DNAEmfietzoglouExcitationModel();
     532          excitationXS->Initialise(G4Electron::ElectronDefinition());
    451533         
    452534          G4double sigmaExcitation=0;
    … …  
    455537 
    456538          G4double partial = PartialCrossSection(k,i,particle);
     539
    457540          if (particle == instance->GetIon("alpha+")) partial = PartialCrossSection(k,i,particle) + sigmaExcitation;
    458541          if (particle == instance->GetIon("helium")) partial = PartialCrossSection(k,i,particle) + 2*sigmaExcitation;
     542
    459543          values.push_front(partial);
    460544          value += partial;
    … …  
    474558    }
    475559  }     
     560*/
    476561
    477562  return 0;
    … …  
    552637
    553638  G4double tElectron = 0.511/3728. * t;
    554   G4double value = 2. * tElectron * slaterEffectiveCharge / (energyTransferred * shellNumber);
    555  
     639 
     640  // The following is provided by M. Dingfelder
     641  G4double H = 2.*13.60569172 * eV;
     642  G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) *  (slaterEffectiveCharge/shellNumber);
     643
    556644  return value;
    557645}

  • trunk/source/processes/electromagnetic/lowenergy/src/G4DNARuddIonisationModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DNARuddIonisationModel.cc,v 1.17 2010/04/07 20:08:31 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4DNARuddIonisationModel.cc,v 1.21 2010/11/04 14:52:17 sincerti Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929
    … …  
    159159
    160160    lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
    161     highEnergyLimit[alphaPlusPlus] = 10. * MeV;
     161    highEnergyLimit[alphaPlusPlus] = 400. * MeV;
    162162
    163163    // Cross section
    … …  
    181181
    182182    lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
    183     highEnergyLimit[alphaPlus] = 10. * MeV;
     183    highEnergyLimit[alphaPlus] = 400. * MeV;
    184184
    185185    // Cross section
    … …  
    202202
    203203    lowEnergyLimit[helium] = lowEnergyLimitForZ2;
    204     highEnergyLimit[helium] = 10. * MeV;
     204    highEnergyLimit[helium] = 400. * MeV;
    205205
    206206    // Cross section
    … …  
    363363              sigma = table->FindValue(k);
    364364
    365               // BEGIN ELECTRON CORRECTION
    366               // add ONE or TWO electron-water excitation for alpha+ and helium
    367    
    368               if ( particleDefinition == instance->GetIon("alpha+")
    369                    ||
    370                    particleDefinition == instance->GetIon("helium")
    371                    )
    372               {
    373      
    374                   G4DNACrossSectionDataSet* electronDataset = new G4DNACrossSectionDataSet
    375                     (new G4LogLogInterpolation, eV, (1./3.343e22)*m*m);
    376        
    377                   electronDataset->LoadData("dna/sigma_ionisation_e_born");
    378 
    379                   G4double kElectron = k * 0.511/3728;
    380 
    381                   if ( particleDefinition == instance->GetIon("alpha+") )
    382                   {
    383                       G4double tmp1 = table->FindValue(k) + electronDataset->FindValue(kElectron);
    384                       delete electronDataset;
    385                       if (verboseLevel > 3)
    386                       {
    387                         G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
    388                         G4cout << " - Cross section per water molecule (cm^2)=" << tmp1/cm/cm << G4endl;
    389                         G4cout << " - Cross section per water molecule (cm^-1)=" <<
    390                         tmp1*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
    391                       }
    392                       return tmp1*material->GetAtomicNumDensityVector()[1];
    393                   }
    394 
    395                   if ( particleDefinition == instance->GetIon("helium") )
    396                   {
    397                       G4double tmp2 = table->FindValue(k) +  2. * electronDataset->FindValue(kElectron);
    398                       delete electronDataset;
    399                       if (verboseLevel > 3)
    400                       {
    401                         G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
    402                         G4cout << " - Cross section per water molecule (cm^2)=" << tmp2/cm/cm << G4endl;
    403                         G4cout << " - Cross section per water molecule (cm^-1)=" << tmp2*
    404                         material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
    405                       }
    406                       return tmp2*material->GetAtomicNumDensityVector()[1];
    407                   }
    408               }     
    409 
    410               // END ELECTRON CORRECTION
    411365         }
    412366      }
    … …  
    489443      G4ParticleDefinition* definition = particle->GetDefinition();
    490444      G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
     445      /*
    491446      G4double particleMass = definition->GetPDGMass();
    492447      G4double totalEnergy = k + particleMass;
    493448      G4double pSquare = k*(totalEnergy+particleMass);
    494449      G4double totalMomentum = std::sqrt(pSquare);
     450      */
    495451
    496452      G4int ionizationShell = RandomSelect(k,particleName);
    … …  
    511467      deltaDirection.rotateUz(primaryDirection);
    512468
     469      // Ignored for ions on electrons
     470      /*
    513471      G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
    514472
    … …  
    525483
    526484      fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
     485      */
     486      fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
     487
    527488      fParticleChangeForGamma->SetProposedKineticEnergy(k-bindingEnergy-secondaryKinetic);
    528489      fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
    … …  
    575536  }
    576537 
     538 
    577539  G4double secElecKinetic = 0.;
    578540 
    … …  
    616578 
    617579  phi = twopi * G4UniformRand();
    618   cosTheta = std::sqrt(secKinetic / maxSecKinetic);
     580 
     581  //cosTheta = std::sqrt(secKinetic / maxSecKinetic);
     582
     583  // Restriction below 100 eV from Emfietzoglou (2000)
     584
     585  if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
     586  else cosTheta = (2.*G4UniformRand())-1.;
     587 
    619588}
    620589
    … …  
    656625  G4double alphaConst ;
    657626
    658   const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
     627//  const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
     628// The following values are provided by M. dingfelder (priv. comm)
     629  const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
    659630
    660631  if (j == 4)
    … …  
    681652      E1 = 0.38;
    682653      A2 = 1.07;
    683       B2 = 14.6;
     654      // Value provided by M. Dingfelder (priv. comm)
     655      B2 = 11.6;
     656      //
    684657      C2 = 0.60;
    685658      D2 = 0.04;
    … …  
    697670
    698671  G4double w = wBig / Bj[ionizationLevelIndex];
     672  // Note that the following (j==4) cases are provided by   M. Dingfelder (priv. comm)
     673  if (j==4) w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
     674
    699675  G4double Ry = 13.6*eV;
    700676
    … …  
    713689      tau = (0.511/3728.) * k ;
    714690  }
    715  
     691
    716692  G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/Bj[ionizationLevelIndex]),2);
     693  if (j==4) S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/waterStructure.IonisationEnergy(ionizationLevelIndex)),2);
     694
    717695  G4double v2 = tau / Bj[ionizationLevelIndex];
     696  if (j==4) v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
     697
    718698  G4double v = std::sqrt(v2);
    719699  G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj[ionizationLevelIndex]));
     700  if (j==4) wc = 4.*v2 - 2.*v - (Ry/(4.*waterStructure.IonisationEnergy(ionizationLevelIndex)));
    720701
    721702  G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
    … …  
    731712    * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
    732713
     714  if (j==4) sigma = CorrectionFactor(particleDefinition, k)
     715    * Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
     716    * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
     717
    733718  if ( (particleDefinition == instance->GetIon("hydrogen")) && (ionizationLevelIndex==4))
    734719
    735     sigma = Gj[j] * (S/Bj[ionizationLevelIndex])
     720//    sigma = Gj[j] * (S/Bj[ionizationLevelIndex])
     721    sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
    736722    * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
    737723
    … …  
    756742  {
    757743      slaterEffectiveCharge[0]=2.0;
    758       slaterEffectiveCharge[1]=1.15;
    759       slaterEffectiveCharge[2]=1.15;
     744      // The following values are provided by M. Dingfelder (priv. comm)   
     745      slaterEffectiveCharge[1]=2.0;
     746      slaterEffectiveCharge[2]=2.0;
     747      //
    760748      sCoefficient[0]=0.7;
    761749      sCoefficient[1]=0.15;
    … …  
    780768      sigma = Gj[j] * (S/Bj[ionizationLevelIndex]) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
    781769   
     770      if (j==4) sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
     771                              * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
     772
    782773      G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
    783774 
    … …  
    852843
    853844  G4double tElectron = 0.511/3728. * t;
    854   G4double value = 2. * tElectron * slaterEffectiveChg / (energyTransferred * shellNumber);
     845  // The following values are provided by M. Dingfelder (priv. comm)   
     846  G4double H = 2.*13.60569172 * eV;
     847  G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) *  (slaterEffectiveChg/shellNumber);
    855848 
    856849  return value;
    … …  
    872865    {
    873866        G4double value = (std::log10(k/eV)-4.2)/0.5;
    874         return((0.8/(1+std::exp(value))) + 0.9);
     867        // The following values are provided by M. Dingfelder (priv. comm)   
     868        return((0.6/(1+std::exp(value))) + 0.9);
    875869    }
    876870    else
    … …  
    891885  G4DNAGenericIonsManager *instance;
    892886  instance = G4DNAGenericIonsManager::Instance();
    893   G4double kElectron(0);
    894 
    895   G4DNACrossSectionDataSet * electronDataset = new G4DNACrossSectionDataSet (new G4LogLogInterpolation, eV, (1./3.343e22)*m*m);
    896887 
    897   if ( particle == instance->GetIon("alpha+")->GetParticleName()
    898        ||
    899        particle == instance->GetIon("helium")->GetParticleName()
    900        )
    901   {     
    902       electronDataset->LoadData("dna/sigma_ionisation_e_born");
    903 
    904       kElectron = k * 0.511/3728;
    905        
    906   }     
    907  
    908   // END PART 1/2 OF ELECTRON CORRECTION
    909  
    910888  G4int level = 0;
    911889
    … …  
    931909              i--;
    932910              valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
    933 
    934               // BEGIN PART 2/2 OF ELECTRON CORRECTION
    935               // Use only electron partial cross sections
    936              
    937               if (particle == instance->GetIon("alpha+")->GetParticleName())
    938                 {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + electronDataset->GetComponent(i)->FindValue(kElectron); }
    939 
    940               if (particle == instance->GetIon("helium")->GetParticleName())
    941                 {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + 2*electronDataset->GetComponent(i)->FindValue(kElectron); }
    942 
    943               // BEGIN PART 2/2 OF ELECTRON CORRECTION
    944 
    945911              value += valuesBuffer[i];
    946912          }
    … …  
    958924              {
    959925                  delete[] valuesBuffer;
    960                  
    961                   if (electronDataset) delete electronDataset;
    962                  
    963926                  return i;
    964927              }
    … …  
    975938  }
    976939
    977   delete electronDataset;
    978      
    979940  return level;
    980941}

  • trunk/source/processes/electromagnetic/lowenergy/src/G4DNAScreenedRutherfordElasticModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DNAScreenedRutherfordElasticModel.cc,v 1.10 2010/01/07 18:10:50 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4DNAScreenedRutherfordElasticModel.cc,v 1.14 2010/09/08 14:07:16 sincerti Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929
    … …  
    4141{
    4242
    43   killBelowEnergy = 8.23*eV; // Minimum e- energy for energy loss by excitation
     43  killBelowEnergy = 0.025*eV; // Minimum e- energy for energy loss by excitation
    4444  lowEnergyLimit = 0 * eV;
    45   lowEnergyLimitOfModel = 7 * eV; // The model lower energy is 7 eV
    4645  intermediateEnergyLimit = 200 * eV; // Switch between two final state models
    47   highEnergyLimit = 10 * MeV;
     46  highEnergyLimit = 1. * MeV;
    4847  SetLowEnergyLimit(lowEnergyLimit);
    4948  SetHighEnergyLimit(highEnergyLimit);
    … …  
    174173  if (ekin < highEnergyLimit)
    175174  {
    176      
    177       //SI : XS must not be zero otherwise sampling of secondaries method ignored
    178       if (ekin < lowEnergyLimitOfModel) ekin = lowEnergyLimitOfModel;
    179       //
    180175     
    181176      G4double z = 10.;

  • trunk/source/processes/electromagnetic/lowenergy/src/G4Generator2BN.cc

    r819 r1340  
    2424// ********************************************************************
    2525//
     26// $Id: G4Generator2BN.cc,v 1.9 2010/10/14 14:01:02 vnivanch Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2628//
    2729// -------------------------------------------------------------------
    … …  
    4547// Class Description:
    4648//
    47 // Concrete base class for Bremsstrahlung Angular Distribution Generation - 2BN Distribution
     49// Concrete base class for Bremsstrahlung Angular Distribution Generation
     50// 2BN Distribution
    4851//
    4952// Class Description: End
    5053//
    5154// -------------------------------------------------------------------
    52 //
    53 //   
     55//   
    5456
    5557#include "G4Generator2BN.hh"
    … …  
    150152
    151153
    152 G4Generator2BN::G4Generator2BN(const G4String& name):G4VBremAngularDistribution(name)
     154G4Generator2BN::G4Generator2BN(const G4String&)
     155 : G4VBremAngularDistribution("AngularGen2BN")
    153156{
    154157  b = 1.2;
    … …  
    172175
    173176G4Generator2BN::~G4Generator2BN()
    174 {;}
     177{}
    175178
    176179//

  • trunk/source/processes/electromagnetic/lowenergy/src/G4Generator2BS.cc

    r819 r1340  
    2424// ********************************************************************
    2525//
     26// $Id: G4Generator2BS.cc,v 1.10 2010/10/14 14:01:02 vnivanch Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2628//
    2729// -------------------------------------------------------------------
    … …  
    3941//
    4042// Modifications:
    41 // 02 Jun 2003                    First implementation acording with new design
    42 // 05 Nov 2003  MGP               Fixed std namespace
    43 // 17 Nov 2003  MGP               Fixed compilation problem on Windows                 
     43// 02 Jun 2003               First implementation acording with new design
     44// 05 Nov 2003  MGP          Fixed std namespace
     45// 17 Nov 2003  MGP          Fixed compilation problem on Windows                 
     46// 12 Oct 2010  V.Ivanchenko Moved RejectionFunction inline, use G4Pow to speadup
    4447//
    4548// Class Description:
    4649//
    47 // Concrete base class for Bremsstrahlung Angular Distribution Generation - 2BS Distribution
     50// Concrete base class for Bremsstrahlung Angular Distribution Generation
     51// 2BS Distribution
    4852//
    4953// Class Description: End
    … …  
    5155// -------------------------------------------------------------------
    5256//
    53 //   
    5457
    5558#include "G4Generator2BS.hh"
    56 #include "Randomize.hh"
    57 //    
     59#include "Randomize.hh"   
     60#include "G4Pow.hh"   
    5861
    59 G4Generator2BS::G4Generator2BS(const G4String& name):G4VBremAngularDistribution(name)
    60 {;}
     62//
     63
     64G4Generator2BS::G4Generator2BS(const G4String&)
     65  : G4VBremAngularDistribution("AngularGen2BS")
     66{
     67  g4pow = G4Pow::GetInstance();
     68}
    6169
    6270//   
    6371
    6472G4Generator2BS::~G4Generator2BS()
    65 {;}
     73{}
    6674
    6775//
    … …  
    7987  // Departement of Medical Physics, Memorial Sloan-Kettering Cancer Center, New York
    8088
    81 
    8289  G4double theta = 0;
    8390
    … …  
    8794  G4double gMaxEnergy = (pi*initialTotalEnergy)*(pi*initialTotalEnergy);
    8895
    89   G4double Zeff = std::sqrt(static_cast<G4double>(Z) * (static_cast<G4double>(Z) + 1.0));
    90   z = (0.00008116224*(std::pow(Zeff,0.3333333)));
     96  //G4double Zeff = std::sqrt(static_cast<G4double>(Z) * (static_cast<G4double>(Z) + 1.0));
     97  //z = (0.00008116224*(std::pow(Zeff,0.3333333)));
     98
     99  // VI speadup
     100  z = 0.00008116224*(g4pow->Z13(Z) + g4pow->Z13(Z+1));
    91101
    92102  // Rejection arguments
    … …  
    97107
    98108  // Calculate rejection function at 0, 1 and Emax
    99   G4double gfunction0 = RejectionFunction(0);
    100   G4double gfunction1 = RejectionFunction(1);
     109  G4double gfunction0 = RejectionFunction(0.0);
     110  G4double gfunction1 = RejectionFunction(1.0);
    101111  G4double gfunctionEmax = RejectionFunction(gMaxEnergy);
    102 
    103112
    104113  // Calculate Maximum value
    … …  
    110119  do{
    111120    rand = G4UniformRand();
    112     rand = rand/(1-rand+1.0/gMaxEnergy);
     121    rand /= (1 - rand + 1.0/gMaxEnergy);
    113122    gfunctionTest = RejectionFunction(rand);
    114123    randTest = G4UniformRand();
    115124
    116   }while(randTest > (gfunctionTest/gMaximum));
     125  } while(randTest*gMaximum > gfunctionTest);
    117126
    118127  theta = std::sqrt(rand)/initialTotalEnergy;
    119128
    120 
    121129  return theta;
    122130}
     131
    123132//
    124 
    125 G4double G4Generator2BS::RejectionFunction(G4double value) const
    126 {
    127 
    128   G4double argument = (1+value)*(1+value);
    129 
    130   G4double gfunction = (4+std::log(rejection_argument3+(z/argument)))*
    131     ((4*EnergyRatio*value/argument)-rejection_argument1)+rejection_argument2;
    132 
    133   return gfunction;
    134 
    135 }
    136133
    137134void G4Generator2BS::PrintGeneratorInformation() const
    138135{
    139136  G4cout << "\n" << G4endl;
    140   G4cout << "Bremsstrahlung Angular Generator is 2BS Generator from 2BS Koch & Motz distribution (Rev Mod Phys 31(4), 920 (1959))" << G4endl;
     137  G4cout << "Bremsstrahlung Angular Generator is 2BS Generator "
     138         << "from 2BS Koch & Motz distribution (Rev Mod Phys 31(4), 920 (1959))" << G4endl;
    141139  G4cout << "Sampling algorithm adapted from PIRS-0203" << G4endl;
    142140  G4cout << "\n" << G4endl;

  • trunk/source/processes/electromagnetic/lowenergy/src/G4IonParametrisedLossModel.cc

    r1228 r1340  
    2424// ********************************************************************
    2525//
    26 //
     26// $Id: G4IonParametrisedLossModel.cc,v 1.10 2010/11/04 12:21:48 vnivanch Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2728//
    2829// ===========================================================================
    … …  
    6465//                               functions accordingly, added new cache param.)
    6566//                             - Removed GetRange function (AL) 
     67//                04. 11. 2010 - Moved virtual methods to the source (VI)
    6668//
    6769//
    … …  
    107109    nmbSubBins(100),
    108110    particleChangeLoss(0),
    109     modelIsInitialised(false),
    110     corrections(0),
    111111    corrFactor(1.0),
    112112    energyLossLimit(0.01),
    113     cutEnergies(0) {
    114 
     113    cutEnergies(0)
     114{
    115115  genericIon = G4GenericIon::Definition();
    116116  genericIonPDGMass = genericIon -> GetPDGMass();
     117  corrections = G4LossTableManager::Instance() -> EmCorrections();
    117118 
    118119  // The upper limit of the current model is set to 100 TeV
    … …  
    192193  return couple -> GetMaterial() -> GetIonisation() ->
    193194                                                  GetMeanExcitationEnergy();
     195}
     196
     197// #########################################################################
     198
     199G4double G4IonParametrisedLossModel::MaxSecondaryEnergy(
     200                             const G4ParticleDefinition* particle,
     201                             G4double kineticEnergy) {
     202
     203  // ############## Maximum energy of secondaries ##########################
     204  // Function computes maximum energy of secondary electrons which are
     205  // released by an ion
     206  //
     207  // See Geant4 physics reference manual (version 9.1), section 9.1.1
     208  //
     209  // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
     210  //       C.Caso et al. (Part. Data Group), Europ. Phys. Jour. C 3 1 (1998).
     211  //       B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
     212  //
     213  // (Implementation adapted from G4BraggIonModel)
     214
     215  if(particle != cacheParticle) UpdateCache(particle);
     216
     217  G4double tau  = kineticEnergy/cacheMass;
     218  G4double tmax = 2.0 * electron_mass_c2 * tau * (tau + 2.) /
     219                  (1. + 2.0 * (tau + 1.) * cacheElecMassRatio +
     220                  cacheElecMassRatio * cacheElecMassRatio);
     221
     222  return tmax;
     223}
     224
     225// #########################################################################
     226
     227G4double G4IonParametrisedLossModel::GetChargeSquareRatio(
     228                             const G4ParticleDefinition* particle,
     229                             const G4Material* material,
     230                             G4double kineticEnergy) {    // Kinetic energy
     231
     232  G4double chargeSquareRatio = corrections ->
     233                                     EffectiveChargeSquareRatio(particle,
     234                                                                material,
     235                                                                kineticEnergy);
     236  corrFactor = chargeSquareRatio *
     237                       corrections -> EffectiveChargeCorrection(particle,
     238                                                                material,
     239                                                                kineticEnergy);
     240  return corrFactor;
     241}
     242
     243// #########################################################################
     244
     245G4double G4IonParametrisedLossModel::GetParticleCharge(
     246                             const G4ParticleDefinition* particle,
     247                             const G4Material* material,
     248                             G4double kineticEnergy) {   // Kinetic energy
     249
     250  return corrections -> GetParticleCharge(particle, material, kineticEnergy);
    194251}
    195252
    … …  
    295352  }
    296353
    297   // The particle change object is cast to G4ParticleChangeForLoss
    298   if(! modelIsInitialised) {
    299 
    300      modelIsInitialised = true;
    301      corrections = G4LossTableManager::Instance() -> EmCorrections();
    302 
    303      if(!particleChangeLoss) {
    304         if(pParticleChange) {
    305 
    306            particleChangeLoss = reinterpret_cast<G4ParticleChangeForLoss*>
    307                (pParticleChange);
    308         }
    309         else {
    310           particleChangeLoss = new G4ParticleChangeForLoss();
    311         }
    312      }
     354  // The particle change object
     355  if(! particleChangeLoss) {
     356    particleChangeLoss = GetParticleChangeForLoss();
    313357  }
    314358 

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreIonisationModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LivermoreIonisationModel.cc,v 1.8 2010/03/26 09:32:50 pandola Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4LivermoreIonisationModel.cc,v 1.9 2010/10/13 07:28:47 pandola Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929// Author: Luciano Pandola
    … …  
    4444//                  - atomic deexcitation managed via G4VEmModel::DeexcitationFlag() is
    4545//                    set as "true" (default would be false)
     46// 12 Oct 2010   L Pandola
     47//                  - add debugging information about energy in
     48//                    SampleDeexcitationAlongStep()
     49//                  - generate fluorescence SampleDeexcitationAlongStep() only above
     50//                    the cuts.
     51//
    4652//
    4753
    … …  
    125131    }
    126132  energySpectrum = new G4eIonisationSpectrum();
    127   if (verboseLevel > 0)
     133  if (verboseLevel > 3)
    128134    G4cout << "G4VEnergySpectrum is initialized" << G4endl;
    129135
    … …  
    159165    }
    160166
    161   if (verboseLevel > 1)
     167  if (verboseLevel > 3)
    162168    {
    163169      G4cout << "Cross section data: " << G4endl;
    … …  
    440446  //secondaries. The eloss value is updated.
    441447  G4double energyLossBefore = eloss;
     448
    442449  if (verboseLevel > 2)
    443     G4cout << "Energy loss along step before deexcitation : " << energyLossBefore/keV <<
    444       " keV" << G4endl;
    445 
     450    {
     451      G4cout << "-----------------------------------------------------------" << G4endl;
     452      G4cout << " SampleDeexcitationAlongStep() from G4LivermoreIonisation" << G4endl;
     453      G4cout << "Energy loss along step before deexcitation : " << energyLossBefore/keV <<
     454        " keV" << G4endl;
     455    }
    446456  G4double incidentEnergy = theTrack.GetDynamicParticle()->GetKineticEnergy();
    447457
    … …  
    453463  G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
    454464 
    455   //Notice: in LowEnergyIonisation, fluorescence is always generated above 250 eV
    456   //not above the tracking cut.
    457   //G4double cutForLowEnergySecondaryParticles = 250.0*eV;
    458465
    459466  std::vector<G4DynamicParticle*>* deexcitationProducts =
    … …  
    505512                              {
    506513                                e = aSecondary->GetKineticEnergy();
    507                                 if ( eTot + e <= eloss )
     514                                G4double itsCut = cutg;
     515                                if (aSecondary->GetParticleDefinition() == G4Electron::Electron())
     516                                  itsCut = cute;
     517                                if ( eTot + e <= eloss && e > itsCut )
    508518                                  {
    509519                                    eTot += e;
    … …  
    524534    }
    525535
     536  G4double energyLossInFluorescence = 0.0;
    526537  size_t nSecondaries = deexcitationProducts->size();
    527538  if (nSecondaries > 0)
    528539    {
    529       fParticleChange->SetNumberOfSecondaries(nSecondaries);
     540      //You may have already secondaries produced by SampleSubCutSecondaries()
     541      //at the process G4VEnergyLossProcess
     542      G4int secondariesBefore = fParticleChange->GetNumberOfSecondaries();
     543      fParticleChange->SetNumberOfSecondaries(nSecondaries+secondariesBefore);
    530544      const G4StepPoint* preStep = theTrack.GetStep()->GetPreStepPoint();
    531545      const G4StepPoint* postStep = theTrack.GetStep()->GetPostStepPoint();
    … …  
    538552      G4double time, q;
    539553      G4ThreeVector position;
    540      
     554
    541555      for (size_t i=0; i<nSecondaries; i++)
    542556        {
    … …  
    553567                  position += r;
    554568                  G4Track* newTrack = new G4Track(part, time, position);
     569                  energyLossInFluorescence += eSecondary;
    555570                  pParticleChange->AddSecondary(newTrack);
    556571                }
    … …  
    566581  delete deexcitationProducts;
    567582 
     583  //Check and verbosities. Ensure energy conservation
    568584  if (verboseLevel > 2)
    569     G4cout << "Energy loss along step after deexcitation : " << eloss/keV << 
    570       " keV" << G4endl;
     585    {
     586      G4cout << "Energy loss along step after deexcitation : " << eloss/keV << 
     587        " keV" << G4endl;
     588    }
     589  if (verboseLevel > 1)
     590    {
     591      G4cout << "------------------------------------------------------------------" << G4endl;
     592      G4cout << "Energy in fluorescence: " << energyLossInFluorescence/keV << " keV" << G4endl;
     593      G4cout << "Residual energy loss: " << eloss/keV << " keV " << G4endl;
     594      G4cout << "Total final: " << (energyLossInFluorescence+eloss)/keV << " keV" << G4endl;
     595      G4cout << "Total initial: " << energyLossBefore/keV << " keV" << G4endl; 
     596      G4cout << "------------------------------------------------------------------" << G4endl;
     597    }
     598  if (verboseLevel > 0)
     599    {
     600      if (std::fabs(energyLossBefore-energyLossInFluorescence-eloss)>10*eV)
     601        {
     602          G4cout << "Found energy non-conservation at SampleDeexcitationAlongStep() " << G4endl;
     603          G4cout << "Energy in fluorescence: " << energyLossInFluorescence/keV << " keV" << G4endl;
     604          G4cout << "Residual energy loss: " << eloss/keV << " keV " << G4endl;
     605          G4cout << "Total final: " << (energyLossInFluorescence+eloss)/keV << " keV" << G4endl;
     606          G4cout << "Total initial: " << energyLossBefore/keV << " keV" << G4endl;     
     607        }
     608    }
    571609}
    572610
    … …  
    655693           xsis->AddComponent(p);
    656694         }
    657        if(verboseLevel>0) xsis->PrintData();
     695       if(verboseLevel>3) xsis->PrintData();
    658696       shellVacancy->AddXsiTable(xsis);
    659697    }

  • trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePhotoElectricModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4LivermorePhotoElectricModel.cc,v 1.11 2010/03/26 09:32:50 pandola Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4LivermorePhotoElectricModel.cc,v 1.12 2010/10/13 07:15:42 pandola Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929//
    … …  
    304304  if (!DeexcitationFlag() && augerbool)
    305305    {
    306       G4cout << "WARNING - G4PenelopePhotoElectricModel" << G4endl;
     306      G4cout << "WARNING - G4LivermorePhotoElectricModel" << G4endl;
    307307      G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
    308308      G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;

  • trunk/source/processes/electromagnetic/lowenergy/src/G4Penelope08ComptonModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4Penelope08ComptonModel.cc,v 1.6 2010/04/12 13:53:29 pandola Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4Penelope08ComptonModel.cc,v 1.7 2010/07/28 07:09:16 pandola Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929// Author: Luciano Pandola
    … …  
    4242#include "G4VEMDataSet.hh"
    4343#include "G4PhysicsTable.hh"
    44 #include "G4ElementTable.hh"
    45 #include "G4Element.hh"
    4644#include "G4PhysicsLogVector.hh"
    4745#include "G4AtomicTransitionManager.hh"

  • trunk/source/processes/electromagnetic/lowenergy/src/G4Penelope08GammaConversionModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4Penelope08GammaConversionModel.cc,v 1.3 2010/06/25 09:41:17 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4Penelope08GammaConversionModel.cc,v 1.4 2010/07/28 07:09:16 pandola Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929// Author: Luciano Pandola
    … …  
    8787  if (fScreeningFunction)
    8888    delete fScreeningFunction;
    89 
    9089}
    9190

  • trunk/source/processes/electromagnetic/lowenergy/src/G4Penelope08PhotoElectricModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4Penelope08PhotoElectricModel.cc,v 1.4 2010/06/25 09:41:19 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4Penelope08PhotoElectricModel.cc,v 1.5 2010/07/28 07:09:16 pandola Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929// Author: Luciano Pandola
    … …  
    4242#include "G4ElementTable.hh"
    4343#include "G4Element.hh"
    44 #include "G4CrossSectionHandler.hh"
    4544#include "G4AtomicTransitionManager.hh"
    4645#include "G4AtomicShell.hh"
    … …  
    487486  //(theTable)[ishell] --> cross section for shell (ishell-1)
    488487
    489 
    490 
    491488  //reserve space for the vectors
    492489  //everything is log-log

  • trunk/source/processes/electromagnetic/lowenergy/src/G4Penelope08RayleighModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4Penelope08RayleighModel.cc,v 1.2 2010/06/25 09:41:28 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4Penelope08RayleighModel.cc,v 1.3 2010/07/28 07:09:16 pandola Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929// Author: Luciano Pandola
    … …  
    4343#include "G4ElementTable.hh"
    4444#include "G4Element.hh"
    45 #include "G4PenelopeIntegrator.hh"
    4645#include "G4PhysicsFreeVector.hh"
    4746
    … …  
    224223
    225224   if (verboseLevel > 2)
    226     G4cout << "Compton cross Section at " << energy/keV << " keV for Z=" << Z <<
     225    G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z=" << Z <<
    227226      " = " << cross/barn << " barn" << G4endl;
    228227    return cross;
    … …  
    414413  fParticleChange->SetProposedKineticEnergy(photonEnergy0) ;
    415414 
    416 
    417415  return;
    418416}

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeOscillatorManager.cc

    r1337 r1340  
    6666{
    6767  Clear();
     68  delete instance;
    6869}
    6970 

  • trunk/source/processes/electromagnetic/lowenergy/src/G4PhotoElectricAngularGeneratorSauterGavrila.cc

    r1337 r1340  
    7979
    8080  if (gamma > 5.) {
    81     return G4ThreeVector(sinteta*cosphi, sinteta*sinphi, costeta);
     81    G4ThreeVector direction (sinteta*cosphi, sinteta*sinphi, costeta);
     82    return direction;
     83    // Bugzilla 1120
     84    // SI on 05/09/2010 as suggested by JG 04/09/10
    8285  }
    8386

  • trunk/source/processes/electromagnetic/standard/History

    r1315 r1340  
    1 $Id: History,v 1.504 2010/06/04 09:25:12 vnivanch Exp $
     1$Id: History,v 1.514 2010/11/04 17:31:43 vnivanch Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     2004 November 2010: V.Ivant (emstand-V09-03-25)
     21- Fixed warnings of the Coverity tool (pedantic initialisation);
     22
     2326 October 2010: V.Ivant (emstand-V09-03-24)
     24- Fixed problem reported by the Coverity tool (mainly pedantic
     25  initialisation);
     26- G4BetheHeitlerModel, G4PairProductionRelModel - removed unused
     27  internal table of cross section
     28- G4WaterStopping - fixed bug in index of data for Fe ion
     29
     3014 October 2010: V.Ivant (emstand-V09-03-23)
     31- G4ModifiedTsai - added (moved from lowenergy)
     32- G4PairProductionRelModel, G4eBremsstrahlungRelModel - return back
     33     the version of 9.4beta (disable tag 21); use general interface
     34     to sample polar angular distribution (G4ModifiedTsai - default)
     35
     3628 September 2010: V.Ivant (emstand-V09-03-22)
     37- G4ionIonisation, G4alphaIonisation - removed obsolete nuclear stopping
     38      flag and Get/Set methods (nuclear stopping is simulated by the
     39      dedicated G4NuclearStopping process)
     40
     4112 September 2010: V.Ivant (emstand-V09-03-21)
     42- G4BetheBlochModel - fixed shell corrections
     43- G4PairProductionRelModel, G4eBremsstrahlungRelModel - fixed
     44     SetCurrentElement method, fix is important for unit tests
     45
     4612 August 2010: V.Ivant (emstand-V09-03-20)
     47- G4UniversalFluctuation - L.Urban revision of width correction
     48  providing better results for thin targets and a good tail 
     49
     5014 June 2010:  V.Ivant (emstand-V09-03-19)
     51- G4UniversalFluctuation93 new class keeping version of the release 9.3
     52  for the G4UniversalFluctuation
     53- G4UniversalFluctuation - improved comments
     54
     5511 June 2010:  V.Ivant (emstand-V09-03-18)
     56- G4UniversalFluctuation - L.Urban fixed anomaly at the tail of
     57                           distribution for thin targets
    1958
    20594 June 2010:  V.Ivant (emstand-V09-03-17)

  • trunk/source/processes/electromagnetic/standard/include/G4BetheHeitlerModel.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BetheHeitlerModel.hh,v 1.6 2007/05/22 17:34:36 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4BetheHeitlerModel.hh,v 1.9 2010/10/26 10:35:22 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6363
    6464  G4BetheHeitlerModel(const G4ParticleDefinition* p = 0,
    65                       const G4String& nam = "Bethe-Heitler");
     65                      const G4String& nam = "BetheHeitler");
    6666 
    6767  virtual ~G4BetheHeitlerModel();
    … …  
    9797  G4ParticleDefinition*     thePositron;
    9898  G4ParticleChangeForGamma* fParticleChange;
    99   G4PhysicsTable*           theCrossSectionTable;
    100 
    101   G4double                  lowGammaEnergy;
    102   G4double                  highGammaEnergy;
    103 
    104   G4int                     nbins;
    105   size_t                    indexZ[120];
    106  
    10799};
    108100

  • trunk/source/processes/electromagnetic/standard/include/G4CoulombScattering.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4CoulombScattering.hh,v 1.14 2010/02/17 18:59:22 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4CoulombScattering.hh,v 1.15 2010/10/25 19:13:23 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9393  G4double thEnergyElec;
    9494  G4bool isInitialised;
    95   G4bool buildElmTableFlag;
    9695  const G4ParticleDefinition* aParticle;
    9796

  • trunk/source/processes/electromagnetic/standard/include/G4PAIySection.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4PAIySection.hh,v 1.1 2007/10/01 17:45:14 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4PAIySection.hh,v 1.2 2010/11/04 17:30:31 vnivanch Exp $
     28// GEANT4 tag $Name: emstand-V09-03-25 $
    2929//
    3030//
    … …  
    110110  // Inline access functions
    111111
    112   G4int GetNumberOfGammas() const { return fNumberOfGammas ; }
    113          
    114   G4int GetSplineSize() const { return fSplineNumber ; }
    115          
    116   G4int GetIntervalNumber() const { return fIntervalNumber ; }
    117 
    118   G4double GetEnergyInterval(G4int i){ return fEnergyInterval[i] ; }
    119 
    120   G4double GetDifPAIySection(G4int i){ return fDifPAIySection[i] ; }
    121   G4double GetPAIdNdxCrenkov(G4int i){ return fdNdxCerenkov[i] ; }
    122   G4double GetPAIdNdxPlasmon(G4int i){ return fdNdxPlasmon[i] ; }
    123          
    124   G4double GetMeanEnergyLoss() const {return fIntegralPAIySection[0] ; }
    125   G4double GetMeanCerenkovLoss() const {return fIntegralCerenkov[0] ; }
    126   G4double GetMeanPlasmonLoss() const {return fIntegralPlasmon[0] ; }
    127 
    128   G4double GetNormalizationCof() const { return fNormalizationCof ; }
     112  inline G4int GetNumberOfGammas() const { return fNumberOfGammas ; }
     113         
     114  inline G4int GetSplineSize() const { return fSplineNumber ; }
     115         
     116  inline G4int GetIntervalNumber() const { return fIntervalNumber ; }
     117
     118  inline G4double GetEnergyInterval(G4int i){ return fEnergyInterval[i] ; }
     119
     120  inline G4double GetDifPAIySection(G4int i){ return fDifPAIySection[i] ; }
     121  inline G4double GetPAIdNdxCrenkov(G4int i){ return fdNdxCerenkov[i] ; }
     122  inline G4double GetPAIdNdxPlasmon(G4int i){ return fdNdxPlasmon[i] ; }
     123         
     124  inline G4double GetMeanEnergyLoss() const {return fIntegralPAIySection[0] ; }
     125  inline G4double GetMeanCerenkovLoss() const {return fIntegralCerenkov[0] ; }
     126  inline G4double GetMeanPlasmonLoss() const {return fIntegralPlasmon[0] ; }
     127
     128  inline G4double GetNormalizationCof() const { return fNormalizationCof ; }
    129129         
    130130  inline G4double GetPAItable(G4int i,G4int j) const ;
    131131
    132   inline G4double    GetLorentzFactor(G4int i) const ;
     132  inline G4double GetLorentzFactor(G4int i) const ;
    133133                 
    134134  inline G4double GetSplineEnergy(G4int i) const ;
    … …  
    185185  G4double   fPAItable[500][112] ;         // Output array
    186186
    187 } ;   
     187};
    188188
    189189////////////////  Inline methods //////////////////////////////////

  • trunk/source/processes/electromagnetic/standard/include/G4PairProductionRelModel.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PairProductionRelModel.hh,v 1.3 2009/06/04 13:45:53 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4PairProductionRelModel.hh,v 1.9 2010/10/26 10:35:22 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6464
    6565  G4PairProductionRelModel(const G4ParticleDefinition* p = 0,
    66                       const G4String& nam = "Bethe-Heitler");
     66                           const G4String& nam = "BetheHeitlerLPM");
    6767 
    6868  virtual ~G4PairProductionRelModel();
    … …  
    8888
    8989  // * fast inline functions *
    90   inline void SetCurrentElement(const G4double);
     90  inline void SetCurrentElement(G4double /*Z*/);
    9191
    9292  // set / get methods
    … …  
    110110  G4double ScreenFunction2(G4double ScreenVariable);
    111111
    112 
    113 
    114 
    115112  G4double ComputeXSectionPerAtom(G4double totalEnergy, G4double Z);
    116113
    … …  
    128125  G4ParticleDefinition*     thePositron;
    129126  G4ParticleChangeForGamma* fParticleChange;
    130   G4PhysicsTable*           theCrossSectionTable;
    131 
    132   G4double                  lowGammaEnergy;
    133   G4double                  highGammaEnergy;
    134 
    135   G4int                     nbins;
    136   size_t                    indexZ[120];
    137127
    138128  G4double fLPMconstant;
    … …  
    196186//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    197187
    198 inline void G4PairProductionRelModel::SetCurrentElement(const G4double Z)
     188inline void G4PairProductionRelModel::SetCurrentElement(G4double Z)
    199189{
    200190  if(Z != currentZ) {
    … …  
    214204      Finel = facFinel - 2.*lnZ/3. ;
    215205    }
    216 
    217206    fCoulomb=GetCurrentElement()->GetfCoulomb();
    218207  }
    219208}
     209
    220210//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    221211

  • trunk/source/processes/electromagnetic/standard/include/G4UniversalFluctuation.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4UniversalFluctuation.hh,v 1.10 2009/03/19 14:15:17 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4UniversalFluctuation.hh,v 1.12 2010/08/08 08:19:59 urban Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    3434// File name:     G4UniversalFluctuation
    3535//
    36 // Author:        Vladimir Ivanchenko
     36// Author:        V.Ivanchenko make a class with the Laszlo Urban model
    3737//
    3838// Creation date: 03.01.2002
    … …  
    110110  G4double ipotLogFluct;
    111111  G4double e0;
     112  G4double esmall;
    112113
    113114  G4double e1,e2;
    … …  
    116117  G4double theBohrBeta2;
    117118  G4double minLoss;
    118   G4double nmaxCont1;
    119   G4double nmaxCont2;
     119  G4double nmaxCont;
     120  G4double rate,fw;
     121
    120122
    121123};

  • trunk/source/processes/electromagnetic/standard/include/G4alphaIonisation.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4alphaIonisation.hh,v 1.1 2009/11/10 11:50:30 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4alphaIonisation.hh,v 1.3 2010/10/26 10:06:12 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    7171  virtual void PrintInfo();
    7272
    73   void ActivateNuclearStopping(G4bool);
    74 
    7573protected:
    7674
    … …  
    8280
    8381  inline G4double BetheBlochEnergyThreshold();
    84 
    85   inline G4bool NuclearStoppingFlag();
    8682
    8783private:
    … …  
    9894
    9995  G4bool     isInitialised;
    100   G4bool     nuclearStopping;
    10196};
    10297
    10398//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    105 
    106 inline void G4alphaIonisation::ActivateNuclearStopping(G4bool val)
    107 {
    108   nuclearStopping = val;
    109 }
    110 
    11199//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    112100
    … …  
    118106//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    119107
    120 inline G4bool G4alphaIonisation::NuclearStoppingFlag()
    121 {
    122   return nuclearStopping;
    123 }
    124 
    125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    126 
    127108#endif

  • trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlungRelModel.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4eBremsstrahlungRelModel.hh,v 1.11 2009/02/20 12:06:37 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4eBremsstrahlungRelModel.hh,v 1.14 2010/10/26 10:35:22 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6565
    6666  G4eBremsstrahlungRelModel(const G4ParticleDefinition* p = 0,
    67                             const G4String& nam = "eBremRel");
     67                            const G4String& nam = "eBremLPM");
    6868
    6969  virtual ~G4eBremsstrahlungRelModel();
    … …  
    159159
    160160private:
     161
    161162  // consts
    162   G4double highKinEnergy;
    163163  G4double lowKinEnergy;
    164164  G4double fMigdalConstant;
    … …  
    195195    }
    196196
    197     fCoulomb=GetCurrentElement()->GetfCoulomb();
    198     fMax   = Fel-fCoulomb + Finel/currentZ  +  (1.+1./currentZ)/12.;
     197    fCoulomb = GetCurrentElement()->GetfCoulomb();
     198    fMax = Fel-fCoulomb + Finel/currentZ  +  (1.+1./currentZ)/12.;
    199199  }
    200200}

  • trunk/source/processes/electromagnetic/standard/include/G4eMultipleScattering.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4eMultipleScattering.hh,v 1.3 2009/11/01 13:04:12 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4eMultipleScattering.hh,v 1.4 2010/10/26 10:39:02 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -----------------------------------------------------------------------------
    … …  
    6464public:    // with description
    6565
    66   G4eMultipleScattering(const G4String& processName="msc");
     66  G4eMultipleScattering(const G4String& processName = "msc");
    6767
    6868  virtual ~G4eMultipleScattering();
    … …  
    7575
    7676  // geom. step length distribution should be sampled or not
    77   void Setsamplez(G4bool value) { samplez = value;};
     77  //void Setsamplez(G4bool value) { samplez = value;};
    7878
    7979  // to reduce the energy/step dependence
    80   void Setdtrl(G4double value) { dtrl = value;};
     80  //void Setdtrl(G4double value) { dtrl = value;};
    8181
    8282  // 'soften' step limitation above lambdalimit
    83   void SetLambdalimit(G4double value) { lambdalimit = value;};
     83  //void SetLambdalimit(G4double value) { lambdalimit = value;};
    8484
    8585protected:
    … …  
    9090private:        // data members
    9191
    92   G4double lambdalimit;
    93   G4double dtrl;
     92  //  G4double lambdalimit;
     93  //G4double dtrl;
    9494
    95   G4bool   samplez;
     95  //G4bool   samplez;
    9696  G4bool   isInitialized;
    9797

  • trunk/source/processes/electromagnetic/standard/include/G4ionIonisation.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4ionIonisation.hh,v 1.57 2009/02/20 12:06:37 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4ionIonisation.hh,v 1.58 2010/09/28 15:50:00 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9595  void ActivateStoppingData(G4bool);
    9696
    97   void ActivateNuclearStopping(G4bool);
    98 
    9997protected:
    10098
    … …  
    106104
    107105  inline G4double BetheBlochEnergyThreshold();
    108 
    109   inline G4bool NuclearStoppingFlag();
    110106
    111107private:
    … …  
    123119  G4bool     isInitialised;
    124120  G4bool     stopDataActive;
    125   G4bool     nuclearStopping;
    126121};
    127122
    … …  
    136131//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    137132
    138 inline void G4ionIonisation::ActivateNuclearStopping(G4bool val)
    139 {
    140   nuclearStopping = val;
    141 }
    142 
    143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    144 
    145133inline G4double G4ionIonisation::BetheBlochEnergyThreshold()
    146134{
    … …  
    150138//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    151139
    152 inline G4bool G4ionIonisation::NuclearStoppingFlag()
    153 {
    154   return nuclearStopping;
    155 }
    156 
    157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    158 
    159140#endif

  • trunk/source/processes/electromagnetic/standard/src/G4BetheBlochModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BetheBlochModel.cc,v 1.37 2010/05/27 10:25:59 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4BetheBlochModel.cc,v 1.40 2010/11/04 17:30:31 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-25 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8282{
    8383  fParticleChange = 0;
     84  theElectron = G4Electron::Electron();
    8485  if(p) {
    8586    SetGenericIon(p);
    8687    SetParticle(p);
    87   }
    88   theElectron = G4Electron::Electron();
     88  } else {
     89    SetParticle(theElectron);
     90  }
    8991  corr = G4LossTableManager::Instance()->EmCorrections(); 
    9092  nist = G4NistManager::Instance();
    … …  
    117119  //     << G4endl;
    118120
    119   corrFactor = chargeSquare;
    120121  // always false before the run
    121122  SetDeexcitationFlag(false);
    … …  
    157158  G4double q = particle->GetPDGCharge()/eplus;
    158159  chargeSquare = q*q;
     160  corrFactor = chargeSquare;
    159161  ratio = electron_mass_c2/mass;
    160162  G4double magmom =
    … …  
    275277
    276278  // shell correction
    277   dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
     279  //dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
     280  dedx -= corr->ShellCorrection(p,material,kineticEnergy);
    278281
    279282  // now compute the total ionization loss
    280 
    281   if (dedx < 0.0) dedx = 0.0 ;
    282 
    283283  dedx *= twopi_mc2_rcl2*chargeSquare*eDensity/beta2;
    284284
    … …  
    289289    dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
    290290  }
     291
     292  if (dedx < 0.0) { dedx = 0.0; }
    291293  return dedx;
    292294}

  • trunk/source/processes/electromagnetic/standard/src/G4BetheHeitlerModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BetheHeitlerModel.cc,v 1.13 2009/04/09 18:41:18 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4BetheHeitlerModel.cc,v 1.15 2010/10/25 19:02:32 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4343// 16-11-05 replace shootBit() by G4UniformRand()  mma
    4444// 04-12-05 SetProposedKineticEnergy(0.) for the killed photon (mma)
    45 // 20-02-20 SelectRandomElement is called for any initial gamma energy
     45// 20-02-07 SelectRandomElement is called for any initial gamma energy
    4646//          in order to have selected element for polarized model (VI)
     47// 25-10-10 Removed unused table, added element selector (VI)
    4748//
    4849// Class Description:
    … …  
    5859#include "G4Gamma.hh"
    5960#include "Randomize.hh"
    60 #include "G4DataVector.hh"
    61 #include "G4PhysicsLogVector.hh"
    6261#include "G4ParticleChangeForGamma.hh"
    63 #include "G4LossTableManager.hh"
    6462
    6563//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    6967G4BetheHeitlerModel::G4BetheHeitlerModel(const G4ParticleDefinition*,
    7068                                         const G4String& nam)
    71   : G4VEmModel(nam),
    72     theCrossSectionTable(0),
    73     nbins(10)
     69  : G4VEmModel(nam)
    7470{
    7571  fParticleChange = 0;
    … …  
    8278
    8379G4BetheHeitlerModel::~G4BetheHeitlerModel()
    84 {
    85   if(theCrossSectionTable) {
    86     theCrossSectionTable->clearAndDestroy();
    87     delete theCrossSectionTable;
    88   }
    89 }
    90 
    91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    92 
    93 void G4BetheHeitlerModel::Initialise(const G4ParticleDefinition*,
    94                                      const G4DataVector&)
    95 {
    96   if(!fParticleChange) fParticleChange = GetParticleChangeForGamma();
    97 
    98   if(theCrossSectionTable) {
    99     theCrossSectionTable->clearAndDestroy();
    100     delete theCrossSectionTable;
    101   }
    102 
    103   const G4ElementTable* theElementTable = G4Element::GetElementTable();
    104   size_t nvect = G4Element::GetNumberOfElements();
    105   theCrossSectionTable = new G4PhysicsTable(nvect);
    106   G4PhysicsLogVector* ptrVector;
    107   G4double emin = LowEnergyLimit();
    108   G4double emax = HighEnergyLimit();
    109   G4int n = nbins*G4int(log10(emax/emin));
    110   G4bool spline = G4LossTableManager::Instance()->SplineFlag();
    111   G4double e, value;
    112 
    113   for(size_t j=0; j<nvect ; j++) {
    114 
    115     ptrVector  = new G4PhysicsLogVector(emin, emax, n);
    116     ptrVector->SetSpline(spline);
    117     G4double Z = (*theElementTable)[j]->GetZ();
    118     G4int   iz = G4int(Z);
    119     indexZ[iz] = j;
    120  
    121     for(G4int i=0; i<nbins; i++) {
    122       e = ptrVector->GetLowEdgeEnergy( i ) ;
    123       value = ComputeCrossSectionPerAtom(theGamma, e, Z); 
    124       ptrVector->PutValue( i, value );
    125     }
    126 
    127     theCrossSectionTable->insert(ptrVector);
    128   }
    129 }
    130 
    131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    132 
    133 G4double G4BetheHeitlerModel::ComputeCrossSectionPerAtom(
    134                                                    const G4ParticleDefinition*,
    135                                               G4double GammaEnergy, G4double Z,
    136                                               G4double, G4double, G4double)
     80{}
     81
     82//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     83
     84void G4BetheHeitlerModel::Initialise(const G4ParticleDefinition* p,
     85                                     const G4DataVector& cuts)
     86{
     87  if(!fParticleChange) { fParticleChange = GetParticleChangeForGamma(); }
     88  InitialiseElementSelectors(p, cuts);
     89}
     90
     91//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     92
     93G4double
     94G4BetheHeitlerModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
     95                                                G4double GammaEnergy, G4double Z,
     96                                                G4double, G4double, G4double)
    13797// Calculates the microscopic cross section in GEANT4 internal units.
    13898// A parametrized formula from L. Urban is used to estimate
    … …  
    144104  static const G4double GammaEnergyLimit = 1.5*MeV;
    145105  G4double CrossSection = 0.0 ;
    146   if ( Z < 1. ) return CrossSection;
    147   if ( GammaEnergy <= 2.0*electron_mass_c2 ) return CrossSection;
     106  if ( Z < 0.9 || GammaEnergy <= 2.0*electron_mass_c2 ) { return CrossSection; }
    148107
    149108  static const G4double
    … …  
    160119
    161120  G4double GammaEnergySave = GammaEnergy;
    162   if (GammaEnergy < GammaEnergyLimit) GammaEnergy = GammaEnergyLimit ;
     121  if (GammaEnergy < GammaEnergyLimit) { GammaEnergy = GammaEnergyLimit; }
    163122
    164123  G4double X=log(GammaEnergy/electron_mass_c2), X2=X*X, X3=X2*X, X4=X3*X, X5=X4*X;
    … …  
    177136  }
    178137
    179   if (CrossSection < 0.) CrossSection = 0.;
     138  if (CrossSection < 0.) { CrossSection = 0.; }
    180139  return CrossSection;
    181140}
    … …  
    208167  G4double epsil ;
    209168  G4double epsil0 = electron_mass_c2/GammaEnergy ;
    210   if(epsil0 > 1.0) return;
     169  if(epsil0 > 1.0) { return; }
    211170
    212171  // do it fast if GammaEnergy < 2. MeV
    … …  
    225184    // Extract Coulomb factor for this Element
    226185    G4double FZ = 8.*(anElement->GetIonisation()->GetlogZ3());
    227     if (GammaEnergy > 50.*MeV) FZ += 8.*(anElement->GetfCoulomb());
     186    if (GammaEnergy > 50.*MeV) { FZ += 8.*(anElement->GetfCoulomb()); }
    228187
    229188    // limits of the screening variable
    … …  
    334293
    335294//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    336 
    337 

  • trunk/source/processes/electromagnetic/standard/src/G4BohrFluctuations.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BohrFluctuations.cc,v 1.8 2009/09/29 11:33:22 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4BohrFluctuations.cc,v 1.9 2010/10/25 18:23:36 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6767  xmin(0.2),
    6868  minLoss(0.001*eV)
    69 {}
     69{
     70  particleMass   = proton_mass_c2;
     71  chargeSquare   = 1.0;
     72  kineticEnergy  = 0.0;
     73  beta2          = 0.0;
     74}
    7075
    7176//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    9297                                                      G4double& meanLoss)
    9398{
    94   if(meanLoss <= minLoss) return meanLoss;
     99  if(meanLoss <= minLoss) { return meanLoss; }
    95100  G4double siga = Dispersion(material,dp,tmax,length);
    96101  G4double loss = meanLoss;
    … …  
    142147                                        G4double& length)
    143148{
    144   if(!particle) InitialiseMe(dp->GetDefinition());
     149  if(!particle) { InitialiseMe(dp->GetDefinition()); }
    145150
    146151  G4double electronDensity = material->GetElectronDensity();

  • trunk/source/processes/electromagnetic/standard/src/G4BraggIonModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BraggIonModel.cc,v 1.27 2009/11/22 18:00:23 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4BraggIonModel.cc,v 1.30 2010/11/04 17:30:31 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-25 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8181    isInitialised(false)
    8282{
    83   if(p) SetParticle(p);
    8483  SetHighEnergyLimit(2.0*MeV);
    8584
    … …  
    9089  theZieglerFactor = eV*cm2*1.0e-15;
    9190  theElectron      = G4Electron::Electron();
     91  corrFactor       = 1.0;
     92  if(p) { SetParticle(p); }
     93  else  { SetParticle(theElectron); }
    9294}
    9395
    … …  
    110112                                 const G4DataVector&)
    111113{
    112   if(p != particle) SetParticle(p);
     114  if(p != particle) { SetParticle(p); }
    113115
    114116  corrFactor = chargeSquare;
    … …  
    122124    G4String pname = particle->GetParticleName();
    123125    if(particle->GetParticleType() == "nucleus" &&
    124        pname != "deuteron" && pname != "triton") isIon = true;
     126       pname != "deuteron" && pname != "triton") { isIon = true; }
    125127
    126128    corr = G4LossTableManager::Instance()->EmCorrections();

  • trunk/source/processes/electromagnetic/standard/src/G4BraggModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BraggModel.cc,v 1.26 2010/06/04 09:08:43 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4BraggModel.cc,v 1.28 2010/11/04 17:30:31 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-25 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8484    isInitialised(false)
    8585{
    86   if(p) SetParticle(p);
    8786  SetHighEnergyLimit(2.0*MeV);
    8887
    … …  
    9089  theZieglerFactor = eV*cm2*1.0e-15;
    9190  theElectron = G4Electron::Electron();
     91  expStopPower125 = 0.0;
    9292
    9393  corr = G4LossTableManager::Instance()->EmCorrections();
     94  if(p) { SetParticle(p); }
     95  else  { SetParticle(theElectron); }
    9496}
    9597

  • trunk/source/processes/electromagnetic/standard/src/G4IonFluctuations.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4IonFluctuations.cc,v 1.26 2009/03/31 13:24:40 toshito Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4IonFluctuations.cc,v 1.27 2010/10/25 19:13:23 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8080    xmin(0.2),
    8181    minLoss(0.001*eV)
    82 {}
     82{
     83  kineticEnergy = 0.0;
     84  beta2 = 0.0;
     85}
    8386
    8487//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

  • trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PAIModel.cc,v 1.52 2010/06/03 07:28:39 grichine Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4PAIModel.cc,v 1.54 2010/11/04 17:30:32 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-25 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8181  fBg2lim(0.0169),
    8282  fTaulim(8.4146e-3)
    83 {
    84   if(p) SetParticle(p);
    85  
     83{ 
    8684  fElectron = G4Electron::Electron();
    8785  fPositron = G4Positron::Positron();
    … …  
    9391  fLambdaVector      = 0;
    9492  fdNdxCutVector     = 0;
     93  fParticleEnergyVector = 0;
     94  fSandiaIntervalNumber = 0;
     95  fMatIndex = 0;
     96  fDeltaCutInKinEnergy = 0.0;
     97
     98  if(p) { SetParticle(p); }
     99  else  { SetParticle(fElectron); }
    95100
    96101  isInitialised      = false;
    … …  
    132137void G4PAIModel::SetParticle(const G4ParticleDefinition* p)
    133138{
    134   if(fParticle == p) return;
     139  if(fParticle == p) { return; }
    135140  fParticle = p;
    136141  fMass = fParticle->GetPDGMass();
    … …  
    141146  fRatio = electron_mass_c2/fMass;
    142147  fQc = fMass/fRatio;
     148  fLowestKineticEnergy  = fMass*(fLowestGamma  - 1.0);
     149  fHighestKineticEnergy = fMass*(fHighestGamma - 1.0);
    143150}
    144151
    … …  
    148155                            const G4DataVector&)
    149156{
    150   G4cout<<"G4PAIModel::Initialise for "<<p->GetParticleName()<<G4endl;
    151   if(isInitialised) return;
     157  //G4cout<<"G4PAIModel::Initialise for "<<p->GetParticleName()<<G4endl;
     158  if(isInitialised) { return; }
    152159  isInitialised = true;
    153160
    154161  SetParticle(p);
    155   fLowestKineticEnergy  = fMass*(fLowestGamma  - 1.0);
    156   fHighestKineticEnergy = fMass*(fHighestGamma - 1.0);
    157162
    158163  fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
    … …  
    213218void G4PAIModel::ComputeSandiaPhotoAbsCof()
    214219{
    215   G4int i, j, numberOfElements ;
    216   static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
     220  G4int i, j;
     221  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
    217222
    218223  G4SandiaTable thisMaterialSandiaTable(fMatIndex) ;
    219   numberOfElements = (*theMaterialTable)[fMatIndex]->
    220                                               GetNumberOfElements();
     224  G4int numberOfElements =
     225    (*theMaterialTable)[fMatIndex]->GetNumberOfElements();
     226
    221227  G4int* thisMaterialZ = new G4int[numberOfElements] ;
    222228
    … …  
    236242  fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber] ;
    237243
    238   for(i=0;i<fSandiaIntervalNumber;i++)  fSandiaPhotoAbsCof[i] = new G4double[5] ;
     244  for(i=0; i<fSandiaIntervalNumber; i++) 
     245  {
     246    fSandiaPhotoAbsCof[i] = new G4double[5];
     247  }
    239248   
    240249  for( i = 0 ; i < fSandiaIntervalNumber ; i++ )
    241250  {
    242     fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0) ;
     251    fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0);
    243252
    244253    for( j = 1; j < 5 ; j++ )
    245254    {
    246       fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.
    247                                       GetPhotoAbsorpCof(i+1,j)*
     255      fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,j)*
    248256                 (*theMaterialTable)[fMatIndex]->GetDensity() ;
    249257    }
    250258  }
    251   // delete[] thisMaterialZ ;
     259  delete[] thisMaterialZ;
    252260}
    253261
    … …  
    307315    ionloss = fPAIySection.GetMeanEnergyLoss() ;   //  total <dE/dx>
    308316
    309     if ( ionloss < DBL_MIN)  ionloss = DBL_MIN;
     317    if ( ionloss < DBL_MIN) { ionloss = 0.0; }
    310318    fdEdxVector->PutValue(i,ionloss) ;
    311319

  • trunk/source/processes/electromagnetic/standard/src/G4PAIPhotonModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PAIPhotonModel.cc,v 1.24 2010/06/03 07:28:39 grichine Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4PAIPhotonModel.cc,v 1.25 2010/10/26 09:16:50 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8383  fTaulim(8.4146e-3)
    8484{
    85   if(p) SetParticle(p);
    86 
    8785  fVerbose  = 0;
    8886  fElectron = G4Electron::Electron();
    … …  
    9088
    9189  fProtonEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
    92                                                            fHighestKineticEnergy,
    93                                                            fTotBin);
     90                                               fHighestKineticEnergy,
     91                                               fTotBin);
    9492  fPAItransferTable     = 0;
    9593  fPAIphotonTable       = 0;
    … …  
    104102  fdNdxCutPhotonVector  = 0;
    105103  fdNdxCutPlasmonVector = 0;
     104
     105  fSandiaIntervalNumber = 0;
     106  fMatIndex = 0;
     107
     108  if(p) { SetParticle(p); }
     109  else  { SetParticle(fElectron); }
    106110
    107111  isInitialised      = false;
    … …  
    161165                                   const G4DataVector&)
    162166{
    163   G4cout<<"G4PAIPhotonModel::Initialise for "<<p->GetParticleName()<<G4endl;
    164   if(isInitialised) return;
     167  //  G4cout<<"G4PAIPhotonModel::Initialise for "<<p->GetParticleName()<<G4endl;
     168  if(isInitialised) { return; }
    165169  isInitialised = true;
    166170
    … …  
    229233{
    230234  G4int i, j, numberOfElements ;
    231   static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
     235  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
    232236
    233237  G4SandiaTable thisMaterialSandiaTable(fMatIndex) ;
    … …  
    264268    }
    265269  }
    266   // delete[] thisMaterialZ ;
     270  delete[] thisMaterialZ ;
    267271}
    268272

  • trunk/source/processes/electromagnetic/standard/src/G4PAIySection.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4PAIySection.cc,v 1.4 2009/07/26 15:51:01 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4PAIySection.cc,v 1.6 2010/11/04 17:30:32 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-25 $
    2928//
    3029//
    … …  
    7069
    7170G4PAIySection::G4PAIySection()
    72 {}
     71{
     72  fSandia = 0;
     73  fDensity = fElectronDensity = fNormalizationCof = 0.0;
     74  fIntervalNumber = fSplineNumber = 0;
     75}
    7376
    7477////////////////////////////////////////////////////////////////////////////

  • trunk/source/processes/electromagnetic/standard/src/G4PairProductionRelModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PairProductionRelModel.cc,v 1.3 2009/05/15 17:12:33 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4PairProductionRelModel.cc,v 1.4 2010/10/26 09:06:04 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8585                                         const G4String& nam)
    8686  : G4VEmModel(nam),
    87     theCrossSectionTable(0),
    88     nbins(10),
    8987    fLPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)*0.5),
    9088    fLPMflag(true),
    … …  
    9997  nist = G4NistManager::Instance(); 
    10098
     99  currentZ = z13 = z23 = lnZ = Fel = Finel = fCoulomb = phiLPM = gLPM = xiLPM = 0;
     100
    101101}
    102102
    … …  
    104104
    105105G4PairProductionRelModel::~G4PairProductionRelModel()
    106 {
    107   if(theCrossSectionTable) {
    108     theCrossSectionTable->clearAndDestroy();
    109     delete theCrossSectionTable;
    110   }
    111 }
    112 
    113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    114 
    115 void G4PairProductionRelModel::Initialise(const G4ParticleDefinition*,
    116                                           const G4DataVector&)
    117 {
    118   fParticleChange = GetParticleChangeForGamma();
    119 
    120   if(theCrossSectionTable) {
    121     theCrossSectionTable->clearAndDestroy();
    122     delete theCrossSectionTable;
    123   }
    124 
    125   const G4ElementTable* theElementTable = G4Element::GetElementTable();
    126   size_t nvect = G4Element::GetNumberOfElements();
    127   theCrossSectionTable = new G4PhysicsTable(nvect);
    128   G4PhysicsLogVector* ptrVector;
    129   G4double emin = LowEnergyLimit();
    130   G4double emax = HighEnergyLimit();
    131   G4int n = nbins*G4int(log10(emax/emin));
    132   G4bool spline = G4LossTableManager::Instance()->SplineFlag();
    133   G4double e, value;
    134 
    135   for(size_t j=0; j<nvect ; j++) {
    136 
    137     ptrVector  = new G4PhysicsLogVector(emin, emax, n);
    138     ptrVector->SetSpline(spline);
    139     G4double Z = (*theElementTable)[j]->GetZ();
    140     G4VEmModel::SetCurrentElement((*theElementTable)[j]);
    141     G4int   iz = G4int(Z);
    142     indexZ[iz] = j;
    143  
    144     for(G4int i=0; i<nbins; i++) {
    145       e = ptrVector->GetLowEdgeEnergy( i ) ;
    146       value = ComputeCrossSectionPerAtom(theGamma, e, Z); 
    147       ptrVector->PutValue( i, value );
    148     }
    149 
    150     theCrossSectionTable->insert(ptrVector);
    151   }
    152 }
    153 
    154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    155 /*
    156 G4double G4PairProductionRelModel::ComputeRelXSectionPerAtom(G4double k, G4double Z)
    157 {
    158 
    159   G4double cross = 0.0;
    160  
    161 
    162 
    163 
    164 }
    165 */
     106{}
     107
     108//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     109
     110void G4PairProductionRelModel::Initialise(const G4ParticleDefinition* p,
     111                                          const G4DataVector& cuts)
     112{
     113  if(!fParticleChange) { fParticleChange = GetParticleChangeForGamma(); }
     114  InitialiseElementSelectors(p, cuts);
     115}
    166116
    167117//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    340290  //  static const G4double gammaEnergyLimit = 1.5*MeV;
    341291  G4double crossSection = 0.0 ;
    342   if ( Z < 1. ) return crossSection;
     292  if ( Z < 0.9 ) return crossSection;
    343293  if ( gammaEnergy <= 2.0*electron_mass_c2 ) return crossSection;
    344294

  • trunk/source/processes/electromagnetic/standard/src/G4UniversalFluctuation.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4UniversalFluctuation.cc,v 1.22 2009/03/20 18:11:23 urban Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4UniversalFluctuation.cc,v 1.28 2010/10/26 10:06:12 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    3434// File name:     G4UniversalFluctuation
    3535//
    36 // Author:        Vladimir Ivanchenko
     36// Author:        Laszlo Urban
    3737//
    3838// Creation date: 03.01.2002
    … …  
    5757// 03-04-07 correction to get better width of eloss distr.(L.Urban)
    5858// 13-07-07 add protection for very small step or low-density material (VI)
    59 // 19-03-09 new width correction (does not depend on previous steps) (L.Urban)         
     59// 19-03-09 new width correction (does not depend on previous steps) (L.Urban)
    6060// 20-03-09 modification in the width correction (L.Urban)
     61// 14-06-10 fixed tail distribution - do not use uniform function (L.Urban)
     62// 08-08-10 width correction algorithm has bee modified -->
     63//          better results for thin targets (L.Urban)
    6164//
    6265
    … …  
    8285  theBohrBeta2(50.0*keV/proton_mass_c2),
    8386  minLoss(10.*eV),
    84   nmaxCont1(4.),
    85   nmaxCont2(16.)
     87  nmaxCont(16.),
     88  rate(0.55),
     89  fw(4.)
    8690{
    8791  lastMaterial = 0;
     92
     93  particleMass = chargeSquare = ipotFluct = electronDensity = f1Fluct = f2Fluct
     94    = e1Fluct = e2Fluct = e1LogFluct = e2LogFluct = ipotLogFluct = e0 = esmall
     95    = e1 = e2 = 0;
    8896}
    8997
    … …  
    118126  // shortcut for very very small loss (out of validity of the model)
    119127  //
    120   if (meanLoss < minLoss)
    121     return meanLoss;
    122 
    123   if(!particle) InitialiseMe(dp->GetDefinition());
     128  if (meanLoss < minLoss) { return meanLoss; }
     129
     130  if(!particle) { InitialiseMe(dp->GetDefinition()); }
    124131
    125132  G4double tau   = dp->GetKineticEnergy()/particleMass;
    … …  
    174181    ipotLogFluct = material->GetIonisation()->GetLogMeanExcEnergy();
    175182    e0 = material->GetIonisation()->GetEnergy0fluct();
     183    esmall = 0.5*sqrt(e0*ipotFluct); 
    176184    lastMaterial = material;
    177  
    178     // modification of some model parameters
    179     // (this part should go to materials later)
    180     G4double p = 1.40;
    181     f2Fluct *= p;
    182     f1Fluct = 1.-f2Fluct;
    183     G4double q = 1.00;
    184     e2Fluct *= q;
    185     e2LogFluct = log(e2Fluct);
    186     e1LogFluct = (ipotLogFluct-f2Fluct*e2LogFluct)/f1Fluct;
    187     e1Fluct = exp(e1LogFluct);
     185   
    188186  }
    189187
    … …  
    193191  G4double a1 = 0. , a2 = 0., a3 = 0. ;
    194192
    195   // cut and material dependent rate
    196   G4double rate = 1.0;
    197193  if(tmax > ipotFluct) {
    198194    G4double w2 = log(2.*electron_mass_c2*beta2*gam2)-beta2;
    199195
    200     if(w2 > ipotLogFluct && w2 > e2LogFluct) {
    201 
    202       rate = 0.03+0.23*log(log(tmax/ipotFluct));
     196    if(w2 > ipotLogFluct && w2 > e2LogFluct && tmax> ipotFluct)  {
    203197      G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
    204198      a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
    205199      a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
    206       // correction in order to get better FWHM values
    207       // ( scale parameters a1 and e1)
    208       G4double width = 1.;
    209       if(meanLoss > 10.*e1Fluct)
     200         
     201
     202  if(a1 < nmaxCont)
     203  {
     204    //small energy loss
     205    G4double sa1 = sqrt(a1);
     206    if(G4UniformRand() < exp(-sa1))
    210207      {
    211         width = 3.1623/sqrt(meanLoss/e1Fluct);
    212         if(width < a2/a1)
    213         width = a2/a1;
     208       e1 = esmall;
     209       a1 = meanLoss*(1.-rate)/e1;
     210       a2 = 0.;
     211       e2 = e2Fluct;
    214212      }
    215       a1 *= width;
    216       e1 = e1Fluct/width;
     213      else
     214      {
     215       a1 = sa1 ;   
     216       e1 = sa1*e1Fluct;
     217       e2 = e2Fluct;
     218      }
     219    }
     220    else
     221    {
     222      //not small energy loss
     223      //correction to get better fwhm value
     224      a1 /= fw;
     225      e1 = fw*e1Fluct;
    217226      e2 = e2Fluct;
    218227    }
    219   }
     228  }   
     229 }
    220230
    221231  G4double w1 = tmax/e0;
    … …  
    223233    a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*log(w1));
    224234
    225   //'nearly' Gaussian fluctuation if a1>nmaxCont2&&a2>nmaxCont2&&a3>nmaxCont2 
     235  //'nearly' Gaussian fluctuation if a1>nmaxCont&&a2>nmaxCont&&a3>nmaxCont 
    226236  G4double emean = 0.;
    227237  G4double sig2e = 0., sige = 0.;
    … …  
    229239 
    230240  // excitation of type 1
    231   if(a1 > nmaxCont2)
     241  if(a1 > nmaxCont)
    232242  {
    233243    emean += a1*e1;
    … …  
    243253
    244254  // excitation of type 2
    245   if(a2 > nmaxCont2)
     255  if(a2 > nmaxCont)
    246256  {
    247257    emean += a2*e2;
    … …  
    262272    p3 = a3;
    263273    G4double alfa = 1.;
    264     if(a3 > nmaxCont2)
     274    if(a3 > nmaxCont)
    265275    {
    266        alfa            = w1*(nmaxCont2+a3)/(w1*nmaxCont2+a3);
     276       alfa            = w1*(nmaxCont+a3)/(w1*nmaxCont+a3);
    267277       G4double alfa1  = alfa*log(alfa)/(alfa-1.);
    268278       G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);

  • trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel90.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4UrbanMscModel90.cc,v 1.13 2009/04/10 18:10:58 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4UrbanMscModel90.cc,v 1.14 2010/10/26 10:06:12 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9393  inside        = false; 
    9494  insideskin    = false;
     95
     96  skindepth = skin*stepmin;
     97
     98  mass = proton_mass_c2;
     99  charge = 1.0;
     100  currentKinEnergy = currentRange = currentRadLength = masslimite = masslimitmu
     101    = lambda0 = lambdaeff = tPathLength = zPathLength = par1 = par2 = par3 = 0;
     102
     103  currentMaterialIndex = 0;
    95104}
    96105
    … …  
    105114                                   const G4DataVector&)
    106115{
    107   skindepth     = skin*stepmin;
     116  skindepth = skin*stepmin;
    108117  if(isInitialized) return;
    109118

  • trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel92.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4UrbanMscModel92.cc,v 1.1 2009/11/01 13:05:01 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4UrbanMscModel92.cc,v 1.2 2010/10/26 10:06:12 vnivanch Exp $
     28// GEANT4 tag $Name: emstand-V09-03-24 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    166166  insideskin    = false;
    167167
     168  skindepth = skin*stepmin;
     169
     170  mass = proton_mass_c2;
     171  charge = ChargeSquare = 1.0;
     172  currentKinEnergy = currentRadLength = lambda0 = lambdaeff = tPathLength
     173    = zPathLength = par1 = par2 = par3 = 0;
     174
     175  currentMaterialIndex = 0;
    168176}
    169177

  • trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel93.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4UrbanMscModel93.cc,v 1.5 2010/06/25 09:41:44 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4UrbanMscModel93.cc,v 1.7 2010/11/04 17:30:32 vnivanch Exp $
     28// GEANT4 tag $Name: emstand-V09-03-25 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    166166  insideskin    = false;
    167167
     168  skindepth = skin*stepmin;
     169
     170  mass = proton_mass_c2;
     171  charge = ChargeSquare = 1.0;
     172  currentKinEnergy = currentRadLength = lambda0 = lambdaeff = tPathLength
     173    = zPathLength = par1 = par2 = par3 = 0;
     174
     175  currentMaterialIndex = -1;
    168176}
    169177
    … …  
    179187{
    180188  skindepth = skin*stepmin;
    181   if(isInitialized) return;
     189
     190  if(isInitialized) { return; }
    182191  // set values of some data members
    183192  SetParticle(p);
    … …  
    460469
    461470  // stop here if small range particle
    462   if(inside) return tPathLength;           
     471  if(inside) { return tPathLength; }           
    463472 
    464   if(tPathLength > currentRange) tPathLength = currentRange;
     473  if(tPathLength > currentRange) { tPathLength = currentRange; }
    465474
    466475  presafety = sp->GetSafety();
    467476
    468  // G4cout << "G4Urban2::StepLimit tPathLength= "
    469  //      <<tPathLength<<" safety= " << presafety
    470  //        << " range= " <<currentRange<< " lambda= "<<lambda0
    471  //      << " Alg: " << steppingAlgorithm <<G4endl;
     477  // G4cout << "G4Urban2::StepLimit tPathLength= "
     478  //     <<tPathLength<<" safety= " << presafety
     479  //        << " range= " <<currentRange<< " lambda= "<<lambda0
     480  //     << " Alg: " << steppingAlgorithm <<G4endl;
    472481
    473482  // far from geometry boundary

  • trunk/source/processes/electromagnetic/standard/src/G4WaterStopping.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4WaterStopping.cc,v 1.21 2010/04/26 17:44:34 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4WaterStopping.cc,v 1.22 2010/10/26 10:06:12 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828
    2929//---------------------------------------------------------------------------
    … …  
    6868{
    6969  G4double res = 0.0;
    70   if((iz > 26) || (iz < 3) || (iz > 18 && iz < 26)) { return res; }
    7170  G4int idx = iz - 3;
     71
     72  if(iz == 26) { idx = 16; }
     73  else if (iz < 3 || iz > 18) { return res; }
     74
    7275  G4double scaledEnergy = energy/A[idx];
    7376  if(scaledEnergy < emin) {

  • trunk/source/processes/electromagnetic/standard/src/G4WentzelOKandVIxSection.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4WentzelOKandVIxSection.cc,v 1.10 2010/06/01 13:34:21 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4WentzelOKandVIxSection.cc,v 1.12 2010/11/04 17:30:32 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-25 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9393    }
    9494  }
     95  currentMaterial = 0;
     96  elecXSRatio = factB = formfactA = screenZ = 0.0;
     97  cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = 1.0;
     98
     99  Initialise(theElectron, 1.0);
     100  SetTargetMass(proton_mass_c2);
    95101}
    96102

  • trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4WentzelVIModel.cc,v 1.60 2010/06/01 11:13:31 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4WentzelVIModel.cc,v 1.61 2010/10/26 10:06:12 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    8989  fG4pow = G4Pow::GetInstance();
    9090  wokvi = new G4WentzelOKandVIxSection();
     91
     92  preKinEnergy = tPathLength = zPathLength = lambdaeff = currentRange = xtsec = 0;
     93  currentMaterialIndex = 0;
     94  cosThetaMax = cosTetMaxNuc = 1.0;
    9195}
    9296

  • trunk/source/processes/electromagnetic/standard/src/G4alphaIonisation.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4alphaIonisation.cc,v 1.1 2009/11/10 11:50:30 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4alphaIonisation.cc,v 1.3 2010/10/26 10:06:12 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6464  : G4VEnergyLossProcess(name),
    6565    theParticle(0),
    66     isInitialised(false),
    67     nuclearStopping(true)
     66    isInitialised(false)
    6867{
    6968  //  SetLinearLossLimit(0.15);
    … …  
    7473  mass = 0.0;
    7574  ratio = 0.0;
     75  eth = 8*MeV;
    7676}
    7777
    … …  
    8585G4bool G4alphaIonisation::IsApplicable(const G4ParticleDefinition& p)
    8686{
    87   return (p.GetPDGCharge() == 2*eplus);
     87  return (!p.IsShortLived() &&
     88          std::fabs(p.GetPDGCharge() - 2*CLHEP::eplus) < 0.01);
    8889}
    8990
    … …  
    122123    SetSecondaryParticle(G4Electron::Electron());
    123124
    124     if (!EmModel(1)) SetEmModel(new G4BraggIonModel(), 1);
     125    if (!EmModel(1)) { SetEmModel(new G4BraggIonModel(), 1); }
    125126    EmModel(1)->SetLowEnergyLimit(MinKinEnergy());
    126127
    … …  
    129130    EmModel(1)->SetHighEnergyLimit(eth);
    130131
    131     if (!FluctModel()) SetFluctModel(new G4UniversalFluctuation());
     132    if (!FluctModel()) { SetFluctModel(new G4UniversalFluctuation()); }
    132133    AddEmModel(1, EmModel(1), new G4IonFluctuations());
    133134
    134     if (!EmModel(2)) SetEmModel(new G4BetheBlochModel(),2); 
     135    if (!EmModel(2)) { SetEmModel(new G4BetheBlochModel(),2); } 
    135136    EmModel(2)->SetLowEnergyLimit(eth);
    136137    EmModel(2)->SetHighEnergyLimit(MaxKinEnergy());
    … …  
    139140    isInitialised = true;
    140141  }
    141   // reinitialisation of corrections for the new run
    142   EmModel(1)->ActivateNuclearStopping(nuclearStopping);
    143   EmModel(2)->ActivateNuclearStopping(nuclearStopping);
    144142}
    145143
    … …  
    147145
    148146void G4alphaIonisation::PrintInfo()
    149 {
    150   if (G4Alpha::Alpha() == theParticle) {
    151     if(EmModel(1) && EmModel(2)) {
    152       G4cout << "      NuclearStopping= " << nuclearStopping
    153              << G4endl;
    154     }
    155   }
    156 }
     147{}
    157148
    158149//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

  • trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4eBremsstrahlungModel.cc,v 1.46 2010/04/28 18:39:40 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4eBremsstrahlungModel.cc,v 1.48 2010/10/26 10:35:22 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5555// 15-02-07  correct LPMconstant by a factor 2, thanks to G. Depaola (mma)
    5656// 09-09-08  MigdalConstant increased in (2pi)^2 times (A.Schaelicke)
     57// 13-10-10  Add angular distributon interface (VI)
    5758//
    5859// Class Description:
    … …  
    7576#include "G4DataVector.hh"
    7677#include "G4ParticleChangeForLoss.hh"
     78#include "G4ModifiedTsai.hh"
    7779
    7880//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    9092    isInitialised(false)
    9193{
    92   if(p) SetParticle(p);
     94  if(p) { SetParticle(p); }
    9395  theGamma = G4Gamma::Gamma();
    9496  minThreshold = 0.1*keV;
     97  SetAngularDistribution(new G4ModifiedTsai());
     98  highKinEnergy = HighEnergyLimit();
     99  lowKinEnergy  = LowEnergyLimit();
    95100}
    96101
    … …  
    129134                                        const G4DataVector& cuts)
    130135{
    131   if(p) SetParticle(p);
     136  if(p) { SetParticle(p); }
    132137  highKinEnergy = HighEnergyLimit();
    133138  lowKinEnergy  = LowEnergyLimit();
    … …  
    172177                                                   G4double cutEnergy)
    173178{
    174   if(!particle) SetParticle(p);
    175   if(kineticEnergy < lowKinEnergy) return 0.0;
     179  if(!particle) { SetParticle(p); }
     180  if(kineticEnergy < lowKinEnergy) { return 0.0; }
    176181
    177182  const G4double thigh = 100.*GeV;
    … …  
    272277    dedx += loss;
    273278  }
    274   if(dedx < 0.) dedx = 0.;
     279  if(dedx < 0.) { dedx = 0.; }
    275280  return dedx;
    276281}
    … …  
    413418                                                    G4double maxEnergy)
    414419{
    415   if(!particle) SetParticle(p);
     420  if(!particle) { SetParticle(p); }
    416421  G4double cross = 0.0;
    417422  G4double tmax = min(maxEnergy, kineticEnergy);
    418423  G4double cut  = max(cutEnergy, minThreshold);
    419   if(cut >= tmax) return cross;
     424  if(cut >= tmax) { return cross; }
    420425
    421426  const G4ElementVector* theElementVector = material->GetElementVector();
    … …  
    495500{
    496501  G4double cross = 0.0 ;
    497   if ( kineticEnergy < 1*keV || kineticEnergy < cut) return cross;
     502  if ( kineticEnergy < 1*keV || kineticEnergy < cut) { return cross; }
    498503
    499504  static const G4double ksi=2.0, alfa=1.00;
    … …  
    664669  G4double kineticEnergy = dp->GetKineticEnergy();
    665670  G4double tmax = min(maxEnergy, kineticEnergy);
    666   if(tmin >= tmax) return;
     671  if(tmin >= tmax) { return; }
    667672
    668673//
    … …  
    709714  G4double xmax     = tmax/kineticEnergy;
    710715  G4double kappa    = 0.0;
    711   if(xmax >= 1.) xmax = 1.;
    712   else     kappa    = log(xmax)/log(xmin);
     716  if(xmax >= 1.) { xmax = 1.; }
     717  else   { kappa    = log(xmax)/log(xmin); }
    713718  G4double epsilmin = tmin/totalEnergy;
    714719  G4double epsilmax = tmax/totalEnergy;
    … …  
    812817      } while( greject < G4UniformRand()*grejmax );
    813818    }
    814     /*
    815     if(x > 0.999) {
    816       G4cout << "### G4eBremsstrahlungModel Warning: e= " << kineticEnergy
    817              << " tlow= " << tlow
    818              << " x= " << x
    819              << " greject= " << greject
    820              << " grejmax= " << grejmax
    821              << " migdal= " << migdal
    822              << G4endl;
    823       //      if(x >= 1.0) G4Exception("X=1");
    824     }
    825     */
    826819    gammaEnergy = x*kineticEnergy;
    827820
    … …  
    837830
    838831  //
    839   //  angles of the emitted gamma. ( Z - axis along the parent particle)
     832  // angles of the emitted gamma. ( Z - axis along the parent particle)
     833  // use general interface
    840834  //
    841   //  universal distribution suggested by L. Urban
    842   // (Geant3 manual (1993) Phys211),
    843   //  derived from Tsai distribution (Rev Mod Phys 49,421(1977))
    844 
    845   G4double u;
    846   const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
    847 
    848   if (9./(9.+d) > G4UniformRand()) u = - log(G4UniformRand()*G4UniformRand())/a1;
    849      else                          u = - log(G4UniformRand()*G4UniformRand())/a2;
    850 
    851   G4double theta = u*electron_mass_c2/totalEnergy;
     835  G4double theta = GetAngularDistribution()->PolarAngle(totalEnergy,
     836                                                        totalEnergy-gammaEnergy,
     837                                                        (G4int)anElement->GetZ());
    852838
    853839  G4double sint = sin(theta);

  • trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4eBremsstrahlungRelModel.cc,v 1.15 2010/04/06 17:02:23 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4eBremsstrahlungRelModel.cc,v 1.18 2010/11/04 17:30:32 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-25 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4242// 13.11.08    add SetLPMflag and SetLPMconstant methods
    4343// 13.11.08    change default LPMconstant value
     44// 13.10.10    add angular distributon interface (VI)
    4445//
    4546// Main References:
    … …  
    6566#include "G4ParticleChangeForLoss.hh"
    6667#include "G4LossTableManager.hh"
    67 
     68#include "G4ModifiedTsai.hh"
    6869
    6970//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    9091    use_completescreening(true),isInitialised(false)
    9192{
    92   if(p) SetParticle(p);
    9393  theGamma = G4Gamma::Gamma();
    9494
    9595  minThreshold = 0.1*keV;
    96   SetLowEnergyLimit(GeV); 
     96  lowKinEnergy = GeV;
     97  SetLowEnergyLimit(lowKinEnergy); 
    9798
    9899  nist = G4NistManager::Instance(); 
     100
     101  SetLPMFlag(true);
     102  SetAngularDistribution(new G4ModifiedTsai());
     103
     104  particleMass = kinEnergy = totalEnergy = currentZ = z13 = z23 = lnZ = Fel
     105    = Finel = fCoulomb = fMax = densityFactor = densityCorr = lpmEnergy
     106    = xiLPM = phiLPM = gLPM = klpm = kp = 0.0;
     107
     108  energyThresholdLPM = 1.e39;
     109
    99110  InitialiseConstants();
    100 
    101   SetLPMFlag(true);
     111  if(p) { SetParticle(p); }
    102112}
    103113
    … …  
    125135  particle = p;
    126136  particleMass = p->GetPDGMass();
    127   if(p == G4Electron::Electron()) isElectron = true;
    128   else                            isElectron = false;
     137  if(p == G4Electron::Electron()) { isElectron = true; }
     138  else                            { isElectron = false;}
    129139}
    130140
    … …  
    140150
    141151void G4eBremsstrahlungRelModel::SetupForMaterial(const G4ParticleDefinition*,
    142                                                  const G4Material* mat, G4double kineticEnergy)
     152                                                 const G4Material* mat,
     153                                                 G4double kineticEnergy)
    143154{
    144155  densityFactor = mat->GetElectronDensity()*fMigdalConstant;
    … …  
    146157
    147158  // Threshold for LPM effect (i.e. below which LPM hidden by density effect)
    148   if (LPMFlag())
     159  if (LPMFlag()) {
    149160    energyThresholdLPM=sqrt(densityFactor)*lpmEnergy;
    150   else
     161  } else {
    151162     energyThresholdLPM=1.e39;   // i.e. do not use LPM effect
    152 
     163  }
    153164  // calculate threshold for density effect
    154165  kinEnergy   = kineticEnergy;
    … …  
    168179                                           const G4DataVector& cuts)
    169180{
    170   if(p) SetParticle(p);
    171 
    172   highKinEnergy = HighEnergyLimit();
     181  if(p) { SetParticle(p); }
     182
    173183  lowKinEnergy  = LowEnergyLimit();
    174184
    … …  
    177187  InitialiseElementSelectors(p, cuts);
    178188
    179   if(isInitialised) return;
     189  if(isInitialised) { return; }
    180190  fParticleChange = GetParticleChangeForLoss();
    181191  isInitialised = true;
    … …  
    190200                                                   G4double cutEnergy)
    191201{
    192   if(!particle) SetParticle(p);
    193   if(kineticEnergy < lowKinEnergy) return 0.0;
     202  if(!particle) { SetParticle(p); }
     203  if(kineticEnergy < lowKinEnergy) { return 0.0; }
    194204  G4double cut = std::min(cutEnergy, kineticEnergy);
    195   if(cut == 0.0) return 0.0;
     205  if(cut == 0.0) { return 0.0; }
    196206
    197207  SetupForMaterial(particle, material,kineticEnergy);
    … …  
    261271                                              G4double maxEnergy)
    262272{
    263   if(!particle) SetParticle(p);
    264   if(kineticEnergy < lowKinEnergy) return 0.0;
     273  if(!particle) { SetParticle(p); }
     274  if(kineticEnergy < lowKinEnergy) { return 0.0; }
    265275  G4double cut  = std::min(cutEnergy, kineticEnergy);
    266276  G4double tmax = std::min(maxEnergy, kineticEnergy);
    267277
    268   if(cut >= tmax) return 0.0;
     278  if(cut >= tmax) { return 0.0; }
    269279
    270280  SetCurrentElement(Z);
    … …  
    273283
    274284  // allow partial integration
    275   if(tmax < kinEnergy) cross -= ComputeXSectionPerAtom(tmax);
     285  if(tmax < kinEnergy) { cross -= ComputeXSectionPerAtom(tmax); }
    276286 
    277287  cross *= Z*Z*bremFactor;
    … …  
    390400  // *** make sure suppression is smaller than 1 ***
    391401  // *** caused by Migdal approximation in xi    ***
    392   if (xiLPM*phiLPM>1. || s>0.57)  xiLPM=1./phiLPM;
     402  if (xiLPM*phiLPM>1. || s>0.57)  { xiLPM=1./phiLPM; }
    393403}
    394404
    … …  
    401411//    * complete screening
    402412{
    403   if(gammaEnergy < 0.0) return 0.0;
     413  if(gammaEnergy < 0.0) { return 0.0; }
    404414
    405415  G4double y = gammaEnergy/totalEnergy;
    … …  
    429439{
    430440
    431   if(gammaEnergy < 0.0) return 0.0;
     441  if(gammaEnergy < 0.0) { return 0.0; }
    432442
    433443  G4double y = gammaEnergy/totalEnergy;
    … …  
    465475{
    466476  G4double kineticEnergy = dp->GetKineticEnergy();
    467   if(kineticEnergy < lowKinEnergy) return;
     477  if(kineticEnergy < lowKinEnergy) { return; }
    468478  G4double cut  = std::min(cutEnergy, kineticEnergy);
    469479  G4double emax = std::min(maxEnergy, kineticEnergy);
    470   if(cut >= emax) return;
     480  if(cut >= emax) { return; }
    471481
    472482  SetupForMaterial(particle, couple->GetMaterial(),kineticEnergy);
    … …  
    483493  //  G4double fmax= fMax;
    484494  G4bool highe = true;
    485   if(totalEnergy < energyThresholdLPM) highe = false;
     495  if(totalEnergy < energyThresholdLPM) { highe = false; }
    486496 
    487497  G4double xmin = log(cut*cut + densityCorr);
    … …  
    507517
    508518  //
    509   //  angles of the emitted gamma. ( Z - axis along the parent particle)
     519  // angles of the emitted gamma. ( Z - axis along the parent particle)
     520  // use general interface
    510521  //
    511   //  universal distribution suggested by L. Urban
    512   // (Geant3 manual (1993) Phys211),
    513   //  derived from Tsai distribution (Rev Mod Phys 49,421(1977))
    514 
    515   G4double u;
    516   const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
    517 
    518   if (9./(9.+d) > G4UniformRand()) u = - log(G4UniformRand()*G4UniformRand())/a1;
    519      else                          u = - log(G4UniformRand()*G4UniformRand())/a2;
    520 
    521   G4double theta = u*particleMass/totalEnergy;
     522  G4double theta = GetAngularDistribution()->PolarAngle(totalEnergy,
     523                                                        totalEnergy-gammaEnergy,
     524                                                        (G4int)currentZ);
     525
    522526  G4double sint = sin(theta);
    523527  G4double phi = twopi * G4UniformRand();

  • trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4eCoulombScatteringModel.cc,v 1.89 2010/05/27 14:22:05 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4eCoulombScatteringModel.cc,v 1.90 2010/10/26 10:35:22 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    9191  currentCouple = 0;
    9292  wokvi = new G4WentzelOKandVIxSection();
     93
     94  currentMaterialIndex = 0;
     95
     96  cosTetMinNuc = 1.0;
     97  cosTetMaxNuc = -1.0;
     98  elecRatio = 0.0;
     99  mass = proton_mass_c2;
    93100}
    94101

  • trunk/source/processes/electromagnetic/standard/src/G4hIonisation.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4hIonisation.cc,v 1.85 2010/06/04 09:22:14 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4hIonisation.cc,v 1.86 2010/10/26 10:42:04 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    116116  mass = 0.0;
    117117  ratio = 0.0;
     118  eth = 2*MeV;
    118119}
    119120

  • trunk/source/processes/electromagnetic/standard/src/G4ionIonisation.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4ionIonisation.cc,v 1.70 2009/11/27 20:06:32 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4ionIonisation.cc,v 1.72 2010/10/26 10:42:04 vnivanch Exp $
     27// GEANT4 tag $Name: emstand-V09-03-24 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    6666#include "G4Electron.hh"
    6767#include "G4Proton.hh"
    68 //#include "G4Alpha.hh"
    6968#include "G4GenericIon.hh"
    7069#include "G4BraggModel.hh"
    … …  
    8382G4ionIonisation::G4ionIonisation(const G4String& name)
    8483  : G4VEnergyLossProcess(name),
    85     corr(0),
    8684    theParticle(0),
    8785    isInitialised(false),
    88     stopDataActive(true),
    89     nuclearStopping(true)
     86    stopDataActive(true)
    9087{
    9188  SetLinearLossLimit(0.02);
    … …  
    9592  //  SetVerboseLevel(1);
    9693  corr = G4LossTableManager::Instance()->EmCorrections();
     94  eth = 2*MeV;
    9795}
    9896
    … …  
    136134    const G4ParticleDefinition* theBaseParticle = 0;
    137135
    138     if(part == ion)     theBaseParticle = 0;
    139     else if(bpart == 0) theBaseParticle = ion;
    140     else                theBaseParticle = bpart;
     136    if(part == ion)     { theBaseParticle = 0; }
     137    else if(bpart == 0) { theBaseParticle = ion; }
     138    else                { theBaseParticle = bpart; }
    141139
    142140    SetBaseParticle(theBaseParticle);
    143141    SetSecondaryParticle(G4Electron::Electron());
    144142
    145     if (!EmModel(1)) SetEmModel(new G4BraggIonModel(), 1);
     143    if (!EmModel(1)) { SetEmModel(new G4BraggIonModel(), 1); }
    146144    EmModel(1)->SetLowEnergyLimit(MinKinEnergy());
    147145
    … …  
    150148    EmModel(1)->SetHighEnergyLimit(eth);
    151149
    152     if (!FluctModel()) SetFluctModel(new G4IonFluctuations());
     150    if (!FluctModel()) { SetFluctModel(new G4IonFluctuations()); }
    153151    AddEmModel(1, EmModel(1), FluctModel());
    154152
    155     if (!EmModel(2)) SetEmModel(new G4BetheBlochModel(),2); 
     153    if (!EmModel(2)) { SetEmModel(new G4BetheBlochModel(),2); } 
    156154    EmModel(2)->SetLowEnergyLimit(eth);
    157155    EmModel(2)->SetHighEnergyLimit(MaxKinEnergy());
    … …  
    166164  }
    167165  // reinitialisation of corrections for the new run
    168   EmModel(1)->ActivateNuclearStopping(nuclearStopping);
    169   EmModel(2)->ActivateNuclearStopping(nuclearStopping);
    170   if(part == ion) corr->InitialiseForNewRun();
     166  if(part == ion) { corr->InitialiseForNewRun(); }
    171167}
    172168
    … …  
    178174    G4cout << "      Stopping Power data for "
    179175           << corr->GetNumberOfStoppingVectors()
    180            << " ion/material pairs, nuclearStopping: " << nuclearStopping
     176           << " ion/material pairs "
    181177           << G4endl;
    182178  }

  • 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

  • trunk/source/processes/electromagnetic/xrays/History

    r1337 r1340  
    1 $Id: History,v 1.72 2010/06/16 15:34:31 gcosmo Exp $
     1$Id: History,v 1.77 2010/11/03 19:22:54 gum Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     2003 November 10: P. Gumplinger (xrays-V09-03-05)
     21- more corrections to G4Scintillation; now all particles are assigned the
     22  ELECTRONSCINTILLATIONYIELD unless the user specifies otherwise
     23  sort out where Ion and Neutron (recoil ions below tracking cut) should go
     24
     2528 October 10: P. Gumplinger (xrays-V09-03-04)
     26- several small corrections to the previous implementation of G4Scintillation
     27
     2819 October 10: P. Gumplinger (xrays-V09-03-03)
     29- G4Scintillation: allow for the light yield to be a function of
     30  particle type and deposited energy in case of non-linear light
     31  emission in scintillators. Thanks to Zach Hartwig (Department
     32  of Nuclear Science and Engineeering - MIT)
     33
     3414 October 10: V. Ivanchenko (xrays-V09-03-02)
     35- G4SynchrotronRadiationInMat, G4SynchrotronRadiation,
     36  G4TransitionRadiation - added process sub-types
    1937
    203816 June 10: G. Cosmo (xrays-V09-03-01)

  • trunk/source/processes/electromagnetic/xrays/include/G4Scintillation.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4Scintillation.hh,v 1.18 2010/06/25 09:41:46 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4Scintillation.hh,v 1.21 2010/10/28 23:29:21 gum Exp $
     28// GEANT4 tag $Name: xrays-V09-03-05 $
    2929//
    3030//
    … …  
    3838// Created:     1998-11-07
    3939// Author:      Peter Gumplinger
    40 // Updated:     2005-07-28 add G4ProcessType to constructor
     40// Updated:     2010-10-20 Allow the scintillation yield to be a function
     41//                         of energy deposited by particle type
     42//                         Thanks to Zach Hartwig (Department of Nuclear
     43//                         Science and Engineeering - MIT)
     44//              2005-07-28 add G4ProcessType to constructor
    4145//              2002-11-21 change to user G4Poisson for small MeanNumPotons
    4246//              2002-11-07 allow for fast and slow scintillation
    … …  
    182186        // Adds Birks Saturation to the process.
    183187
     188        void RemoveSaturation() { emSaturation = NULL; }
     189        // Removes the Birks Saturation from the process.
     190
    184191        G4EmSaturation* GetSaturation() const { return emSaturation; }
    185192        // Returns the Birks Saturation.
     193
     194        void SetScintillationByParticleType(const G4bool );
     195        // Called by the user to set the scintillation yield as a function
     196        // of energy deposited by particle type
     197
     198        G4bool GetScintillationByParticleType() const
     199        { return scintillationByParticleType; }
     200        // Return the boolean that determines the method of scintillation
     201        // production
    186202
    187203        void DumpPhysicsTable() const;
    … …  
    210226
    211227        G4double ExcitationRatio;
     228
     229        G4bool scintillationByParticleType;
    212230
    213231private:

  • trunk/source/processes/electromagnetic/xrays/src/G4Scintillation.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4Scintillation.cc,v 1.32 2010/06/16 15:34:15 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4Scintillation.cc,v 1.36 2010/11/03 19:23:48 gum Exp $
     28// GEANT4 tag $Name: xrays-V09-03-05 $
    2929//
    3030////////////////////////////////////////////////////////////////////////
    … …  
    3232////////////////////////////////////////////////////////////////////////
    3333//
    34 // File:        G4Scintillation.cc 
     34// File:        G4Scintillation.cc
    3535// Description: RestDiscrete Process - Generation of Scintillation Photons
    3636// Version:     1.0
    37 // Created:     1998-11-07 
     37// Created:     1998-11-07
    3838// Author:      Peter Gumplinger
    39 // Updated:     2010-92-22 by Peter Gumplinger
     39// Updated:     2010-10-20 Allow the scintillation yield to be a function
     40//              of energy deposited by particle type
     41//              Thanks to Zach Hartwig (Department of Nuclear
     42//              Science and Engineeering - MIT)
     43//              2010-09-22 by Peter Gumplinger
    4044//              > scintillation rise time included, thanks to
    4145//              > Martin Goettlich/DESY
    … …  
    5963//              > change: aSecondaryPosition=x0+rand*aStep.GetDeltaPosition();
    6064//                        aSecondaryTrack->SetTouchable(0);
    61 //              2001-09-17, migration of Materials to pure STL (mma) 
     65//              2001-09-17, migration of Materials to pure STL (mma)
    6266//              2003-06-03, V.Ivanchenko fix compilation warnings
    6367//
    … …  
    6771
    6872#include "G4ios.hh"
     73#include "G4ParticleTypes.hh"
    6974#include "G4EmProcessSubType.hh"
    7075
    … …  
    7277
    7378/////////////////////////
    74 // Class Implementation 
     79// Class Implementation
    7580/////////////////////////
    7681
    … …  
    9398        SetProcessSubType(fScintillation);
    9499
    95         fTrackSecondariesFirst = false;
     100        fTrackSecondariesFirst = false;
    96101        fFiniteRiseTime = false;
    97102
    … …  
    99104        ExcitationRatio = 1.0;
    100105
     106        scintillationByParticleType = false;
     107
    101108        theFastIntegralTable = NULL;
    102109        theSlowIntegralTable = NULL;
    103110
    104         if (verboseLevel>0) {
     111        if (verboseLevel>0) {
    105112           G4cout << GetProcessName() << " is created " << G4endl;
    106         }
    107 
    108         BuildThePhysicsTable();
     113        }
     114
     115        BuildThePhysicsTable();
    109116
    110117        emSaturation = NULL;
    … …  
    115122        ////////////////
    116123
    117 G4Scintillation::~G4Scintillation() 
     124G4Scintillation::~G4Scintillation()
    118125{
    119126        if (theFastIntegralTable != NULL) {
    120            theFastIntegralTable->clearAndDestroy();
     127           theFastIntegralTable->clearAndDestroy();
    121128           delete theFastIntegralTable;
    122         }
     129        }
    123130        if (theSlowIntegralTable != NULL) {
    124131           theSlowIntegralTable->clearAndDestroy();
    … …  
    161168        const G4Material* aMaterial = aTrack.GetMaterial();
    162169
    163         G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
    164         G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
    165 
    166         G4ThreeVector x0 = pPreStepPoint->GetPosition();
     170        G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
     171        G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
     172
     173        G4ThreeVector x0 = pPreStepPoint->GetPosition();
    167174        G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
    168         G4double      t0 = pPreStepPoint->GetGlobalTime();
     175        G4double      t0 = pPreStepPoint->GetGlobalTime();
    169176
    170177        G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
    … …  
    175182             return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
    176183
    177         const G4MaterialPropertyVector* Fast_Intensity =
     184        const G4MaterialPropertyVector* Fast_Intensity =
    178185                aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
    179186        const G4MaterialPropertyVector* Slow_Intensity =
    … …  
    186193        if (Fast_Intensity && Slow_Intensity) nscnt = 2;
    187194
    188         G4double ScintillationYield = aMaterialPropertiesTable->
     195        G4double ScintillationYield = 0.;
     196
     197        if (scintillationByParticleType) {
     198           // The scintillation response is a function of the energy
     199           // deposited by particle types.
     200
     201           // Get the definition of the current particle
     202           G4ParticleDefinition *pDef = aParticle->GetDefinition();
     203           const G4MaterialPropertyVector *Scint_Yield_Vector = NULL;
     204
     205           // Obtain the G4MaterialPropertyVectory containing the
     206           // scintillation light yield as a function of the deposited
     207           // energy for the current particle type
     208
     209           // Protons
     210           if(pDef==G4Proton::ProtonDefinition())
     211             Scint_Yield_Vector = aMaterialPropertiesTable->
     212               GetProperty("PROTONSCINTILLATIONYIELD");
     213
     214           // Deuterons
     215           else if(pDef==G4Deuteron::DeuteronDefinition())
     216             Scint_Yield_Vector = aMaterialPropertiesTable->
     217               GetProperty("DEUTERONSCINTILLATIONYIELD");
     218
     219           // Tritons
     220           else if(pDef==G4Triton::TritonDefinition())
     221             Scint_Yield_Vector = aMaterialPropertiesTable->
     222               GetProperty("TRITONSCINTILLATIONYIELD");
     223
     224           // Alphas
     225           else if(pDef==G4Alpha::AlphaDefinition())
     226             Scint_Yield_Vector = aMaterialPropertiesTable->
     227               GetProperty("ALPHASCINTILLATIONYIELD");
     228         
     229           // Ions (particles derived from G4VIon and G4Ions)
     230           // and recoil ions below tracking cut from neutrons after hElastic
     231           else if(pDef->GetParticleType()== "nucleus" ||
     232                   pDef==G4Neutron::NeutronDefinition())
     233             Scint_Yield_Vector = aMaterialPropertiesTable->
     234               GetProperty("IONSCINTILLATIONYIELD");
     235
     236           // Electrons (must also account for shell-binding energy
     237           // attributed to gamma from standard PhotoElectricEffect)
     238           else if(pDef==G4Electron::ElectronDefinition() ||
     239                   pDef==G4Gamma::GammaDefinition())
     240             Scint_Yield_Vector = aMaterialPropertiesTable->
     241               GetProperty("ELECTRONSCINTILLATIONYIELD");
     242
     243           // Default for particles not enumerated/listed above
     244           else
     245             Scint_Yield_Vector = aMaterialPropertiesTable->
     246               GetProperty("ELECTRONSCINTILLATIONYIELD");
     247           
     248           // If the user has not specified yields for (p,d,t,a,carbon)
     249           // then these unspecified particles will default to the
     250           // electron's scintillation yield
     251           if(!Scint_Yield_Vector){
     252             Scint_Yield_Vector = aMaterialPropertiesTable->
     253               GetProperty("ELECTRONSCINTILLATIONYIELD");
     254           }
     255
     256           // Throw an exception if no scintillation yield is found
     257           if (!Scint_Yield_Vector) {
     258              G4cerr << "\nG4Scintillation::PostStepDoIt(): "
     259                     << "Request for scintillation yield for energy deposit and particle type without correct entry in MaterialPropertiesTable\n"
     260                     << "ScintillationByParticleType requires at minimum that ELECTRONSCINTILLATIONYIELD is set by the user\n"
     261                     << G4endl;
     262             G4Exception("G4Scintillation::PostStepDoIt",
     263                         "No correct entry in MaterialPropertiesTable",
     264                         FatalException,"Missing MaterialPropertiesTable entry.");
     265           }
     266
     267           if (verboseLevel>1) {
     268             G4cout << "\n"
     269                    << "Particle = " << pDef->GetParticleName() << "\n"
     270                    << "Energy Dep. = " << TotalEnergyDeposit/MeV << "\n"
     271                    << "Yield = "
     272                    << Scint_Yield_Vector->GetProperty(TotalEnergyDeposit)
     273                    << "\n" << G4endl;
     274           }
     275
     276           // Obtain the scintillation yield using the total energy
     277           // deposited by the particle in this step.
     278
     279           // Units: [# scintillation photons]
     280           ScintillationYield = Scint_Yield_Vector->
     281                                            GetProperty(TotalEnergyDeposit);
     282        } else {
     283           // The default linear scintillation process
     284           G4double ScintillationYield = aMaterialPropertiesTable->
    189285                                      GetConstProperty("SCINTILLATIONYIELD");
    190         ScintillationYield *= YieldFactor;
     286
     287           // Units: [# scintillation photons / MeV]
     288           ScintillationYield *= YieldFactor;
     289        }
    191290
    192291        G4double ResolutionScale    = aMaterialPropertiesTable->
    … …  
    201300        G4double MeanNumberOfPhotons;
    202301
    203         if (emSaturation) {
     302        // Birk's correction via emSaturation and specifying scintillation by
     303        // by particle type are physically mutually exclusive
     304
     305        if (scintillationByParticleType)
     306           MeanNumberOfPhotons = ScintillationYield;
     307        else if (emSaturation)
    204308           MeanNumberOfPhotons = ScintillationYield*
    205309                              (emSaturation->VisibleEnergyDeposition(&aStep));
    206         } else {
     310        else
    207311           MeanNumberOfPhotons = ScintillationYield*TotalEnergyDeposit;
    208         }
    209312
    210313        G4int NumPhotons;
    … …  
    220323        }
    221324
    222         if (NumPhotons <= 0)
     325        if (NumPhotons <= 0)
    223326        {
    224            // return unchanged particle and no secondaries
    225 
    226            aParticleChange.SetNumberOfSecondaries(0);
     327           // return unchanged particle and no secondaries
     328
     329           aParticleChange.SetNumberOfSecondaries(0);
    227330
    228331           return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
    229         }
    230 
    231         ////////////////////////////////////////////////////////////////
    232 
    233         aParticleChange.SetNumberOfSecondaries(NumPhotons);
    234 
    235         if (fTrackSecondariesFirst) {
     332        }
     333
     334        ////////////////////////////////////////////////////////////////
     335
     336        aParticleChange.SetNumberOfSecondaries(NumPhotons);
     337
     338        if (fTrackSecondariesFirst) {
    236339           if (aTrack.GetTrackStatus() == fAlive )
    237                    aParticleChange.ProposeTrackStatus(fSuspend);
    238         }
    239        
    240         ////////////////////////////////////////////////////////////////
    241 
    242         G4int materialIndex = aMaterial->GetIndex();
    243 
    244         // Retrieve the Scintillation Integral for this material 
    245         // new G4PhysicsOrderedFreeVector allocated to hold CII's
     340                  aParticleChange.ProposeTrackStatus(fSuspend);
     341        }
     342
     343        ////////////////////////////////////////////////////////////////
     344
     345        G4int materialIndex = aMaterial->GetIndex();
     346
     347        // Retrieve the Scintillation Integral for this material 
     348        // new G4PhysicsOrderedFreeVector allocated to hold CII's
    246349
    247350        G4int Num = NumPhotons;
    … …  
    270373                   if (fFiniteRiseTime) {
    271374                      ScintillationRiseTime = aMaterialPropertiesTable->
    272                                   GetConstProperty("SLOWSCINTILLATIONRISETIME");                   }
     375                                  GetConstProperty("SLOWSCINTILLATIONRISETIME");
     376                   }
    273377                   ScintillationIntegral =
    274378                   (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
    … …  
    310414            // Max Scintillation Integral
    311415 
    312             G4double CIImax = ScintillationIntegral->GetMaxValue();
    313                
    314             for (G4int i = 0; i < Num; i++) {
    315 
    316                 // Determine photon energy
     416            G4double CIImax = ScintillationIntegral->GetMaxValue();
     417
     418            for (G4int i = 0; i < Num; i++) {
     419
     420                // Determine photon energy
    317421
    318422                G4double CIIvalue = G4UniformRand()*CIImax;
    319                 G4double sampledEnergy =
     423                G4double sampledEnergy =
    320424                              ScintillationIntegral->GetEnergy(CIIvalue);
    321425
    322                 if (verboseLevel>1) {
     426                if (verboseLevel>1) {
    323427                   G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
    324                    G4cout << "CIIvalue =        " << CIIvalue << G4endl;
    325                 }
    326 
    327                 // Generate random photon direction
     428                   G4cout << "CIIvalue =        " << CIIvalue << G4endl;
     429                }
     430
     431                // Generate random photon direction
    328432
    329433                G4double cost = 1. - 2.*G4UniformRand();
    330434                G4double sint = std::sqrt((1.-cost)*(1.+cost));
    331435
    332                 G4double phi = twopi*G4UniformRand();
    333                 G4double sinp = std::sin(phi);
    334                 G4double cosp = std::cos(phi);
    335 
    336                 G4double px = sint*cosp;
    337                 G4double py = sint*sinp;
    338                 G4double pz = cost;
    339 
    340                 // Create photon momentum direction vector
    341 
    342                 G4ParticleMomentum photonMomentum(px, py, pz);
    343 
    344                 // Determine polarization of new photon
    345 
    346                 G4double sx = cost*cosp;
    347                 G4double sy = cost*sinp;
    348                 G4double sz = -sint;
    349 
    350                 G4ThreeVector photonPolarization(sx, sy, sz);
     436                G4double phi = twopi*G4UniformRand();
     437                G4double sinp = std::sin(phi);
     438                G4double cosp = std::cos(phi);
     439
     440                G4double px = sint*cosp;
     441                G4double py = sint*sinp;
     442                G4double pz = cost;
     443
     444                // Create photon momentum direction vector
     445
     446                G4ParticleMomentum photonMomentum(px, py, pz);
     447
     448                // Determine polarization of new photon
     449
     450                G4double sx = cost*cosp;
     451                G4double sy = cost*sinp;
     452                G4double sz = -sint;
     453
     454                G4ThreeVector photonPolarization(sx, sy, sz);
    351455
    352456                G4ThreeVector perp = photonMomentum.cross(photonPolarization);
    353457
    354                 phi = twopi*G4UniformRand();
    355                 sinp = std::sin(phi);
    356                 cosp = std::cos(phi);
     458                phi = twopi*G4UniformRand();
     459                sinp = std::sin(phi);
     460                cosp = std::cos(phi);
    357461
    358462                photonPolarization = cosp * photonPolarization + sinp * perp;
    … …  
    364468                G4DynamicParticle* aScintillationPhoton =
    365469                  new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(),
    366                                                          photonMomentum);
    367                 aScintillationPhoton->SetPolarization
    368                                      (photonPolarization.x(),
    369                                       photonPolarization.y(),
    370                                       photonPolarization.z());
    371 
    372                 aScintillationPhoton->SetKineticEnergy(sampledEnergy);
     470                                                         photonMomentum);
     471                aScintillationPhoton->SetPolarization
     472                                     (photonPolarization.x(),
     473                                      photonPolarization.y(),
     474                                      photonPolarization.z());
     475
     476                aScintillationPhoton->SetKineticEnergy(sampledEnergy);
    373477
    374478                // Generate new G4Track object:
    … …  
    401505                                    x0 + rand * aStep.GetDeltaPosition();
    402506
    403                 G4Track* aSecondaryTrack =
    404                 new G4Track(aScintillationPhoton,aSecondaryTime,aSecondaryPosition);
     507                G4Track* aSecondaryTrack =
     508                new G4Track(aScintillationPhoton,aSecondaryTime,aSecondaryPosition);
    405509
    406510                aSecondaryTrack->SetTouchableHandle(
    … …  
    410514                aSecondaryTrack->SetParentID(aTrack.GetTrackID());
    411515
    412                 aParticleChange.AddSecondary(aSecondaryTrack);
    413 
    414             }
    415         }
    416 
    417         if (verboseLevel>0) {
    418         G4cout << "\n Exiting from G4Scintillation::DoIt -- NumberOfSecondaries = "
    419              << aParticleChange.GetNumberOfSecondaries() << G4endl;
    420         }
    421 
    422         return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
     516                aParticleChange.AddSecondary(aSecondaryTrack);
     517
     518            }
     519        }
     520
     521        if (verboseLevel>0) {
     522        G4cout << "\n Exiting from G4Scintillation::DoIt -- NumberOfSecondaries = "
     523               << aParticleChange.GetNumberOfSecondaries() << G4endl;
     524        }
     525
     526        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
    423527}
    424528
    … …  
    429533void G4Scintillation::BuildThePhysicsTable()
    430534{
    431         if (theFastIntegralTable && theSlowIntegralTable) return;
    432 
    433         const G4MaterialTable* theMaterialTable =
     535        if (theFastIntegralTable && theSlowIntegralTable) return;
     536
     537        const G4MaterialTable* theMaterialTable =
    434538                               G4Material::GetMaterialTable();
    435         G4int numOfMaterials = G4Material::GetNumberOfMaterials();
    436 
    437         // create new physics table
     539        G4int numOfMaterials = G4Material::GetNumberOfMaterials();
     540
     541        // create new physics table
    438542       
    439         if(!theFastIntegralTable)theFastIntegralTable = new G4PhysicsTable(numOfMaterials);
     543        if(!theFastIntegralTable)theFastIntegralTable = new G4PhysicsTable(numOfMaterials);
    440544        if(!theSlowIntegralTable)theSlowIntegralTable = new G4PhysicsTable(numOfMaterials);
    441545
    442         // loop for materials
    443 
    444         for (G4int i=0 ; i < numOfMaterials; i++)
    445         {
    446                 G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
     546        // loop for materials
     547
     548        for (G4int i=0 ; i < numOfMaterials; i++)
     549        {
     550                G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
    447551                                        new G4PhysicsOrderedFreeVector();
    448552                G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector =
    449553                                        new G4PhysicsOrderedFreeVector();
    450554
    451                 // Retrieve vector of scintillation wavelength intensity for
     555                // Retrieve vector of scintillation wavelength intensity for
    452556                // the material from the material's optical properties table.
    453557
    454                 G4Material* aMaterial = (*theMaterialTable)[i];
    455 
    456                 G4MaterialPropertiesTable* aMaterialPropertiesTable =
    457                                 aMaterial->GetMaterialPropertiesTable();
    458 
    459                 if (aMaterialPropertiesTable) {
    460 
    461                    G4MaterialPropertyVector* theFastLightVector =
    462                    aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
    463 
    464                    if (theFastLightVector) {
    465                
    466                       // Retrieve the first intensity point in vector
    467                       // of (photon energy, intensity) pairs
    468 
    469                       theFastLightVector->ResetIterator();
    470                       ++(*theFastLightVector);  // advance to 1st entry
    471 
    472                       G4double currentIN = theFastLightVector->
    473                                            GetProperty();
    474 
    475                       if (currentIN >= 0.0) {
    476 
    477                         // Create first (photon energy, Scintillation
     558                G4Material* aMaterial = (*theMaterialTable)[i];
     559
     560                G4MaterialPropertiesTable* aMaterialPropertiesTable =
     561                                aMaterial->GetMaterialPropertiesTable();
     562
     563                if (aMaterialPropertiesTable) {
     564
     565                   G4MaterialPropertyVector* theFastLightVector =
     566                   aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
     567
     568                   if (theFastLightVector) {
     569
     570                      // Retrieve the first intensity point in vector
     571                      // of (photon energy, intensity) pairs
     572
     573                      theFastLightVector->ResetIterator();
     574                      ++(*theFastLightVector);  // advance to 1st entry
     575
     576                      G4double currentIN = theFastLightVector->
     577                                               GetProperty();
     578
     579                      if (currentIN >= 0.0) {
     580
     581                        // Create first (photon energy, Scintillation
    478582                         // Integral pair 
    479583
    480                         G4double currentPM = theFastLightVector->
    481                                                 GetPhotonEnergy();
    482 
    483                         G4double currentCII = 0.0;
    484 
    485                         aPhysicsOrderedFreeVector->
    486                                 InsertValues(currentPM , currentCII);
    487 
    488                         // Set previous values to current ones prior to loop
    489 
    490                         G4double prevPM  = currentPM;
    491                         G4double prevCII = currentCII;
    492                         G4double prevIN  = currentIN;
    493 
    494                         // loop over all (photon energy, intensity)
    495                         // pairs stored for this material 
    496 
    497                         while(++(*theFastLightVector))
    498                         {
    499                                 currentPM = theFastLightVector->
    500                                                 GetPhotonEnergy();
    501 
    502                                 currentIN=theFastLightVector-> 
    503                                                 GetProperty();
    504 
    505                                 currentCII = 0.5 * (prevIN + currentIN);
    506 
    507                                 currentCII = prevCII +
    508                                              (currentPM - prevPM) * currentCII;
    509 
    510                                 aPhysicsOrderedFreeVector->
    511                                     InsertValues(currentPM, currentCII);
    512 
    513                                 prevPM  = currentPM;
    514                                 prevCII = currentCII;
    515                                 prevIN  = currentIN;
    516                         }
    517 
    518                       }
    519                    }
     584                        G4double currentPM = theFastLightVector->
     585                                                  GetPhotonEnergy();
     586
     587                        G4double currentCII = 0.0;
     588
     589                        aPhysicsOrderedFreeVector->
     590                                InsertValues(currentPM , currentCII);
     591
     592                        // Set previous values to current ones prior to loop
     593
     594                        G4double prevPM  = currentPM;
     595                        G4double prevCII = currentCII;
     596                        G4double prevIN  = currentIN;
     597
     598                        // loop over all (photon energy, intensity)
     599                        // pairs stored for this material 
     600
     601                        while(++(*theFastLightVector))
     602                        {
     603                                currentPM = theFastLightVector->
     604                                                GetPhotonEnergy();
     605
     606                                currentIN = theFastLightVector->       
     607                                                GetProperty();
     608
     609                                currentCII = 0.5 * (prevIN + currentIN);
     610
     611                                currentCII = prevCII +
     612                                             (currentPM - prevPM) * currentCII;
     613
     614                                aPhysicsOrderedFreeVector->
     615                                    InsertValues(currentPM, currentCII);
     616
     617                                prevPM  = currentPM;
     618                                prevCII = currentCII;
     619                                prevIN  = currentIN;
     620                        }
     621
     622                      }
     623                   }
    520624
    521625                   G4MaterialPropertyVector* theSlowLightVector =
    … …  
    578682                      }
    579683                   }
    580                 }
    581 
    582         // The scintillation integral(s) for a given material
    583         // will be inserted in the table(s) according to the
    584         // position of the material in the material table.
    585 
    586         theFastIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
     684                }
     685
     686        // The scintillation integral(s) for a given material
     687        // will be inserted in the table(s) according to the
     688        // position of the material in the material table.
     689
     690        theFastIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
    587691        theSlowIntegralTable->insertAt(i,bPhysicsOrderedFreeVector);
    588692
    589         }
     693        }
     694}
     695
     696// Called by the user to set the scintillation yield as a function
     697// of energy deposited by particle type
     698
     699void G4Scintillation::SetScintillationByParticleType(const G4bool scintType)
     700{
     701        if (emSaturation) {
     702           G4Exception("G4Scintillation::SetScintillationByParticleType", "Redefinition",
     703                       JustWarning, "Birks Saturation is replaced by ScintillationByParticleType!");
     704           RemoveSaturation();
     705        }
     706        scintillationByParticleType = scintType;
    590707}
    591708
    … …  
    600717        *condition = StronglyForced;
    601718
    602         return DBL_MAX;
     719        return DBL_MAX;
    603720
    604721}

  • trunk/source/processes/electromagnetic/xrays/src/G4SynchrotronRadiation.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4SynchrotronRadiation.cc,v 1.6 2010/06/16 15:34:15 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4SynchrotronRadiation.cc,v 1.8 2010/10/14 18:38:21 vnivanch Exp $
     28// GEANT4 tag $Name: xrays-V09-03-05 $
    2929//
    3030// --------------------------------------------------------------
    … …  
    4545#include "G4SynchrotronRadiation.hh"
    4646#include "G4UnitsTable.hh"
     47#include "G4EmProcessSubType.hh"
    4748
    4849///////////////////////////////////////////////////////////////////////
    … …  
    6566                           (2.5*fine_structure_const*eplus*c_light) ;
    6667  fEnergyConst = 1.5*c_light*c_light*eplus*hbar_Planck/electron_mass_c2  ;
     68
     69  SetProcessSubType(fSynchrotronRadiation);
    6770  verboseLevel=1;
    6871}
    … …  
    7477
    7578G4SynchrotronRadiation::~G4SynchrotronRadiation()
    76 {
    77      ;
    78 }
     79{}
    7980
    8081/////////////////////////////// METHODS /////////////////////////////////

  • trunk/source/processes/electromagnetic/xrays/src/G4SynchrotronRadiationInMat.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4SynchrotronRadiationInMat.cc,v 1.3 2010/06/16 15:34:15 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4SynchrotronRadiationInMat.cc,v 1.5 2010/10/14 18:38:21 vnivanch Exp $
     28// GEANT4 tag $Name: xrays-V09-03-05 $
    2929//
    3030// --------------------------------------------------------------
    … …  
    4343#include "G4SynchrotronRadiationInMat.hh"
    4444#include "G4Integrator.hh"
     45#include "G4EmProcessSubType.hh"
    4546
    4647////////////////////////////////////////////////////////////////////
    … …  
    132133
    133134  fFieldPropagator = transportMgr->GetPropagatorInField();
     135  SetProcessSubType(fSynchrotronRadiation);
    134136
    135137}
    … …  
    141143 
    142144G4SynchrotronRadiationInMat::~G4SynchrotronRadiationInMat()
    143 {
    144      ;
    145 }
     145{}
    146146 
    147147 

  • trunk/source/processes/electromagnetic/xrays/src/G4TransitionRadiation.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4TransitionRadiation.cc,v 1.8 2010/06/16 15:34:15 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4TransitionRadiation.cc,v 1.10 2010/10/14 18:38:21 vnivanch Exp $
     28// GEANT4 tag $Name: xrays-V09-03-05 $
    2929//
    3030// G4TransitionRadiation class -- implementation file
    … …  
    4545#include "G4TransitionRadiation.hh"
    4646#include "G4Material.hh"
     47#include "G4EmProcessSubType.hh"
    4748
    4849// Local constants
    … …  
    6263  : G4VDiscreteProcess(processName, type)
    6364{
     65  SetProcessSubType(fTransitionRadiation);
    6466  //  fMatIndex1 = pMat1->GetIndex() ;
    6567  //  fMatIndex2 = pMat2->GetIndex() ;
    … …  
    7274
    7375G4TransitionRadiation::~G4TransitionRadiation()
    74 {
    75         ;
    76 }
     76{}
    7777
    7878
Note: See TracChangeset for help on using the changeset viewer.