Ignore:
Timestamp:
Dec 22, 2010, 3:52:27 PM (14 years ago)
Author:
garnier
Message:

geant4 tag 9.4

Location:
trunk/source/processes/hadronic/models/de_excitation/evaporation
Files:
38 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4AlphaEvaporationChannel.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4AlphaEvaporationChannel.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4AlphaEvaporationChannel.hh,v 1.9 2010/11/17 12:14:59 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Nov. 1999)
    3232//
     33// 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup
    3334
    3435#ifndef G4AlphaEvaporationChannel_h
    … …  
    4344public:
    4445  // only available constructor
    45   G4AlphaEvaporationChannel() : G4EvaporationChannel(4,2,"alpha",
    46                                                      &theEvaporationProbability,&theCoulombBarrier) {};
     46  G4AlphaEvaporationChannel();
    4747
    4848  // destructor
    49   ~G4AlphaEvaporationChannel() {};
     49  virtual ~G4AlphaEvaporationChannel();
    5050
    5151private:
    … …  
    5454  G4AlphaEvaporationChannel(const G4AlphaEvaporationChannel & right);
    5555
    56 
    57 public:
    5856  G4bool operator==(const G4AlphaEvaporationChannel & right) const;
    5957  G4bool operator!=(const G4AlphaEvaporationChannel & right) const;
    60 
    61 private:
    6258
    6359  G4AlphaCoulombBarrier theCoulombBarrier;     

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4AlphaEvaporationProbability.hh

    r962 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4AlphaEvaporationProbability.hh,v 1.14 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
     33//
     34// Modified:
     35// 17-11-2010 V.Ivanchenko integer Z and A
    3036//
    3137#ifndef G4AlphaEvaporationProbability_h
    … …  
    4248  G4AlphaEvaporationProbability();
    4349
    44   ~G4AlphaEvaporationProbability() {}
     50  virtual ~G4AlphaEvaporationProbability();
     51
    4552private: 
    4653  // Copy constructor
    … …  
    5158  G4bool operator!=(const G4AlphaEvaporationProbability &right) const;
    5259
    53 
    5460private:
    5561
    56   virtual G4double CrossSection(const  G4Fragment & fragment, const  G4double K);
     62  virtual G4double CrossSection(const  G4Fragment & fragment, G4double K);
    5763
    58   G4double GetOpt0(const G4double K);
    59   G4double GetOpt12(const G4double K);
    60   G4double GetOpt34(const G4double K);
     64  G4double GetOpt0(G4double K);
     65  G4double GetOpt12(G4double K);
     66  G4double GetOpt34(G4double K);
     67 
     68  virtual G4double CalcAlphaParam(const G4Fragment & fragment) ;
     69 
     70  virtual G4double CalcBetaParam(const G4Fragment & fragment)  ;
     71 
     72  G4double CCoeficient(G4int aZ) ;
     73 
     74  //data members
     75   
     76  G4AlphaCoulombBarrier theCoulombBarrier;
    6177
    62  
    63  virtual G4double CalcAlphaParam(const G4Fragment & fragment) ;
    64  
    65  virtual G4double CalcBetaParam(const G4Fragment & fragment)  ;
    66  
    67   G4double CCoeficient(const G4double aZ) ;
    68  
    69 //data members
    70    
    71       G4AlphaCoulombBarrier theCoulombBarrier;
    72 
    73       G4double ResidualA;
    74       G4double ResidualZ;
    75       G4double theA;
    76       G4double theZ;
    77       G4double ResidualAthrd;
    78       G4double FragmentA;
    79       G4double FragmentAthrd;
    80 
     78  G4int ResidualA;
     79  G4int ResidualZ;
     80  G4int theA;
     81  G4int theZ;
     82  G4double ResidualAthrd;
     83  G4int FragmentA;
     84  G4double FragmentAthrd;
    8185
    8286};

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4DeuteronEvaporationChannel.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4DeuteronEvaporationChannel.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4DeuteronEvaporationChannel.hh,v 1.9 2010/11/17 12:14:59 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Nov. 1999)
    3232//
     33// 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup
    3334
    3435#ifndef G4DeuteronEvaporationChannel_h
    … …  
    4344public:
    4445  // only available constructor
    45   G4DeuteronEvaporationChannel() : G4EvaporationChannel(2,1,"deuteron",
    46                                                         &theEvaporationProbability,&theCoulombBarrier) {};
     46  G4DeuteronEvaporationChannel();
    4747
    4848  // destructor
    49   ~G4DeuteronEvaporationChannel() {};
     49  virtual ~G4DeuteronEvaporationChannel();
    5050
    5151private:
    … …  
    5454  G4DeuteronEvaporationChannel(const G4DeuteronEvaporationChannel & right);
    5555
    56 public:
    5756  G4bool operator==(const G4DeuteronEvaporationChannel & right) const;
    5857  G4bool operator!=(const G4DeuteronEvaporationChannel & right) const;
    59 
    60 private:
    6158
    6259  G4DeuteronCoulombBarrier theCoulombBarrier;

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4DeuteronEvaporationProbability.hh

    r962 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4DeuteronEvaporationProbability.hh,v 1.14 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
     33//
     34// Modified:
     35// 17-11-2010 V.Ivanchenko integer Z and A
    3036//
    3137
    … …  
    4349  G4DeuteronEvaporationProbability();
    4450
    45   ~G4DeuteronEvaporationProbability() {}
     51  virtual ~G4DeuteronEvaporationProbability();
     52
    4653private: 
    4754  // Copy constructor
    … …  
    5562private:
    5663
    57   virtual G4double CrossSection(const  G4Fragment & fragment, const  G4double K);
     64  virtual G4double CrossSection(const  G4Fragment & fragment, G4double K);
    5865
    59   G4double GetOpt0(const G4double K);
    60   G4double GetOpt12(const G4double K);
    61   G4double GetOpt34(const G4double K);
     66  G4double GetOpt0(G4double K);
     67  G4double GetOpt12(G4double K);
     68  G4double GetOpt34(G4double K);
    6269
     70  virtual G4double CalcAlphaParam(const G4Fragment & fragment) ;
     71 
     72  virtual G4double CalcBetaParam(const G4Fragment & fragment) ;
     73 
     74  G4double CCoeficient(G4int aZ) ;
    6375 
    64  virtual G4double CalcAlphaParam(const G4Fragment & fragment) ;
    65  
    66  virtual G4double CalcBetaParam(const G4Fragment & fragment) ;
    67  
    68   G4double CCoeficient(const G4double aZ) ;
    69  
    70 //data members
     76  //data members
    7177
    72       G4DeuteronCoulombBarrier theCoulombBarrier;
     78  G4DeuteronCoulombBarrier theCoulombBarrier;
    7379
    74       G4double ResidualA;
    75       G4double ResidualZ;
    76       G4double theA;
    77       G4double theZ;
    78       G4double ResidualAthrd;
    79       G4double FragmentA;
    80       G4double FragmentAthrd;
    81 
    82 
     80  G4int ResidualA;
     81  G4int ResidualZ;
     82  G4int theA;
     83  G4int theZ;
     84  G4double ResidualAthrd;
     85  G4int FragmentA;
     86  G4double FragmentAthrd;
    8387};
    8488

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4Evaporation.hh

    r1340 r1347  
    2626//
    2727// $Id: G4Evaporation.hh,v 1.12 2010/05/11 11:34:09 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4EvaporationChannel.hh

    r962 r1347  
    2424// ********************************************************************
    2525//
     26// $Id: G4EvaporationChannel.hh,v 1.11 2010/11/17 12:19:08 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29//
    2630//J.M. Quesada (August2008). Based on:
    2731//
    … …  
    2933// by V. Lara (Oct 1998)
    3034//
    31 
     35// 17-11-2010 V.Ivanchenko in constructor replace G4VEmissionProbability by
     36//            G4EvaporationProbability and do not new and delete probability
     37//            object at each call; use G4Pow
    3238
    3339#ifndef G4EvaporationChannel_h
    … …  
    3541
    3642#include "G4VEvaporationChannel.hh"
    37 #include "G4VEmissionProbability.hh"
    3843#include "G4EvaporationProbability.hh"
    39 #include "G4NeutronEvaporationProbability.hh"
    40 #include "G4ProtonEvaporationProbability.hh"
    41 #include "G4DeuteronEvaporationProbability.hh"
    42 #include "G4TritonEvaporationProbability.hh"
    43 #include "G4He3EvaporationProbability.hh"
    44 #include "G4AlphaEvaporationProbability.hh"
    45 #include "G4VLevelDensityParameter.hh"
    4644#include "G4VCoulombBarrier.hh"
    47 #include "G4EvaporationLevelDensityParameter.hh"
    48 #include "G4NucleiProperties.hh"
    49 #include "Randomize.hh"
    50 #include "G4ParticleTable.hh"
    51 #include "G4IonTable.hh"
    5245
     46class G4EvaporationLevelDensityParameter;
    5347
    5448class G4EvaporationChannel : public G4VEvaporationChannel
    … …  
    5650public:
    5751  // constructor
    58 
    59 
    60   G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String & aName,
    61                        G4VEmissionProbability * aEmissionStrategy,
     52  G4EvaporationChannel(G4int theA, G4int theZ, const G4String & aName,
     53                       G4EvaporationProbability * aEmissionStrategy,
    6254                       G4VCoulombBarrier * aCoulombBarrier);
    6355public:
    6456  // destructor
    65   ~G4EvaporationChannel();
     57  virtual ~G4EvaporationChannel();
    6658 
    67   void SetEmissionStrategy(G4VEmissionProbability * aEmissionStrategy)
     59  inline void SetEmissionStrategy(G4EvaporationProbability * aEmissionStrategy)
    6860  {theEvaporationProbabilityPtr = aEmissionStrategy;}
    6961
    70   void SetCoulombBarrierStrategy(G4VCoulombBarrier * aCoulombBarrier)
     62  inline void SetCoulombBarrierStrategy(G4VCoulombBarrier * aCoulombBarrier)
    7163  {theCoulombBarrierPtr = aCoulombBarrier;}
    72  
    73 
    74  
     64   
    7565protected:
    7666  // default constructor
    77   G4EvaporationChannel() {};
     67  G4EvaporationChannel();
    7868 
    7969private:
    … …  
    9383  G4FragmentVector * BreakUp(const G4Fragment & theNucleus);
    9484
    95   // void SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity);
    96  
    9785public:
    98 
    9986
    10087  inline G4double GetEmissionProbability(void) const
    10188  {return EmissionProbability;}
    102  
    103  
     89   
    10490  inline G4double GetMaximalKineticEnergy(void) const
    10591  { return MaximalKineticEnergy; }
    … …  
    10894 
    10995  // Calculate Binding Energy for separate fragment from nucleus
    110   G4double CalcBindingEnergy(const G4int anA, const G4int aZ);
     96  G4double CalcBindingEnergy(G4int anA, G4int aZ);
    11197
    11298  // Calculate maximal kinetic energy that can be carried by fragment (in MeV)
    113   G4double CalcMaximalKineticEnergy(const G4double U);
     99  G4double CalcMaximalKineticEnergy(G4double U);
    114100
    115101  // Samples fragment kinetic energy.
    116     G4double  GetKineticEnergy(const G4Fragment & aFragment);
     102  G4double  GetKineticEnergy(const G4Fragment & aFragment);
    117103
    118104  // This has to be removed and put in Random Generator
    119   G4ThreeVector IsotropicVector(const G4double Magnitude  = 1.0);
     105  G4ThreeVector IsotropicVector(G4double Magnitude  = 1.0);
    120106
    121107        // Data Members
    … …  
    132118  G4int theZ;
    133119
     120  G4double EvaporatedMass;
     121  G4double ResidualMass;
    134122
    135123  // For evaporation probability calcualation
    136   G4VEmissionProbability * theEvaporationProbabilityPtr;
     124  G4EvaporationProbability * theEvaporationProbabilityPtr;
    137125
    138126  // For Level Density calculation
    139  // G4bool MyOwnLevelDensity;
     127  // G4bool MyOwnLevelDensity;
    140128  G4VLevelDensityParameter * theLevelDensityPtr;
    141 
    142129
    143130  // For Coulomb Barrier calculation
    144131  G4VCoulombBarrier * theCoulombBarrierPtr;
    145132  G4double CoulombBarrier;
    146  
    147  
     133   
    148134  //---------------------------------------------------
    149135
    … …  
    161147  G4double EmissionProbability;
    162148
    163 
    164149  // Maximal Kinetic Energy that can be carried by fragment
    165150  G4double MaximalKineticEnergy;
    166 
    167151
    168152};

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4EvaporationDefaultGEMFactory.hh

    r1340 r1347  
    2626//
    2727// $Id: G4EvaporationDefaultGEMFactory.hh,v 1.1 2009/07/27 10:20:13 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4EvaporationFactory.hh

    r1340 r1347  
    2626//
    2727// $Id: G4EvaporationFactory.hh,v 1.4 2010/04/27 11:43:16 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4EvaporationProbability.hh

    r962 r1347  
    2424// ********************************************************************
    2525//
     26// $Id: G4EvaporationProbability.hh,v 1.13 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
    2629//J.M. Quesada (August2008). Based on:
    2730//
    … …  
    4447public:
    4548  // Only available constructor
    46   G4EvaporationProbability(const G4int anA, const G4int aZ, const G4double aGamma,G4VCoulombBarrier * aCoulombBarrier) :
    47     theA(anA),
    48     theZ(aZ),
    49     Gamma(aGamma)
    50 ,  theCoulombBarrierptr(aCoulombBarrier)
    51   {
    52     theEvapLDPptr = new G4EvaporationLevelDensityParameter;
     49  G4EvaporationProbability(G4int anA, G4int aZ, G4double aGamma,
     50                           G4VCoulombBarrier * aCoulombBarrier);
    5351
    54    
    55   }
     52  virtual ~G4EvaporationProbability();
    5653
    57   ~G4EvaporationProbability()
    58   {
    59     if (theEvapLDPptr != 0) delete theEvapLDPptr;
    60 
    61   }
    62 
    63 
     54  inline G4int GetZ(void) const { return theZ; }
    6455       
    65   G4double GetZ(void) const { return theZ; }
    66        
    67   G4double GetA(void) const { return theA;}
     56  inline G4int GetA(void) const { return theA;}
    6857
    6958protected:
    7059 
    7160  // Default constructor
    72   G4EvaporationProbability() {}
     61  G4EvaporationProbability();
    7362
    7463private:
    … …  
    8271public:
    8372
    84  G4double ProbabilityDistributionFunction( const G4Fragment & aFragment, const G4double K);
     73  G4double ProbabilityDistributionFunction( const G4Fragment & aFragment, G4double K);
    8574
    86  G4double EmissionProbability(const G4Fragment & fragment, const G4double anEnergy);
     75  G4double EmissionProbability(const G4Fragment & fragment, G4double anEnergy);
    8776
    8877private:
    8978
    90   G4double CalculateProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy );
     79  G4double CalculateProbability(const G4Fragment & fragment, G4double MaximalKineticEnergy );
    9180
    92   G4double IntegrateEmissionProbability(const G4Fragment & aFragment, const G4double & Low, const G4double & Up );
     81  G4double IntegrateEmissionProbability(const G4Fragment & aFragment,
     82                                        const G4double & Low, const G4double & Up );
    9383
    9484protected:
    9585
    96  virtual G4double CrossSection( const  G4Fragment & fragment, const G4double K )= 0; 
     86 virtual G4double CrossSection( const  G4Fragment & fragment, G4double K )= 0; 
    9787
    9888 virtual G4double CalcAlphaParam(const G4Fragment & fragment)=0 ;
    … …  
    10292private:
    10393
    104   // Data Members
    105 
    106   G4VLevelDensityParameter * theEvapLDPptr;
    107        
     94  // Data Members       
    10895  G4int theA;
    10996  G4int theZ;
    … …  
    113100  G4double Gamma;
    114101
    115 //The Coulomb Barrier
    116          G4VCoulombBarrier * theCoulombBarrierptr;
    117 
     102  //The Coulomb Barrier
     103  G4VCoulombBarrier * theCoulombBarrierptr;
    118104
    119105};
    120106
    121 
    122 
    123 
    124107#endif

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4He3EvaporationChannel.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4He3EvaporationChannel.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4He3EvaporationChannel.hh,v 1.9 2010/11/17 12:14:59 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Nov. 1999)
    3232//
     33// 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup
    3334
    3435#ifndef G4He3EvaporationChannel_h
    … …  
    4344public:
    4445  // only available constructor
    45   G4He3EvaporationChannel() : G4EvaporationChannel(3,2,"He3",
    46                                                    &theEvaporationProbability,&theCoulombBarrier) {};
     46  G4He3EvaporationChannel();
    4747
    4848  // destructor
    49   ~G4He3EvaporationChannel() {};
     49  virtual ~G4He3EvaporationChannel();
    5050
    5151private:
     52
    5253  const G4He3EvaporationChannel & operator=(const G4He3EvaporationChannel & right); 
    5354
    5455  G4He3EvaporationChannel(const G4He3EvaporationChannel & right);
    5556
    56 public:
    5757  G4bool operator==(const G4He3EvaporationChannel & right) const;
    5858  G4bool operator!=(const G4He3EvaporationChannel & right) const;
    5959
    60 private:
    61 
    62     G4He3CoulombBarrier theCoulombBarrier;
     60  G4He3CoulombBarrier theCoulombBarrier;
    6361       
    6462  G4He3EvaporationProbability theEvaporationProbability;

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4He3EvaporationProbability.hh

    r962 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4He3EvaporationProbability.hh,v 1.14 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
     33//
     34// Modified:
     35// 17-11-2010 V.Ivanchenko integer Z and A
    3036//
    3137#ifndef G4He3EvaporationProbability_h
    … …  
    4248  G4He3EvaporationProbability();
    4349
    44   ~G4He3EvaporationProbability() {}
     50  virtual ~G4He3EvaporationProbability();
     51
    4552private: 
    4653  // Copy constructor
    … …  
    5158  G4bool operator!=(const G4He3EvaporationProbability &right) const;
    5259
    53 
    5460private:
    5561
    56   virtual G4double CrossSection(const  G4Fragment & fragment, const  G4double K);
     62  virtual G4double CrossSection(const  G4Fragment & fragment, G4double K);
    5763
    58   G4double GetOpt0(const G4double K);
    59   G4double GetOpt12(const G4double K);
    60   G4double GetOpt34(const G4double K);
     64  G4double GetOpt0(G4double K);
     65  G4double GetOpt12(G4double K);
     66  G4double GetOpt34(G4double K);
    6167
     68  virtual G4double CalcAlphaParam(const G4Fragment & fragment) ;
     69 
     70  virtual G4double CalcBetaParam(const G4Fragment & fragment) ;
     71 
     72  G4double CCoeficient(G4int aZ) ;
    6273 
    63  virtual G4double CalcAlphaParam(const G4Fragment & fragment) ;
    64  
    65  virtual G4double CalcBetaParam(const G4Fragment & fragment) ;
    66  
    67   G4double CCoeficient(const G4double aZ) ;
    68  
    69 //data members
     74  //data members
    7075   
    71       G4He3CoulombBarrier theCoulombBarrier;
     76  G4He3CoulombBarrier theCoulombBarrier;
    7277
    73       G4double ResidualA;
    74       G4double ResidualZ;
    75       G4double theA;
    76       G4double theZ;
    77       G4double ResidualAthrd;
    78       G4double FragmentA;
    79       G4double FragmentAthrd;
    80 
    81 
     78  G4int ResidualA;
     79  G4int ResidualZ;
     80  G4int theA;
     81  G4int theZ;
     82  G4double ResidualAthrd;
     83  G4int FragmentA;
     84  G4double FragmentAthrd;
    8285};
    8386

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4NeutronEvaporationChannel.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4NeutronEvaporationChannel.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4NeutronEvaporationChannel.hh,v 1.9 2010/11/17 12:14:59 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Nov. 1999)
    3232//
     33// 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup
    3334
    3435#ifndef G4NeutronEvaporationChannel_h
    … …  
    4344public:
    4445  // only available constructor
    45   G4NeutronEvaporationChannel() : G4EvaporationChannel(1,0,"neutron",
    46                                                        &theEvaporationProbability,&theCoulombBarrier) {};
     46  G4NeutronEvaporationChannel();
    4747
    4848  // destructor
    49   ~G4NeutronEvaporationChannel() {};
     49  virtual ~G4NeutronEvaporationChannel();
    5050
    5151private:
     52
    5253  const G4NeutronEvaporationChannel & operator=(const G4NeutronEvaporationChannel & right); 
    5354
    5455  G4NeutronEvaporationChannel(const G4NeutronEvaporationChannel & right);
    5556
    56 public:
    5757  G4bool operator==(const G4NeutronEvaporationChannel & right) const;
    5858  G4bool operator!=(const G4NeutronEvaporationChannel & right) const;
    5959
    60 private:
    61 
    62  G4NeutronCoulombBarrier theCoulombBarrier;
     60  G4NeutronCoulombBarrier theCoulombBarrier;
    6361       
    6462  G4NeutronEvaporationProbability theEvaporationProbability;

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4NeutronEvaporationProbability.hh

    r962 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4NeutronEvaporationProbability.hh,v 1.15 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
     33//
     34// Modified:
     35// 17-11-2010 V.Ivanchenko integer Z and A
    3036//
    3137
    … …  
    4349  G4NeutronEvaporationProbability();
    4450               
    45   ~G4NeutronEvaporationProbability() {}
     51  virtual ~G4NeutronEvaporationProbability();
     52
    4653private: 
    4754 
    … …  
    5461private:
    5562
    56   virtual G4double CrossSection(const  G4Fragment & fragment, const  G4double K);
     63  virtual G4double CrossSection(const  G4Fragment & fragment, G4double K);
    5764
    58   G4double GetOpt12(const G4double K);
    59   G4double GetOpt34(const G4double K);
     65  G4double GetOpt12(G4double K);
     66  G4double GetOpt34(G4double K);
    6067
    61  virtual G4double CalcAlphaParam(const G4Fragment & fragment);
     68  virtual G4double CalcAlphaParam(const G4Fragment & fragment);
    6269 
    63  virtual G4double CalcBetaParam(const G4Fragment & fragment);
     70  virtual G4double CalcBetaParam(const G4Fragment & fragment);
    6471 
    65  
    66 //data members
     72  //data members
    6773
    68       G4NeutronCoulombBarrier theCoulombBarrier;
     74  G4NeutronCoulombBarrier theCoulombBarrier;
    6975
    70       G4double ResidualA;
    71       G4double ResidualZ;
    72       G4double theA;
    73       G4double theZ;
    74       G4double ResidualAthrd;
    75       G4double FragmentA;
    76       G4double FragmentAthrd;
    77 
    78 
    79 
    80 
     76  G4int ResidualA;
     77  G4int ResidualZ;
     78  G4int theA;
     79  G4int theZ;
     80  G4double ResidualAthrd;
     81  G4int FragmentA;
     82  G4double FragmentAthrd;
    8183};
    8284

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4ProtonEvaporationChannel.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4ProtonEvaporationChannel.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4ProtonEvaporationChannel.hh,v 1.9 2010/11/17 12:14:59 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Nov. 1999)
    3232//
     33// 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup
    3334
    3435#ifndef G4ProtonEvaporationChannel_h
    … …  
    4344public:
    4445  // only available constructor
    45   G4ProtonEvaporationChannel() : G4EvaporationChannel(1,1,"proton",
    46                                                       &theEvaporationProbability,&theCoulombBarrier) {};
     46  G4ProtonEvaporationChannel();
    4747
    4848  // destructor
    49   ~G4ProtonEvaporationChannel() {};
     49  virtual ~G4ProtonEvaporationChannel();
    5050
    5151private:
     52
    5253  const G4ProtonEvaporationChannel & operator=(const G4ProtonEvaporationChannel & right); 
    5354
    5455  G4ProtonEvaporationChannel(const G4ProtonEvaporationChannel & right);
    5556
    56 public:
    5757  G4bool operator==(const G4ProtonEvaporationChannel & right) const;
    5858  G4bool operator!=(const G4ProtonEvaporationChannel & right) const;
    5959
    60 
    61 private:
    62 
    63    G4ProtonCoulombBarrier  theCoulombBarrier;
     60  G4ProtonCoulombBarrier  theCoulombBarrier;
    6461 
    6562  G4ProtonEvaporationProbability  theEvaporationProbability;

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4ProtonEvaporationProbability.hh

    r962 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4ProtonEvaporationProbability.hh,v 1.14 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
     33//
     34// Modified:
     35// 17-11-2010 V.Ivanchenko integer Z and A
    3036//
    3137
    … …  
    4349  G4ProtonEvaporationProbability();
    4450
    45   ~G4ProtonEvaporationProbability() {}
     51  virtual ~G4ProtonEvaporationProbability();
     52
    4653private: 
    4754  // Copy constructor
    … …  
    5259  G4bool operator!=(const G4ProtonEvaporationProbability &right) const;
    5360
    54 
    5561private:
    5662
    57   virtual G4double CrossSection(const  G4Fragment & fragment, const  G4double K);
     63  virtual G4double CrossSection(const  G4Fragment & fragment, G4double K);
    5864
    59   G4double GetOpt0(const G4double K);
    60   G4double GetOpt1(const G4double K);
    61   G4double GetOpt2(const G4double K);
    62   G4double GetOpt3(const G4double K);
     65  G4double GetOpt0(G4double K);
     66  G4double GetOpt1(G4double K);
     67  G4double GetOpt2(G4double K);
     68  G4double GetOpt3(G4double K);
    6369
     70  virtual G4double CalcAlphaParam(const G4Fragment & fragment)  ;
     71 
     72  virtual G4double CalcBetaParam(const G4Fragment & fragment) ;
     73 
     74  G4double CCoeficient(G4int aZ) ;
    6475 
    65  virtual G4double CalcAlphaParam(const G4Fragment & fragment)  ;
    66  
    67  virtual G4double CalcBetaParam(const G4Fragment & fragment) ;
    68  
    69   G4double CCoeficient(const G4double aZ) ;
    70  
    71 //data members
     76  //data members
    7277
    73       G4ProtonCoulombBarrier theCoulombBarrier;
     78  G4ProtonCoulombBarrier theCoulombBarrier;
    7479
    75       G4double ResidualA;
    76       G4double ResidualZ;
    77       G4double theA;
    78       G4double theZ;
    79       G4double ResidualAthrd;
    80       G4double FragmentA;
    81       G4double FragmentAthrd;
    82       G4double U;
    83 
    84 
    85 
     80  G4int ResidualA;
     81  G4int ResidualZ;
     82  G4int theA;
     83  G4int theZ;
     84  G4double ResidualAthrd;
     85  G4int FragmentA;
     86  G4double FragmentAthrd;
     87  G4double U;
    8688
    8789};

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4TritonEvaporationChannel.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4TritonEvaporationChannel.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4TritonEvaporationChannel.hh,v 1.9 2010/11/17 12:14:59 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Nov. 1999)
    3232//
     33// 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup
    3334
    3435#ifndef G4TritonEvaporationChannel_h
    … …  
    4344public:
    4445  // only available constructor
    45   G4TritonEvaporationChannel() : G4EvaporationChannel(3,1,"triton",
    46                                                       &theEvaporationProbability,&theCoulombBarrier) {};
     46  G4TritonEvaporationChannel();
    4747
    4848  // destructor
    49   ~G4TritonEvaporationChannel() {};
     49  virtual ~G4TritonEvaporationChannel();
    5050
    5151private:
    … …  
    5454  G4TritonEvaporationChannel(const G4TritonEvaporationChannel & right);
    5555
    56 public:
    5756  G4bool operator==(const G4TritonEvaporationChannel & right) const;
    5857  G4bool operator!=(const G4TritonEvaporationChannel & right) const;
    59 
    60 private:
    6158
    6259  G4TritonCoulombBarrier theCoulombBarrier;

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4TritonEvaporationProbability.hh

    r962 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4TritonEvaporationProbability.hh,v 1.14 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
    3033//
     34// Modified:
     35// 17-11-2010 V.Ivanchenko integer Z and A
     36
     37
    3138#ifndef G4TritonEvaporationProbability_h
    3239#define G4TritonEvaporationProbability_h 1
    … …  
    4249  G4TritonEvaporationProbability();
    4350
    44   ~G4TritonEvaporationProbability() {}
     51  virtual ~G4TritonEvaporationProbability();
     52
    4553private: 
    4654  // Copy constructor
    … …  
    5462private:
    5563
    56   virtual G4double CrossSection(const  G4Fragment & fragment, const G4double K);
     64  virtual G4double CrossSection(const G4Fragment & fragment, G4double K);
    5765
    58   G4double GetOpt0(const G4double K);
    59   G4double GetOpt12(const G4double K);
    60   G4double GetOpt34(const G4double K);
     66  G4double GetOpt0(G4double K);
     67  G4double GetOpt12(G4double K);
     68  G4double GetOpt34(G4double K);
    6169
    6270 
    63  virtual G4double CalcAlphaParam(const G4Fragment & fragment)  ;
     71  virtual G4double CalcAlphaParam(const G4Fragment & fragment)  ;
    6472 
    65  virtual G4double CalcBetaParam(const G4Fragment & fragment)  ;
     73  virtual G4double CalcBetaParam(const G4Fragment & fragment)  ;
    6674 
    67   G4double CCoeficient(const G4double aZ) ;
     75  G4double CCoeficient(G4int aZ) ;
    6876 
    69 //data members
     77  //data members
     78     
     79  G4TritonCoulombBarrier theCoulombBarrier;
    7080
    71      
    72       G4TritonCoulombBarrier theCoulombBarrier;
    73 
    74       G4double ResidualA;
    75       G4double ResidualZ;
    76       G4double theA;
    77       G4double theZ;
    78       G4double ResidualAthrd;
    79       G4double FragmentA;
    80       G4double FragmentAthrd;
    81 
     81  G4int ResidualA;
     82  G4int ResidualZ;
     83  G4int theA;
     84  G4int theZ;
     85  G4double ResidualAthrd;
     86  G4int FragmentA;
     87  G4double FragmentAthrd;
    8288
    8389};
    8490
    85 
    8691#endif
    8792
    88 

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4UnstableFragmentBreakUp.hh

    r1340 r1347  
    2525//
    2626// $Id: G4UnstableFragmentBreakUp.hh,v 1.2 2010/05/11 11:26:15 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2828//
    2929// -------------------------------------------------------------------

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4VEvaporation.hh

    r1340 r1347  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4VEvaporation.hh,v 1.6 2010/04/27 14:00:40 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     26// $Id: G4VEvaporation.hh,v 1.8 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2928//
    3029// Hadronic Process: Nuclear De-excitations
    … …  
    6665  // for superimposed Coulomb Barrier for inverse cross sections       
    6766  inline void UseSICB(G4bool use) { useSICB = use; }   
     67
    6868protected:
    69    G4int OPTxs;
    70    G4bool useSICB;
     69  G4int OPTxs;
     70  G4bool useSICB;
    7171
    7272};

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationChannel.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4AlphaEvaporationChannel.cc,v 1.4 2006/06/29 20:10:17 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4AlphaEvaporationChannel.cc,v 1.5 2010/11/17 12:14:59 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Nov. 1999)
    3232//
     33// 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup
    3334
    3435#include "G4AlphaEvaporationChannel.hh"
    3536
     37G4AlphaEvaporationChannel::G4AlphaEvaporationChannel()
     38: G4EvaporationChannel(4,2,"alpha",&theEvaporationProbability,&theCoulombBarrier)
     39{}
    3640
    37 const G4AlphaEvaporationChannel & G4AlphaEvaporationChannel::operator=(const G4AlphaEvaporationChannel & )
    38 {
    39   throw G4HadronicException(__FILE__, __LINE__, "G4AlphaEvaporationChannel::operator= meant to not be accessable");
    40   return *this;
    41 }
     41G4AlphaEvaporationChannel::~G4AlphaEvaporationChannel()
     42{}
    4243
    43 G4AlphaEvaporationChannel::G4AlphaEvaporationChannel(const G4AlphaEvaporationChannel & ) : G4EvaporationChannel()
    44 {
    45   throw G4HadronicException(__FILE__, __LINE__, "G4AlphaEvaporationChannel::CopyConstructor meant to not be accessable");
    46 }
    47 
    48 G4bool G4AlphaEvaporationChannel::operator==(const G4AlphaEvaporationChannel &right ) const
    49 {
    50   return (this == (G4AlphaEvaporationChannel *) &right);
    51   //  return false;
    52 }
    53 
    54 G4bool G4AlphaEvaporationChannel::operator!=(const G4AlphaEvaporationChannel &right ) const
    55 {
    56   return (this != (G4AlphaEvaporationChannel *) &right);
    57   //  return true;
    58 }
    59 

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc

    r1315 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4AlphaEvaporationProbability.cc,v 1.18 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
    3033//
    31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse
    32 // cross section option
     34// Modified:
     35// 03-09-2008 J.M. Quesada for external choice of inverse cross section option
     36// 17-11-2010 V.Ivanchenko integer Z and A
    3337
    3438#include "G4AlphaEvaporationProbability.hh"
    … …  
    3640G4AlphaEvaporationProbability::G4AlphaEvaporationProbability() :
    3741    G4EvaporationProbability(4,2,1,&theCoulombBarrier) // A,Z,Gamma,&theCoumlombBarrier
    38 {
     42{}
     43
     44G4AlphaEvaporationProbability::~G4AlphaEvaporationProbability()
     45{}
     46
     47G4double G4AlphaEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
     48  { return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());}
    3949       
    40 }
    41 
    42 
    43 G4AlphaEvaporationProbability::G4AlphaEvaporationProbability(const G4AlphaEvaporationProbability &): G4EvaporationProbability()
    44 {
    45     throw G4HadronicException(__FILE__, __LINE__, "G4AlphaEvaporationProbability::copy_constructor meant to not be accessable");
    46 }
    47 
    48 
    49 
    50 
    51 const G4AlphaEvaporationProbability & G4AlphaEvaporationProbability::
    52 operator=(const G4AlphaEvaporationProbability &)
    53 {
    54     throw G4HadronicException(__FILE__, __LINE__, "G4AlphaEvaporationProbability::operator= meant to not be accessable");
    55     return *this;
    56 }
    57 
    58 
    59 G4bool G4AlphaEvaporationProbability::operator==(const G4AlphaEvaporationProbability &) const
    60 {
    61     return false;
    62 }
    63 
    64 G4bool G4AlphaEvaporationProbability::operator!=(const G4AlphaEvaporationProbability &) const
    65 {
    66     return true;
    67 }
    68 
    69 
    70  G4double G4AlphaEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
    71   { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
     50G4double G4AlphaEvaporationProbability::CalcBetaParam(const G4Fragment &)
     51  { return 0.0; }
     52
     53G4double G4AlphaEvaporationProbability::CCoeficient(G4int aZ)
     54{
     55  // Data comes from
     56  // Dostrovsky, Fraenkel and Friedlander
     57  // Physical Review, vol 116, num. 3 1959
     58  //
     59  // const G4int size = 5;
     60  // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
     61  //    G4double Calpha[5] = { 0.10, 0.10, 0.10, 0.08, 0.06};
     62  G4double C = 0.0;
    7263       
    73  G4double G4AlphaEvaporationProbability::CalcBetaParam(const G4Fragment &)
    74   { return 0.0; }
    75 
    76 G4double G4AlphaEvaporationProbability::CCoeficient(const G4double aZ)
    77 {
    78     // Data comes from
    79     // Dostrovsky, Fraenkel and Friedlander
    80     // Physical Review, vol 116, num. 3 1959
    81     //
    82     // const G4int size = 5;
    83     // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
    84     //  G4double Calpha[5] = { 0.10, 0.10, 0.10, 0.08, 0.06};
    85     G4double C = 0.0;
    86        
    87        
    88     if (aZ <= 30) {
    89         C = 0.10;
    90     } else if (aZ <= 50) {
    91         C = 0.1 + -((aZ-50.)/20.)*0.02;
    92     } else if (aZ < 70) {
    93         C = 0.08 + -((aZ-70.)/20.)*0.02;
    94     } else {
    95         C = 0.06;
    96     }
    97     return C;
     64  if (aZ <= 30)
     65    {
     66      C = 0.10;
     67    }
     68  else if (aZ <= 50)
     69    {
     70      C = 0.1 - (aZ-30)*0.001;
     71    }
     72  else if (aZ < 70)
     73    {
     74      C = 0.08 - (aZ-50)*0.001;
     75    }
     76  else
     77    {
     78      C = 0.06;
     79    }
     80  return C;
    9881}
    9982
    … …  
    10487//OPT=3,4 Kalbach's parameterization
    10588//
    106 G4double G4AlphaEvaporationProbability::CrossSection(const  G4Fragment & fragment,
    107                                                      const G4double K)
     89G4double
     90G4AlphaEvaporationProbability::CrossSection(const G4Fragment & fragment, G4double K)
    10891{
    10992  theA=GetA();
    11093  theZ=GetZ();
    111   ResidualA=fragment.GetA()-theA;
    112   ResidualZ=fragment.GetZ()-theZ;
    113  
    114   ResidualAthrd=std::pow(ResidualA,0.33333);
    115   FragmentA=fragment.GetA();
    116   FragmentAthrd=std::pow(FragmentA,0.33333);
    117  
     94  ResidualA=fragment.GetA_asInt()-theA;
     95  ResidualZ=fragment.GetZ_asInt()-theZ;
     96 
     97  ResidualAthrd=fG4pow->Z13(ResidualA);
     98  FragmentA=fragment.GetA_asInt();
     99  FragmentAthrd=fG4pow->Z13(FragmentA);
    118100 
    119101  if (OPTxs==0) {std::ostringstream errOs;
    … …  
    133115}
    134116
    135 
    136 //
    137 //********************* OPT=1,2 : Chatterjee's cross section ************************
     117//
     118//********************* OPT=1,2 : Chatterjee's cross section ********************
    138119//(fitting to cross section from Bechetti & Greenles OM potential)
    139120
    140 G4double G4AlphaEvaporationProbability::GetOpt12(const  G4double K)
    141 {
    142 // c     ** alpha from huizenga and igo
     121G4double G4AlphaEvaporationProbability::GetOpt12(G4double K)
     122{
    143123  G4double Kc=K;
    144  
    145   // JMQ xsec is set constat above limit of validity
    146   if (K>50) Kc=50;
    147  
     124
     125  // JMQ xsec is set constant above limit of validity
     126  if (K > 50*MeV) { Kc = 50*MeV; }
     127
    148128  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
    149  
     129
    150130  G4double     p0 = 10.95;
    151131  G4double     p1 = -85.2;
    … …  
    160140  G4double     delta=1.2;         
    161141
    162  
    163142  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
    164143  p = p0 + p1/Ec + p2/(Ec*Ec);
    165144  landa = landa0*ResidualA + landa1;
    166   mu = mu0*std::pow(ResidualA,mu1);
    167   nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
     145  G4double resmu1 = fG4pow->powZ(ResidualA,mu1);
     146  mu = mu0*resmu1;
     147  nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
    168148  q = landa - nu/(Ec*Ec) - 2*p*Ec;
    169149  r = mu + 2*nu/Ec + p*(Ec*Ec);
    170  
     150
    171151  ji=std::max(Kc,Ec);
    172152  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
    173153  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
    174                  
     154 
    175155  if (xs <0.0) {xs=0.0;}
    176  
     156             
    177157  return xs;
    178  
    179158}
    180159
    181160// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
    182 G4double G4AlphaEvaporationProbability::GetOpt34(const  G4double K)
     161G4double G4AlphaEvaporationProbability::GetOpt34(G4double K)
    183162// c     ** alpha from huizenga and igo
    184163{
    185  
    186164  G4double landa, mu, nu, p , signor(1.),sig;
    187165  G4double ec,ecsq,xnulam,etest(0.),a;
    188166  G4double b,ecut,cut,ecut2,geom,elab;
    189167
    190 
    191168  G4double     flow = 1.e-18;
    192169  G4double     spill= 1.e+18;
    193170
    194   G4double       p0 = 10.95;
     171  G4double     p0 = 10.95;
    195172  G4double     p1 = -85.2;
    196173  G4double     p2 = 1146.;
    … …  
    204181 
    205182  G4double      ra=1.20;
    206  
     183       
    207184  //JMQ 13/02/09 increase of reduced radius to lower the barrier
    208185  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
    … …  
    211188  p = p0 + p1/ec + p2/ecsq;
    212189  landa = landa0*ResidualA + landa1;
    213   a = std::pow(ResidualA,mu1);
     190  a = fG4pow->powZ(ResidualA,mu1);
    214191  mu = mu0 * a;
    215192  nu = a* (nu0+nu1*ec+nu2*ecsq); 
    216193  xnulam = nu / landa;
    217   if (xnulam > spill) xnulam=0.;
    218   if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
    219  
     194  if (xnulam > spill) { xnulam=0.; }
     195  if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
     196
    220197  a = -2.*p*ec + landa - nu/ecsq;
    221198  b = p*ecsq + mu + 2.*nu/ec;
    222199  ecut = 0.;
    223200  cut = a*a - 4.*p*b;
    224   if (cut > 0.) ecut = std::sqrt(cut);
     201  if (cut > 0.) { ecut = std::sqrt(cut); }
    225202  ecut = (ecut-a) / (p+p);
    226203  ecut2 = ecut;
    227 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
    228 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
    229 //  if (cut < 0.) ecut2 = ecut - 2.;
    230   if (cut < 0.) ecut2 = ecut;
    231   elab = K * FragmentA / ResidualA;
     204  //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     205  // ecut<0 means that there is no cut with energy axis, i.e. xs is set
     206  // to 0 bellow minimum
     207  //  if (cut < 0.) ecut2 = ecut - 2.;
     208  if (cut < 0.) { ecut2 = ecut; }
     209  elab = K * FragmentA / G4double(ResidualA);
    232210  sig = 0.;
    233211 
    234212  if (elab <= ec) { //start for E<Ec
    235     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     213    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
    236214  }           //end for E<Ec
    237215  else {           //start for E>Ec
    238216    sig = (landa*elab+mu+nu/elab) * signor;
    239217    geom = 0.;
    240     if (xnulam < flow || elab < etest) return sig;
     218    if (xnulam < flow || elab < etest) { return sig; }
    241219    geom = std::sqrt(theA*K);
    242220    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
    … …  
    245223  }           //end for E>Ec
    246224  return sig;
    247  
    248 }
    249 
    250 //   ************************** end of cross sections *******************************
     225}
     226

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationChannel.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4DeuteronEvaporationChannel.cc,v 1.4 2006/06/29 20:10:21 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4DeuteronEvaporationChannel.cc,v 1.5 2010/11/17 12:14:59 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Nov. 1999)
    3232//
     33// 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup
    3334
    3435#include "G4DeuteronEvaporationChannel.hh"
    3536
     37G4DeuteronEvaporationChannel::G4DeuteronEvaporationChannel()
     38: G4EvaporationChannel(2,1,"deuteron",&theEvaporationProbability,&theCoulombBarrier)
     39{}
    3640
    37 const G4DeuteronEvaporationChannel & G4DeuteronEvaporationChannel::operator=(const G4DeuteronEvaporationChannel & )
    38 {
    39     throw G4HadronicException(__FILE__, __LINE__, "G4DeuteronEvaporationChannel::operator= meant to not be accessable");
    40     return *this;
    41 }
     41G4DeuteronEvaporationChannel::~G4DeuteronEvaporationChannel()
     42{}
    4243
    43 G4DeuteronEvaporationChannel::G4DeuteronEvaporationChannel(const G4DeuteronEvaporationChannel & ) : G4EvaporationChannel()
    44 {
    45     throw G4HadronicException(__FILE__, __LINE__, "G4DeuteronEvaporationChannel::CopyConstructor meant to not be accessable");
    46 }
    47 
    48 G4bool G4DeuteronEvaporationChannel::operator==(const G4DeuteronEvaporationChannel & right) const
    49 {
    50     return (this == (G4DeuteronEvaporationChannel *) &right);
    51     //  return false;
    52 }
    53 
    54 G4bool G4DeuteronEvaporationChannel::operator!=(const G4DeuteronEvaporationChannel & right) const
    55 {
    56     return (this != (G4DeuteronEvaporationChannel *) &right);
    57     //  return true;
    58 }

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationProbability.cc

    r1315 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4DeuteronEvaporationProbability.cc,v 1.20 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
    3033//
    31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse
    32 // cross section option
     34// Modified:
     35// 03-09-2008 J.M. Quesada for external choice of inverse cross section option
     36// 17-11-2010 V.Ivanchenko integer Z and A
    3337
    3438
    … …  
    3842G4DeuteronEvaporationProbability::G4DeuteronEvaporationProbability() :
    3943    G4EvaporationProbability(2,1,3,&theCoulombBarrier) // A,Z,Gamma (fixed JMQ)
    40 {       }
    41 G4DeuteronEvaporationProbability::G4DeuteronEvaporationProbability(const G4DeuteronEvaporationProbability &) : G4EvaporationProbability()
    42 {
    43     throw G4HadronicException(__FILE__, __LINE__, "G4DeuteronEvaporationProbability::copy_constructor meant to not be accessable");
    44 }
    45 
    46 
    47 const G4DeuteronEvaporationProbability & G4DeuteronEvaporationProbability::
    48 operator=(const G4DeuteronEvaporationProbability &)
    49 {
    50     throw G4HadronicException(__FILE__, __LINE__, "G4DeuteronEvaporationProbability::operator= meant to not be accessable");
    51     return *this;
    52 }
    53 
    54 
    55 G4bool G4DeuteronEvaporationProbability::operator==(const G4DeuteronEvaporationProbability &) const
    56 {
    57     return false;
    58 }
    59 
    60 
    61 G4bool G4DeuteronEvaporationProbability::operator!=(const G4DeuteronEvaporationProbability &) const
    62 {
    63     return true;
    64 }
    65 
     44{}
     45
     46G4DeuteronEvaporationProbability::~G4DeuteronEvaporationProbability()
     47{}
    6648
    6749G4double G4DeuteronEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
    6850{
    69     return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));
    70 }
    71 
     51  return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());
     52}
    7253
    7354G4double G4DeuteronEvaporationProbability::CalcBetaParam(const G4Fragment & )
    7455{
    75     return 0.0;
    76 }
    77 
    78 
    79 G4double G4DeuteronEvaporationProbability::CCoeficient(const G4double aZ)
    80 {
    81     // Data comes from
    82     // Dostrovsky, Fraenkel and Friedlander
    83     // Physical Review, vol 116, num. 3 1959
    84     //
    85     // const G4int size = 5;
    86     // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
    87     // G4double Cp[5] = { 0.50, 0.28, 0.20, 0.15, 0.10};
    88     // C for deuteron is equal to C for protons divided by 2
    89     G4double C = 0.0;
     56  return 0.0;
     57}
     58
     59G4double G4DeuteronEvaporationProbability::CCoeficient(G4int aZ)
     60{
     61  // Data comes from
     62  // Dostrovsky, Fraenkel and Friedlander
     63  // Physical Review, vol 116, num. 3 1959
     64  //
     65  // const G4int size = 5;
     66  // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
     67  // G4double Cp[5] = { 0.50, 0.28, 0.20, 0.15, 0.10};
     68  // C for deuteron is equal to C for protons divided by 2
     69  G4double C = 0.0;
    9070       
    91     if (aZ >= 70) {
    92         C = 0.10;
    93     } else {
    94         C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
    95     }
     71  if (aZ >= 70) {
     72    C = 0.10;
     73  } else {
     74    C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
     75  }
    9676       
    97     return C/2.0;
     77  return C/2.0;
    9878       
    9979}
    … …  
    10686//OPT=3,4 Kalbach's parameterization
    10787//
    108 G4double G4DeuteronEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
    109 {
    110        theA=GetA();
    111        theZ=GetZ();
    112        ResidualA=fragment.GetA()-theA;
    113        ResidualZ=fragment.GetZ()-theZ;
     88G4double
     89G4DeuteronEvaporationProbability::CrossSection(const  G4Fragment & fragment, G4double K)
     90{
     91  theA=GetA();
     92  theZ=GetZ();
     93  ResidualA=fragment.GetA_asInt()-theA;
     94  ResidualZ=fragment.GetZ_asInt()-theZ;
     95 
     96  ResidualAthrd=fG4pow->Z13(ResidualA);
     97  FragmentA=fragment.GetA_asInt();
     98  FragmentAthrd=fG4pow->Z13(FragmentA);
     99
     100  if (OPTxs==0) {std::ostringstream errOs;
     101  errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (deuterons)!!" 
     102        <<G4endl;
     103  throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     104  return 0.;}
     105  if( OPTxs==1 || OPTxs==2) return G4DeuteronEvaporationProbability::GetOpt12( K);
     106  else if (OPTxs==3 || OPTxs==4)  return G4DeuteronEvaporationProbability::GetOpt34( K);
     107  else{
     108    std::ostringstream errOs;
     109    errOs << "BAD Deuteron CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
     110    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     111    return 0.;
     112  }
     113}
     114
     115//
     116//********************* OPT=1,2 : Chatterjee's cross section ********************
     117//(fitting to cross section from Bechetti & Greenles OM potential)
     118
     119G4double G4DeuteronEvaporationProbability::GetOpt12(G4double K)
     120{
     121  G4double Kc = K;
     122
     123  // JMQ xsec is set constat above limit of validity
     124  if (K > 50*MeV) { Kc = 50*MeV; }
     125
     126  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
    114127 
    115        ResidualAthrd=std::pow(ResidualA,0.33333);
    116        FragmentA=fragment.GetA();
    117        FragmentAthrd=std::pow(FragmentA,0.33333);
    118 
    119        if (OPTxs==0) {std::ostringstream errOs;
    120          errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (deuterons)!!"  <<G4endl;
    121          throw G4HadronicException(__FILE__, __LINE__, errOs.str());
    122          return 0.;}
    123        if( OPTxs==1 || OPTxs==2) return G4DeuteronEvaporationProbability::GetOpt12( K);
    124        else if (OPTxs==3 || OPTxs==4)  return G4DeuteronEvaporationProbability::GetOpt34( K);
    125        else{
    126          std::ostringstream errOs;
    127          errOs << "BAD Deuteron CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
    128          throw G4HadronicException(__FILE__, __LINE__, errOs.str());
    129          return 0.;
    130        }
    131 }
    132 
    133 
    134 //********************* OPT=1,2 : Chatterjee's cross section ************************
    135 //(fitting to cross section from Bechetti & Greenles OM potential)
    136 
    137 G4double G4DeuteronEvaporationProbability::GetOpt12(const  G4double K)
    138 {
    139 
    140   G4double Kc=K;
    141 
    142 // JMQ xsec is set constat above limit of validity
    143  if (K>50) Kc=50;
    144 
    145  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
    146 
    147  
    148  G4double    p0 = -38.21;
    149  G4double    p1 = 922.6;
    150  G4double    p2 = -2804.;
    151  G4double    landa0 = -0.0323;
    152  G4double    landa1 = -5.48;
    153  G4double    mu0 = 336.1;
    154  G4double    mu1 = 0.48;
    155  G4double    nu0 = 524.3;
    156  G4double    nu1 = -371.8;
    157  G4double    nu2 = -5.924; 
    158  G4double    delta=1.2;           
    159  
    160 
    161  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
    162  p = p0 + p1/Ec + p2/(Ec*Ec);
    163  landa = landa0*ResidualA + landa1;
    164  mu = mu0*std::pow(ResidualA,mu1);
    165  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
    166  q = landa - nu/(Ec*Ec) - 2*p*Ec;
    167  r = mu + 2*nu/Ec + p*(Ec*Ec);
    168  
    169  ji=std::max(Kc,Ec);
    170  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
    171  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
     128  G4double    p0 = -38.21;
     129  G4double    p1 = 922.6;
     130  G4double    p2 = -2804.;
     131  G4double    landa0 = -0.0323;
     132  G4double    landa1 = -5.48;
     133  G4double    mu0 = 336.1;
     134  G4double    mu1 = 0.48;
     135  G4double    nu0 = 524.3;
     136  G4double    nu1 = -371.8;
     137  G4double    nu2 = -5.924; 
     138  G4double    delta=1.2;           
     139
     140  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
     141  p = p0 + p1/Ec + p2/(Ec*Ec);
     142  landa = landa0*ResidualA + landa1;
     143  G4double resmu1 = fG4pow->powZ(ResidualA,mu1);
     144  mu = mu0*resmu1;
     145  nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
     146  q = landa - nu/(Ec*Ec) - 2*p*Ec;
     147  r = mu + 2*nu/Ec + p*(Ec*Ec);
     148
     149  ji=std::max(Kc,Ec);
     150  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
     151  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
    172152                 
    173  if (xs <0.0) {xs=0.0;}
    174  
    175  return xs;
    176 
    177 }
    178 
     153  if (xs <0.0) {xs=0.0;}
     154             
     155  return xs;
     156}
    179157
    180158// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
    181 G4double G4DeuteronEvaporationProbability::GetOpt34(const  G4double K)
     159G4double G4DeuteronEvaporationProbability::GetOpt34(G4double K)
    182160//     ** d from o.m. of perey and perey
    183161{
    184162
    185   G4double landa, mu, nu, p , signor(1.),sig;
     163  G4double landa, mu, nu, p ,signor(1.),sig;
    186164  G4double ec,ecsq,xnulam,etest(0.),a;
    187165  G4double b,ecut,cut,ecut2,geom,elab;
    188166
    189 
    190167  G4double     flow = 1.e-18;
    191168  G4double     spill= 1.e+18;
    192 
    193169
    194170  G4double     p0 = 0.798;
    … …  
    202178  G4double     nu1 = -474.5;
    203179  G4double     nu2 = -3.592;     
    204   
    205   G4double      ra=0.80;
     180 
     181  G4double     ra=0.80;
    206182       
    207183  //JMQ 13/02/09 increase of reduced radius to lower the barrier
    208   //  ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
    209    ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
    210 
     184  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     185  ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
    211186  ecsq = ec * ec;
    212187  p = p0 + p1/ec + p2/ecsq;
    213188  landa = landa0*ResidualA + landa1;
    214   a = std::pow(ResidualA,mu1);
     189  a = fG4pow->powZ(ResidualA,mu1);
    215190  mu = mu0 * a;
    216191  nu = a* (nu0+nu1*ec+nu2*ecsq); 
    217192  xnulam = nu / landa;
    218   if (xnulam > spill) xnulam=0.;
    219   if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
     193  if (xnulam > spill) { xnulam=0.; }
     194  if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
    220195
    221196  a = -2.*p*ec + landa - nu/ecsq;
    … …  
    223198  ecut = 0.;
    224199  cut = a*a - 4.*p*b;
    225   if (cut > 0.) ecut = std::sqrt(cut);
     200  if (cut > 0.) { ecut = std::sqrt(cut); }
    226201  ecut = (ecut-a) / (p+p);
    227202  ecut2 = ecut;
    228 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
    229 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
    230 //  if (cut < 0.) ecut2 = ecut - 2.;
    231   if (cut < 0.) ecut2 = ecut;
    232   elab = K * FragmentA / ResidualA;
     203  //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     204  //ecut<0 means that there is no cut with energy axis, i.e. xs is set
     205  //to 0 bellow minimum
     206  //  if (cut < 0.) ecut2 = ecut - 2.;
     207  if (cut < 0.) { ecut2 = ecut; }
     208  elab = K * FragmentA / G4double(ResidualA);
    233209  sig = 0.;
    234210
    235211  if (elab <= ec) { //start for E<Ec
    236     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     212    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
    237213  }           //end for E<Ec
    238214  else {           //start for E>Ec
    239215    sig = (landa*elab+mu+nu/elab) * signor;
    240216    geom = 0.;
    241     if (xnulam < flow || elab < etest) return sig;
     217    if (xnulam < flow || elab < etest) { return sig; }
    242218    geom = std::sqrt(theA*K);
    243219    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
    … …  
    246222  }           //end for E>Ec
    247223  return sig;
    248  
    249 }
    250 
    251 //   ************************** end of cross sections *******************************
     224}
     225

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4Evaporation.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4Evaporation.cc,v 1.25 2010/06/09 11:56:47 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4Evaporation.cc,v 1.26 2010/11/23 18:10:10 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    142142    G4int Z = theResidualNucleus->GetZ_asInt();
    143143    G4double abun = nist->GetIsotopeAbundance(Z, A);
    144     /*
    145     G4cout << "### G4Evaporation::BreakItUp step " << ia << " Z= " << Z
    146            << " A= " << A << " Eex(MeV)= "
    147            << theResidualNucleus->GetExcitationEnergy()
    148            << " aban= " << abun << G4endl;
    149     */
     144   
     145    // G4cout << "### G4Evaporation::BreakItUp step " << ia << " Z= " << Z
     146    //     << " A= " << A << " Eex(MeV)= "
     147    //     << theResidualNucleus->GetExcitationEnergy()
     148    //     << " aban= " << abun << G4endl;
     149   
    150150    if(theResidualNucleus->GetExcitationEnergy() <= minExcitation &&
    151151       (abun > 0.0))
    … …  
    184184    // do evaporation chain and reset total probability
    185185    if(0.0 < totprob && probabilities[0] == totprob) {
     186      //G4cout << "Start gamma evaporation" << G4endl;
    186187      theTempResult = (*theChannels)[0]->BreakUpFragment(theResidualNucleus);
    187188      if(theTempResult) {
    … …  
    228229    // single photon evaporation, primary pointer is kept
    229230    if(0 == i) {
     231      //G4cout << "Single gamma" << G4endl;
    230232      G4Fragment* gamma = (*theChannels)[0]->EmittedFragment(theResidualNucleus);
    231233      if(gamma) { theResult->push_back(gamma); }
    … …  
    233235      // fission, return results to the main loop if fission is succesful
    234236    } else if(1 == i) {
     237      //G4cout << "Fission" << G4endl;
    235238      theTempResult = (*theChannels)[1]->BreakUp(*theResidualNucleus);
    236239      if(theTempResult) {
    … …  
    248251      // other channels
    249252    } else {
     253      //G4cout << "Channel # " << i << G4endl;
    250254      theTempResult = (*theChannels)[i]->BreakUp(*theResidualNucleus);
    251255      if(theTempResult) {

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationChannel.cc

    r962 r1347  
    2424// ********************************************************************
    2525//
     26// $Id: G4EvaporationChannel.cc,v 1.19 2010/11/24 15:30:49 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2628//
    2729//J.M. Quesada (August2008). Based on:
    … …  
    3032// by V. Lara (Oct 1998)
    3133//
    32 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse
    33 // cross section option
    34 // JMQ (06 September 2008) Also external choices have been added for
    35 // superimposed Coulomb barrier (if useSICB is set true, by default is false)
    36 
     34// Modified:
     35// 03-09-2008 J.M. Quesada for external choice of inverse cross section option
     36// 06-09-2008 J.M. Quesada Also external choices have been added for superimposed
     37//                 Coulomb barrier (if useSICB is set true, by default is false)
     38// 17-11-2010 V.Ivanchenko in constructor replace G4VEmissionProbability by
     39//            G4EvaporationProbability and do not new and delete probability
     40//            object at each call; use G4Pow
    3741
    3842#include "G4EvaporationChannel.hh"
    3943#include "G4PairingCorrection.hh"
    40 
    41 
    42 
    43 G4EvaporationChannel::G4EvaporationChannel(const G4int anA, const G4int aZ, const G4String & aName,
    44                                            G4VEmissionProbability * aEmissionStrategy,
    45                                            G4VCoulombBarrier * aCoulombBarrier):
     44#include "G4NucleiProperties.hh"
     45#include "G4Pow.hh"
     46#include "G4EvaporationLevelDensityParameter.hh"
     47#include "Randomize.hh"
     48#include "G4Alpha.hh"
     49
     50G4EvaporationChannel::G4EvaporationChannel(G4int anA, G4int aZ,
     51                                           const G4String & aName,
     52                                           G4EvaporationProbability* aEmissionStrategy,
     53                                           G4VCoulombBarrier* aCoulombBarrier):
    4654    G4VEvaporationChannel(aName),
    4755    theA(anA),
    … …  
    5260    MaximalKineticEnergy(-1000.0)
    5361{
    54     theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
    55 //    MyOwnLevelDensity = true; 
     62  ResidualA = 0;
     63  ResidualZ = 0;
     64  ResidualMass = CoulombBarrier=0.0;
     65  EvaporatedMass = G4NucleiProperties::GetNuclearMass(theA, theZ);
     66  theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
     67}
     68
     69G4EvaporationChannel::G4EvaporationChannel():
     70    G4VEvaporationChannel(""),
     71    theA(0),
     72    theZ(0),
     73    theEvaporationProbabilityPtr(0),
     74    theCoulombBarrierPtr(0),
     75    EmissionProbability(0.0),
     76    MaximalKineticEnergy(-1000.0)
     77{
     78  ResidualA = 0;
     79  ResidualZ = 0;
     80  EvaporatedMass = ResidualMass = CoulombBarrier = 0.0;
     81  theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
    5682}
    5783
    … …  
    6086  delete theLevelDensityPtr;
    6187}
    62 
    63 G4EvaporationChannel::G4EvaporationChannel(const G4EvaporationChannel & ) : G4VEvaporationChannel()
    64 {
    65     throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::copy_costructor meant to not be accessable");
    66 }
    67 
    68 const G4EvaporationChannel & G4EvaporationChannel::operator=(const G4EvaporationChannel & )
    69 {
    70     throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::operator= meant to not be accessable");
    71     return *this;
    72 }
    73 
    74 G4bool G4EvaporationChannel::operator==(const G4EvaporationChannel & right) const
    75 {
    76     return (this == (G4EvaporationChannel *) &right);
    77     //  return false;
    78 }
    79 
    80 G4bool G4EvaporationChannel::operator!=(const G4EvaporationChannel & right) const
    81 {
    82     return (this != (G4EvaporationChannel *) &right);
    83     //  return true;
    84 }
    85 
    86 //void G4EvaporationChannel::SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity)
    87 //  {
    88 //    if (MyOwnLevelDensity) delete theLevelDensityPtr;
    89 //    theLevelDensityPtr = aLevelDensity;
    90 //    MyOwnLevelDensity = false;
    91 //  }
    9288
    9389void G4EvaporationChannel::Initialize(const G4Fragment & fragment)
    … …  
    9894  theEvaporationProbabilityPtr->UseSICB(useSICB);
    9995 
    100  
    101   G4int FragmentA = static_cast<G4int>(fragment.GetA());
    102   G4int FragmentZ = static_cast<G4int>(fragment.GetZ());
     96  G4int FragmentA = fragment.GetA_asInt();
     97  G4int FragmentZ = fragment.GetZ_asInt();
    10398  ResidualA = FragmentA - theA;
    10499  ResidualZ = FragmentZ - theZ;
     100  //G4cout << "G4EvaporationChannel::Initialize Z= " << theZ << " A= " << theA
     101  //     << " FragZ= " << FragmentZ << " FragA= " << FragmentA << G4endl;
    105102 
    106103  //Effective excitation energy
    … …  
    108105    G4PairingCorrection::GetInstance()->GetPairingCorrection(FragmentA,FragmentZ);
    109106 
    110  
    111107  // Only channels which are physically allowed are taken into account
    112108  if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ ||
    113109      (ResidualA == ResidualZ && ResidualA > 1) || ExEnergy <= 0.0) {
    114     CoulombBarrier=0.0;
     110    CoulombBarrier = ResidualMass = 0.0;
    115111    MaximalKineticEnergy = -1000.0*MeV;
    116112    EmissionProbability = 0.0;
    117113  } else {
     114    ResidualMass = G4NucleiProperties::GetNuclearMass(ResidualA, ResidualZ);
     115    G4double FragmentMass = fragment.GetGroundStateMass();
    118116    CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy);
    119117    // Maximal Kinetic Energy
    120     MaximalKineticEnergy = CalcMaximalKineticEnergy
    121       (G4ParticleTable::GetParticleTable()->
    122        GetIonTable()->GetNucleusMass(FragmentZ,FragmentA)+ExEnergy);
     118    MaximalKineticEnergy = CalcMaximalKineticEnergy(FragmentMass + ExEnergy);
     119    //MaximalKineticEnergy = ExEnergy + fragment.GetGroundStateMass()
     120    //  - EvaporatedMass - ResidualMass;
    123121   
    124122    // Emission probability
    … …  
    131129      limit= CoulombBarrier;
    132130    else limit=0.;
     131    limit = std::max(limit, FragmentMass - ResidualMass - EvaporatedMass);
    133132 
    134133    // The threshold for charged particle emission must be  set to 0 if Coulomb
    135134    //cutoff  is included in the cross sections
    136     //if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0; 
    137135    if (MaximalKineticEnergy <= limit) EmissionProbability = 0.0;
    138136    else {
    … …  
    142140    }
    143141  }
     142  //G4cout << "G4EvaporationChannel:: probability= " << EmissionProbability << G4endl;
    144143 
    145144  return;
    … …  
    148147G4FragmentVector * G4EvaporationChannel::BreakUp(const G4Fragment & theNucleus)
    149148{
    150   G4double EvaporatedKineticEnergy=GetKineticEnergy(theNucleus);
    151  
    152   G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
    153     GetIonMass(theZ,theA);
    154   G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
    155  
     149  /*
     150  G4double Ecm = GetKineticEnergy(theNucleus) + ResidualMass + EvaporatedMass;
     151 
     152  G4double EvaporatedEnergy =
     153    ((Ecm-ResidualMass)*(Ecm+ResidualMass) + EvaporatedMass*EvaporatedMass)/(2*Ecm);
     154  */
     155  G4double EvaporatedEnergy = GetKineticEnergy(theNucleus) + EvaporatedMass;
     156
    156157  G4ThreeVector momentum(IsotropicVector
    157                          (std::sqrt(EvaporatedKineticEnergy*
    158                                     (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
    159  
    160   momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
     158                         (std::sqrt((EvaporatedEnergy - EvaporatedMass)*
     159                                    (EvaporatedEnergy + EvaporatedMass))));
    161160 
    162161  G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
    163   EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
     162  G4LorentzVector ResidualMomentum = theNucleus.GetMomentum();
     163  EvaporatedMomentum.boost(ResidualMomentum.boostVector());
    164164 
    165165  G4Fragment * EvaporatedFragment = new G4Fragment(theA,theZ,EvaporatedMomentum);
    166 #ifdef PRECOMPOUND_TEST
    167   EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
    168 #endif
    169   // ** And now the residual nucleus **
    170   G4double theExEnergy = theNucleus.GetExcitationEnergy();
    171   G4double theMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
    172     GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
    173                    static_cast<G4int>(theNucleus.GetA()));
    174   G4double ResidualEnergy = theMass +
    175     (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
    176  
    177   G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
    178   ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
    179  
    180   G4Fragment * ResidualFragment = new G4Fragment(ResidualA, ResidualZ, ResidualMomentum );
    181  
    182 #ifdef PRECOMPOUND_TEST
    183   ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
    184 #endif
     166  ResidualMomentum -= EvaporatedMomentum;
     167
     168  G4Fragment * ResidualFragment = new G4Fragment(ResidualA, ResidualZ, ResidualMomentum);
     169 
    185170  G4FragmentVector * theResult = new G4FragmentVector;
    186171 
    … …  
    211196/////////////////////////////////////////
    212197// Calculates the maximal kinetic energy that can be carried by fragment.
    213 G4double G4EvaporationChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
    214 {
    215   G4double ResidualMass = G4ParticleTable::GetParticleTable()->
    216     GetIonTable()->GetNucleusMass( ResidualZ, ResidualA );
    217   G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->
    218     GetIonTable()->GetNucleusMass( theZ, theA );
    219 
     198G4double G4EvaporationChannel::CalcMaximalKineticEnergy(G4double NucleusTotalE)
     199{
    220200  // This is the "true" assimptotic kinetic energy (from energy conservation)   
    221   G4double Tmax = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass -               
    222                    ResidualMass*ResidualMass)/(2.0*NucleusTotalE) - EvaporatedMass;
     201  G4double Tmax =
     202    ((NucleusTotalE-ResidualMass)*(NucleusTotalE+ResidualMass) + EvaporatedMass*EvaporatedMass)
     203    /(2.0*NucleusTotalE) - EvaporatedMass;
    223204 
    224205  //JMQ (13-09-08) bug fixed: in the original version the Tmax is calculated
    … …  
    245226   
    246227   
    247     if (MaximalKineticEnergy < 0.0)  
     228    if (MaximalKineticEnergy < 0.0) {
    248229      throw G4HadronicException(__FILE__, __LINE__,
    249230                                "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic at the Coulomb barrier is less than 0");
    250    
     231    }
    251232    G4double Rb = 4.0*theLevelDensityPtr->
    252233      LevelDensityParameter(ResidualA+theA,ResidualZ+theZ,MaximalKineticEnergy)*
    … …  
    263244      G4double Q2 = 1.0;
    264245      if (theZ == 0) { // for emitted neutron
    265         G4double Beta = (2.12/std::pow(G4double(ResidualA),2./3.) - 0.05)*MeV/
    266           (0.76 + 2.2/std::pow(G4double(ResidualA),1./3.));
     246        G4double Beta = (2.12/G4Pow::GetInstance()->Z23(ResidualA) - 0.05)*MeV/
     247          (0.76 + 2.2/G4Pow::GetInstance()->Z13(ResidualA));
    267248        Q1 = 1.0 + Beta/(MaximalKineticEnergy);
    268249        Q2 = Q1*std::sqrt(Q1);
    … …  
    277258  } else if (OPTxs==1 || OPTxs==2 || OPTxs==3 || OPTxs==4) {
    278259   
    279     G4double V;
    280     if(useSICB) V= CoulombBarrier;
    281     else V=0.;
    282     //If Coulomb barrier is just included  in the cross sections
    283     //  G4double V=0.;
     260    // Coulomb barrier is just included  in the cross sections
     261    G4double V = 0;
     262    if(useSICB) { V= CoulombBarrier; }
     263
     264    V = std::max(V, aFragment.GetGroundStateMass()-EvaporatedMass-ResidualMass);
    284265
    285266    G4double Tmax=MaximalKineticEnergy;
    286267    G4double T(0.0);
    287268    G4double NormalizedProbability(1.0);
    288    
     269
     270    // VI: This is very ineffective - create new objects at each call to the method   
     271    /*
    289272    // A pointer is created in order to access the distribution function.
    290     G4EvaporationProbability * G4EPtemp;
     273    G4EvaporationProbability * G4EPtemp = 0;
    291274   
    292275    if (theA==1 && theZ==0) G4EPtemp=new G4NeutronEvaporationProbability();
    … …  
    305288      G4EPtemp->SetOPTxs(OPTxs);
    306289      G4EPtemp->UseSICB(useSICB);
    307 
     290    */
     291
     292    // use local pointer and not create a new one
    308293    do
    309294      { 
    310295        T=V+G4UniformRand()*(Tmax-V);
    311         NormalizedProbability=G4EPtemp->ProbabilityDistributionFunction(aFragment,T)/
    312           (this->GetEmissionProbability());
     296        NormalizedProbability =
     297          theEvaporationProbabilityPtr->ProbabilityDistributionFunction(aFragment,T)/
     298          GetEmissionProbability();
    313299       
    314300      }
    315301    while (G4UniformRand() > NormalizedProbability);
    316    delete G4EPtemp;
    317    return T;
     302    //   delete G4EPtemp;
     303    return T;
    318304  } else{
    319305    std::ostringstream errOs;
    … …  
    323309}
    324310
    325 G4ThreeVector G4EvaporationChannel::IsotropicVector(const G4double Magnitude)
     311G4ThreeVector G4EvaporationChannel::IsotropicVector(G4double Magnitude)
    326312    // Samples a isotropic random vectorwith a magnitud given by Magnitude.
    327313    // By default Magnitude = 1.0
    328314{
    329     G4double CosTheta = 1.0 - 2.0*G4UniformRand();
    330     G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
    331     G4double Phi = twopi*G4UniformRand();
    332     G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
    333                         Magnitude*std::sin(Phi)*SinTheta,
    334                         Magnitude*CosTheta);
    335     return Vector;
    336             }
     315  G4double CosTheta = 1.0 - 2.0*G4UniformRand();
     316  G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
     317  G4double Phi = twopi*G4UniformRand();
     318  G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
     319                      Magnitude*std::sin(Phi)*SinTheta,
     320                      Magnitude*CosTheta);
     321  return Vector;
     322}
    337323
    338324

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationDefaultGEMFactory.cc

    r1340 r1347  
    2626//
    2727// $Id: G4EvaporationDefaultGEMFactory.cc,v 1.2 2010/04/27 11:43:16 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationFactory.cc

    r1340 r1347  
    2626//
    2727// $Id: G4EvaporationFactory.cc,v 1.5 2010/04/27 11:43:16 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.cc

    r1315 r1347  
    4444#include "G4IonTable.hh"
    4545
    46 
    47 G4EvaporationProbability::G4EvaporationProbability(const G4EvaporationProbability &) : G4VEmissionProbability()
     46G4EvaporationProbability::G4EvaporationProbability(G4int anA, G4int aZ,
     47                                                   G4double aGamma,
     48                                                   G4VCoulombBarrier * aCoulombBarrier)
     49  : theA(anA),
     50    theZ(aZ),
     51    Gamma(aGamma),
     52    theCoulombBarrierptr(aCoulombBarrier)
     53{}
     54
     55G4EvaporationProbability::G4EvaporationProbability()
     56  : theA(0),
     57    theZ(0),
     58    Gamma(0.0),
     59    theCoulombBarrierptr(0)
     60{}
     61
     62G4EvaporationProbability::~G4EvaporationProbability()
     63{}
     64 
     65G4double
     66G4EvaporationProbability::EmissionProbability(const G4Fragment & fragment, G4double anEnergy)
    4867{
    49     throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationProbability::copy_constructor meant to not be accessable");
    50 }
    51 
    52 
    53 
    54 
    55 const G4EvaporationProbability & G4EvaporationProbability::
    56 operator=(const G4EvaporationProbability &)
     68  G4double EmissionProbability = 0.0;
     69  G4double MaximalKineticEnergy = anEnergy;
     70
     71  if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
     72    EmissionProbability = CalculateProbability(fragment, MaximalKineticEnergy);
     73
     74  }
     75  return EmissionProbability;
     76}
     77
     78////////////////////////////////////
     79
     80// Computes the integrated probability of evaporation channel
     81G4double
     82G4EvaporationProbability::CalculateProbability(const G4Fragment & fragment,
     83                                               G4double MaximalKineticEnergy)
    5784{
    58     throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationProbability::operator= meant to not be accessable");
    59     return *this;
    60 }
    61 
    62 
    63 G4bool G4EvaporationProbability::operator==(const G4EvaporationProbability &) const
    64 {
    65     return false;
    66 }
    67 
    68 G4bool G4EvaporationProbability::operator!=(const G4EvaporationProbability &) const
    69 {
    70     return true;
    71 }
    72  
    73 G4double G4EvaporationProbability::EmissionProbability(const G4Fragment & fragment, const G4double anEnergy)
    74 {
    75     G4double EmissionProbability = 0.0;
    76     G4double MaximalKineticEnergy = anEnergy;
    77 
    78     if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
    79         EmissionProbability = CalculateProbability(fragment, MaximalKineticEnergy);
    80 
    81     }
    82     return EmissionProbability;
    83 }
    84 
    85 ////////////////////////////////////
    86 
    87 // Computes the integrated probability of evaporation channel
    88 G4double G4EvaporationProbability::CalculateProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy)
    89 {
    90     G4double ResidualA = fragment.GetA() - theA;
    91     G4double ResidualZ = fragment.GetZ() - theZ;
    92     G4double U = fragment.GetExcitationEnergy();
     85  G4int ResidualA = fragment.GetA_asInt() - theA;
     86  G4int ResidualZ = fragment.GetZ_asInt() - theZ;
     87  G4double U = fragment.GetExcitationEnergy();
    9388   
    94  if (OPTxs==0) {
    95 
    96        
    97     G4double NuclearMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(theZ,theA);
    98 
    99 
    100     G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()),
    101                                                                                static_cast<G4int>(fragment.GetZ()));
    102 
    103     G4double SystemEntropy = 2.0*std::sqrt(theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),
    104                                                                            static_cast<G4int>(fragment.GetZ()),U)*
    105                                       (U-delta0));
     89  if (OPTxs==0) {
     90
     91    G4double NuclearMass = fragment.ComputeGroundStateMass(theZ,theA);
     92
     93    G4double delta0 = fPairCorr->GetPairingCorrection(fragment.GetA_asInt(),
     94                                                      fragment.GetZ_asInt());
     95
     96    G4double SystemEntropy = 2.0*std::sqrt(
     97      theEvapLDPptr->LevelDensityParameter(fragment.GetA_asInt(),fragment.GetZ_asInt(),U)*
     98      (U-delta0));
    10699                                                                 
    107 
    108     G4double RN = 1.5*fermi;
     100    const G4double RN = 1.5*fermi;
    109101
    110102    G4double Alpha = CalcAlphaParam(fragment);
    … …  
    112104       
    113105    G4double Rmax = MaximalKineticEnergy;
    114     G4double a = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA),
    115                                                       static_cast<G4int>(ResidualZ),
    116                                                       Rmax);
    117     G4double GlobalFactor = static_cast<G4double>(Gamma) * (Alpha/(a*a)) *
    118         (NuclearMass*RN*RN*std::pow(ResidualA,2./3.))/
    119         (2.*pi* hbar_Planck*hbar_Planck);
     106    G4double a = theEvapLDPptr->LevelDensityParameter(ResidualA,ResidualZ,Rmax);
     107    G4double GlobalFactor = Gamma * Alpha/(a*a) *
     108        (NuclearMass*RN*RN*fG4pow->Z23(ResidualA))/
     109        (twopi* hbar_Planck*hbar_Planck);
    120110    G4double Term1 = (2.0*Beta*a-3.0)/2.0 + Rmax*a;
    121111    G4double Term2 = (2.0*Beta*a-3.0)*std::sqrt(Rmax*a) + 2.0*a*Rmax;
    122112       
    123113    G4double ExpTerm1 = 0.0;
    124     if (SystemEntropy <= 600.0) ExpTerm1 = std::exp(-SystemEntropy);
     114    if (SystemEntropy <= 600.0) { ExpTerm1 = std::exp(-SystemEntropy); }
    125115       
    126116    G4double ExpTerm2 = 2.*std::sqrt(a*Rmax) - SystemEntropy;
    127     if (ExpTerm2 > 700.0) ExpTerm2 = 700.0;
     117    if (ExpTerm2 > 700.0) { ExpTerm2 = 700.0; }
    128118    ExpTerm2 = std::exp(ExpTerm2);
    129119       
    … …  
    134124 } else if (OPTxs==1 || OPTxs==2 ||OPTxs==3 || OPTxs==4) {
    135125
    136    G4double limit;
    137    if (useSICB) limit=theCoulombBarrierptr->GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U);
    138    else limit=0.;
    139 
    140    if (MaximalKineticEnergy <= limit)  return 0.0;
    141 
    142 
    143    // if Coulomb barrier cutoff is superimposed for all cross sections the limit is the Coulomb Barrier
     126   G4double EvaporatedMass = fragment.ComputeGroundStateMass(theZ,theA);
     127   G4double ResidulalMass = fragment.ComputeGroundStateMass(ResidualZ,ResidualA);
     128   G4double limit = std::max(0.0,fragment.GetGroundStateMass()-EvaporatedMass-ResidulalMass);
     129   if (useSICB) {
     130     limit = std::max(limit,theCoulombBarrierptr->GetCoulombBarrier(ResidualA,ResidualZ,U));
     131   }
     132
     133   if (MaximalKineticEnergy <= limit) { return 0.0; }
     134
     135   // if Coulomb barrier cutoff is superimposed for all cross sections
     136   // then the limit is the Coulomb Barrier
    144137   G4double LowerLimit= limit;
    145138
    146    //  MaximalKineticEnergy: asimptotic value (already accounted for in G4EvaporationChannel)     
     139   //MaximalKineticEnergy: asimptotic value (already accounted for in G4EvaporationChannel)     
    147140
    148141   G4double UpperLimit = MaximalKineticEnergy;
    149142
    150 
    151143   G4double Width = IntegrateEmissionProbability(fragment,LowerLimit,UpperLimit);
    152144
    153145   return Width;
    154  } else{
     146 } else {
    155147   std::ostringstream errOs;
    156148   errOs << "Bad option for cross sections at evaporation"  <<G4endl;
    … …  
    163155
    164156G4double G4EvaporationProbability::
    165 IntegrateEmissionProbability(const G4Fragment & fragment, const G4double & Low, const G4double & Up )
     157IntegrateEmissionProbability(const G4Fragment & fragment,
     158                             const G4double & Low, const G4double & Up )
    166159{
    167 
    168160  static const G4int N = 10;
    169161  // 10-Points Gauss-Legendre abcisas and weights
    … …  
    212204//New method (OPT=1,2,3,4)
    213205
    214 G4double G4EvaporationProbability::ProbabilityDistributionFunction( const G4Fragment & fragment, const G4double K)
     206G4double
     207G4EvaporationProbability::ProbabilityDistributionFunction( const G4Fragment & fragment,
     208                                                           G4double K)
    215209{
    216 
     210  G4int ResidualA = fragment.GetA_asInt() - theA;
     211  G4int ResidualZ = fragment.GetZ_asInt() - theZ; 
     212  G4double U = fragment.GetExcitationEnergy();
     213  //G4cout << "### G4EvaporationProbability::ProbabilityDistributionFunction" << G4endl;
     214  //G4cout << "FragZ= " << fragment.GetZ_asInt() << " FragA= " << fragment.GetA_asInt()
     215  //     << " Z= " << theZ << "  A= " << theA << G4endl;
     216  //G4cout << "PC " << fPairCorr << "   DP " << theEvapLDPptr << G4endl;
     217
     218  // if(K <= theCoulombBarrierptr->GetCoulombBarrier(ResidualA,ResidualZ,U)) return 0.0;
     219
     220  G4double delta1 = fPairCorr->GetPairingCorrection(ResidualA,ResidualZ);
    217221 
    218 
    219 
    220   G4double ResidualA = fragment.GetA() - theA;
    221   G4double ResidualZ = fragment.GetZ() - theZ; 
    222   G4double U = fragment.GetExcitationEnergy();
    223 
    224   //        if(K <= theCoulombBarrierptr->GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U)) return 0.0;   
    225 
    226   G4double delta1 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ));
    227 
    228  
    229   G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()));
    230 
    231  
    232   G4double ParticleMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(theZ,theA);
    233 
    234   G4double theSeparationEnergy= G4NucleiProperties::GetMassExcess(static_cast<G4int>(theA),static_cast<G4int>(theZ)) +
    235     G4NucleiProperties::GetMassExcess(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ)) -
    236     G4NucleiProperties::GetMassExcess(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()));
    237 
    238   G4double a0 = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()),U - delta0);
    239 
    240   G4double a1 = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ),U - theSeparationEnergy - delta1);
    241  
    242  
    243   G4double E0=U-delta0;
    244 
    245   G4double E1=U-theSeparationEnergy-delta1-K;
    246 
    247   if (E1<0.) return 0.;
     222  G4double delta0 = fPairCorr->GetPairingCorrection(fragment.GetA_asInt(),
     223                                                    fragment.GetZ_asInt());
     224
     225 
     226  G4double ParticleMass = fragment.ComputeGroundStateMass(theZ,theA);
     227  G4double ResidualMass = fragment.ComputeGroundStateMass(ResidualZ,ResidualA);
     228
     229  G4double theSeparationEnergy = ParticleMass + ResidualMass
     230    - fragment.GetGroundStateMass();
     231
     232  G4double a0 = theEvapLDPptr->LevelDensityParameter(fragment.GetA_asInt(),
     233                                                     fragment.GetZ_asInt(),
     234                                                     U - delta0);
     235
     236  G4double a1 = theEvapLDPptr->LevelDensityParameter(ResidualA, ResidualZ,
     237                                                     U - theSeparationEnergy - delta1);
     238 
     239 
     240  G4double E0 = U - delta0;
     241
     242  G4double E1 = U - theSeparationEnergy - delta1 - K;
     243
     244  if (E1<0.) { return 0.; }
    248245
    249246  //JMQ 14/02/09 BUG fixed: hbarc should be in the denominator instead of hbar_Planck
    250247  //Without 1/hbar_Panck remains as a width
    251   //  G4double  Prob=Gamma*ParticleMass/((pi*hbar_Planck)*(pi*hbar_Planck)*
    252   //std::exp(2*std::sqrt(a0*E0)))*K*CrossSection(fragment,K)*std::exp(2*std::sqrt(a1*E1))*millibarn;
    253248
    254249  G4double Prob=Gamma*ParticleMass/((pi*hbarc)*(pi*hbarc)*std::exp(2*std::sqrt(a0*E0)))

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationChannel.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4He3EvaporationChannel.cc,v 1.4 2006/06/29 20:10:33 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4He3EvaporationChannel.cc,v 1.5 2010/11/17 12:14:59 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Nov. 1999)
    3232//
     33// 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup
    3334
    3435#include "G4He3EvaporationChannel.hh"
    3536
     37G4He3EvaporationChannel::G4He3EvaporationChannel()
     38: G4EvaporationChannel(3,2,"He3",&theEvaporationProbability,&theCoulombBarrier)
     39{}
    3640
    37 const G4He3EvaporationChannel & G4He3EvaporationChannel::operator=(const G4He3EvaporationChannel & )
    38 {
    39     throw G4HadronicException(__FILE__, __LINE__, "G4He3EvaporationChannel::operator= meant to not be accessable");
    40     return *this;
    41 }
     41G4He3EvaporationChannel::~G4He3EvaporationChannel()
     42{}
    4243
    43 G4He3EvaporationChannel::G4He3EvaporationChannel(const G4He3EvaporationChannel & ) : G4EvaporationChannel()
    44 {
    45     throw G4HadronicException(__FILE__, __LINE__, "G4He3EvaporationChannel::CopyConstructor meant to not be accessable");
    46 }
    47 
    48 G4bool G4He3EvaporationChannel::operator==(const G4He3EvaporationChannel & right) const
    49 {
    50     return (this == (G4He3EvaporationChannel *) &right);
    51     //  return false;
    52 }
    53 
    54 G4bool G4He3EvaporationChannel::operator!=(const G4He3EvaporationChannel & right) const
    55 {
    56     return (this != (G4He3EvaporationChannel *) &right);
    57     //  return true;
    58 }
    59 

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationProbability.cc

    r1315 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4He3EvaporationProbability.cc,v 1.18 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
    3033//
    31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse
    32 // cross section option
     34// Modified:
     35// 03-09-2008 J.M. Quesada for external choice of inverse cross section option
     36// 17-11-2010 V.Ivanchenko integer Z and A
    3337
    3438#include "G4He3EvaporationProbability.hh"
    … …  
    3640G4He3EvaporationProbability::G4He3EvaporationProbability() :
    3741   G4EvaporationProbability(3,2,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
    38 {
    39 
    40 }
    41 
    42 G4He3EvaporationProbability::G4He3EvaporationProbability(const G4He3EvaporationProbability &) : G4EvaporationProbability()
    43 {
    44     throw G4HadronicException(__FILE__, __LINE__, "G4He3EvaporationProbability::copy_constructor meant to not be accessable");
    45 }
    46 
    47 
    48 
    49 
    50 const G4He3EvaporationProbability & G4He3EvaporationProbability::
    51 operator=(const G4He3EvaporationProbability &)
    52 {
    53     throw G4HadronicException(__FILE__, __LINE__, "G4He3EvaporationProbability::operator= meant to not be accessable");
    54     return *this;
    55 }
    56 
    57 
    58 G4bool G4He3EvaporationProbability::operator==(const G4He3EvaporationProbability &) const
    59 {
    60     return false;
    61 }
    62 
    63 G4bool G4He3EvaporationProbability::operator!=(const G4He3EvaporationProbability &) const
    64 {
    65     return true;
    66 }
    67 
    68   G4double G4He3EvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
    69   { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
     42{}
     43
     44G4He3EvaporationProbability::~G4He3EvaporationProbability()
     45{}
     46
     47G4double G4He3EvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
     48  { return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());}
    7049       
    71   G4double G4He3EvaporationProbability::CalcBetaParam(const G4Fragment & ) 
     50G4double G4He3EvaporationProbability::CalcBetaParam(const G4Fragment & ) 
    7251  { return 0.0; }
    7352
    7453
    75 G4double G4He3EvaporationProbability::CCoeficient(const G4double aZ)
    76 {
    77     // Data comes from
    78     // Dostrovsky, Fraenkel and Friedlander
    79     // Physical Review, vol 116, num. 3 1959
    80     //
    81     // const G4int size = 5;
    82     // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
    83     //  G4double Calpha[5] = { 0.10, 0.10, 0.10, 0.08, 0.06};
    84     // C for He3 is equal to C for alpha times 4/3
    85     G4double C = 0.0;
     54G4double G4He3EvaporationProbability::CCoeficient(G4int aZ)
     55{
     56  // Data comes from
     57  // Dostrovsky, Fraenkel and Friedlander
     58  // Physical Review, vol 116, num. 3 1959
     59  //
     60  // const G4int size = 5;
     61  // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
     62  //    G4double Calpha[5] = { 0.10, 0.10, 0.10, 0.08, 0.06};
     63  // C for He3 is equal to C for alpha times 4/3
     64  G4double C = 0.0;
    8665       
    87        
    88     if (aZ <= 30) {
    89         C = 0.10;
    90     } else if (aZ <= 50) {
    91         C = 0.1 + -((aZ-50.)/20.)*0.02;
    92     } else if (aZ < 70) {
    93         C = 0.08 + -((aZ-70.)/20.)*0.02;
    94     } else {
    95         C = 0.06;
    96     }
    97     return C*(4.0/3.0);
     66  if (aZ <= 30)
     67    {
     68      C = 0.10;
     69    }
     70  else if (aZ <= 50)
     71    {
     72      C = 0.1 - (aZ - 30)*0.001;
     73    }
     74  else if (aZ < 70)
     75    {
     76      C = 0.08 - (aZ - 50)*0.001;
     77    }
     78  else
     79    {
     80      C = 0.06;
     81    }
     82  return C*(4.0/3.0);
    9883}
    9984
    … …  
    10489//OPT=3,4 Kalbach's parameterization
    10590//
    106 G4double G4He3EvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
    107 {
    108        theA=GetA();
    109        theZ=GetZ();
    110        ResidualA=fragment.GetA()-theA;
    111        ResidualZ=fragment.GetZ()-theZ;
    112  
    113        ResidualAthrd=std::pow(ResidualA,0.33333);
    114        FragmentA=fragment.GetA();
    115        FragmentAthrd=std::pow(FragmentA,0.33333);
    116 
    117 
    118        if (OPTxs==0) {std::ostringstream errOs;
    119          errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (He3's)!!" 
    120                <<G4endl;
    121          throw G4HadronicException(__FILE__, __LINE__, errOs.str());
    122          return 0.;}
    123        if( OPTxs==1 || OPTxs==2) return G4He3EvaporationProbability::GetOpt12( K);
    124        else if (OPTxs==3 || OPTxs==4)  return G4He3EvaporationProbability::GetOpt34( K);
    125        else{
    126          std::ostringstream errOs;
    127          errOs << "BAD He3's CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
    128          throw G4HadronicException(__FILE__, __LINE__, errOs.str());
    129          return 0.;
    130        }
    131 }
    132 //
    133 //********************* OPT=1,2 : Chatterjee's cross section ************************
     91G4double
     92G4He3EvaporationProbability::CrossSection(const  G4Fragment & fragment, G4double K)
     93{
     94
     95  theA=GetA();
     96  theZ=GetZ();
     97  ResidualA=fragment.GetA_asInt()-theA;
     98  ResidualZ=fragment.GetZ_asInt()-theZ;
     99 
     100  ResidualAthrd=fG4pow->Z13(ResidualA);
     101  FragmentA=fragment.GetA_asInt();
     102  FragmentAthrd=fG4pow->Z13(FragmentA);
     103
     104  if (OPTxs==0) {std::ostringstream errOs;
     105  errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (He3's)!!" 
     106        <<G4endl;
     107  throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     108  return 0.;}
     109  if( OPTxs==1 || OPTxs==2) return G4He3EvaporationProbability::GetOpt12( K);
     110  else if (OPTxs==3 || OPTxs==4)  return G4He3EvaporationProbability::GetOpt34( K);
     111  else{
     112    std::ostringstream errOs;
     113    errOs << "BAD He3's CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
     114    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
     115    return 0.;
     116  }
     117}
     118
     119//********************* OPT=1,2 : Chatterjee's cross section *****************
    134120//(fitting to cross section from Bechetti & Greenles OM potential)
    135121
    136122G4double G4He3EvaporationProbability::GetOpt12(const  G4double K)
    137123{
    138 
    139 G4double Kc=K;
    140 
    141 // JMQ xsec is set constat above limit of validity
    142  if (K>50) Kc=50;
    143 
    144  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
    145 
    146  G4double     p0 = -3.06;
    147  G4double     p1 = 278.5;
    148  G4double     p2 = -1389.;
    149  G4double     landa0 = -0.00535;
    150  G4double     landa1 = -11.16;
    151  G4double     mu0 = 555.5;
    152  G4double     mu1 = 0.40;
    153  G4double     nu0 = 687.4;
    154  G4double     nu1 = -476.3;
    155  G4double     nu2 = 0.509;   
    156  G4double     delta=1.2;             
    157 
    158       Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
    159       p = p0 + p1/Ec + p2/(Ec*Ec);
    160       landa = landa0*ResidualA + landa1;
    161       mu = mu0*std::pow(ResidualA,mu1);
    162       nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
    163       q = landa - nu/(Ec*Ec) - 2*p*Ec;
    164       r = mu + 2*nu/Ec + p*(Ec*Ec);
    165 
    166    ji=std::max(Kc,Ec);
    167    if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
    168    else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
    169    
    170    if (xs <0.0) {xs=0.0;}
     124  G4double Kc = K;
     125
     126  // JMQ xsec is set constat above limit of validity
     127  if (K > 50*MeV) { Kc = 50*MeV; }
     128
     129  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
     130
     131  G4double     p0 = -3.06;
     132  G4double     p1 = 278.5;
     133  G4double     p2 = -1389.;
     134  G4double     landa0 = -0.00535;
     135  G4double     landa1 = -11.16;
     136  G4double     mu0 = 555.5;
     137  G4double     mu1 = 0.40;
     138  G4double     nu0 = 687.4;
     139  G4double     nu1 = -476.3;
     140  G4double     nu2 = 0.509;   
     141  G4double     delta=1.2;             
     142
     143  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
     144  p = p0 + p1/Ec + p2/(Ec*Ec);
     145  landa = landa0*ResidualA + landa1;
     146
     147  G4double resmu1 = fG4pow->powZ(ResidualA,mu1);
     148  mu = mu0*resmu1;
     149  nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
     150  q = landa - nu/(Ec*Ec) - 2*p*Ec;
     151  r = mu + 2*nu/Ec + p*(Ec*Ec);
     152 
     153  ji=std::max(Kc,Ec);
     154  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
     155  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
     156 
     157  if (xs <0.0) {xs=0.0;}
    171158             
    172    return xs;
     159  return xs;
    173160
    174161}
    … …  
    178165//c     ** 3he from o.m. of gibson et al
    179166{
    180  
    181167  G4double landa, mu, nu, p , signor(1.),sig;
    182168  G4double ec,ecsq,xnulam,etest(0.),a;
    183169  G4double b,ecut,cut,ecut2,geom,elab;
    184  
    185  
     170
    186171  G4double     flow = 1.e-18;
    187172  G4double     spill= 1.e+18;
    188 
    189173
    190174  G4double     p0 = -2.88;
    … …  
    200184 
    201185  G4double      ra=0.80;
    202  
     186       
    203187  //JMQ 13/02/09 increase of reduced radius to lower the barrier
    204188  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
    … …  
    207191  p = p0 + p1/ec + p2/ecsq;
    208192  landa = landa0*ResidualA + landa1;
    209   a = std::pow(ResidualA,mu1);
     193  a = fG4pow->powZ(ResidualA,mu1);
    210194  mu = mu0 * a;
    211195  nu = a* (nu0+nu1*ec+nu2*ecsq); 
    212196  xnulam = nu / landa;
    213   if (xnulam > spill) xnulam=0.;
    214   if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
     197  if (xnulam > spill) { xnulam=0.; }
     198  if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
    215199 
    216200  a = -2.*p*ec + landa - nu/ecsq;
    … …  
    221205  ecut = (ecut-a) / (p+p);
    222206  ecut2 = ecut;
    223 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
    224 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
    225 //  if (cut < 0.) ecut2 = ecut - 2.;
    226   if (cut < 0.) ecut2 = ecut;
    227   elab = K * FragmentA / ResidualA;
     207  //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     208  // ecut<0 means that there is no cut with energy axis, i.e. xs is set
     209  // to 0 bellow minimum
     210  //  if (cut < 0.) ecut2 = ecut - 2.;
     211  if (cut < 0.) { ecut2 = ecut; }
     212  elab = K * FragmentA /G4double(ResidualA);
    228213  sig = 0.;
    229 
     214 
    230215  if (elab <= ec) { //start for E<Ec
    231     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     216    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
    232217  }           //end for E<Ec
    233218  else {           //start for E>Ec
    234219    sig = (landa*elab+mu+nu/elab) * signor;
    235220    geom = 0.;
    236     if (xnulam < flow || elab < etest) return sig;
     221    if (xnulam < flow || elab < etest) { return sig; }
    237222    geom = std::sqrt(theA*K);
    238223    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
    … …  
    243228 
    244229}
    245 
    246 //   ************************** end of cross sections *******************************
    247 

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationChannel.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4NeutronEvaporationChannel.cc,v 1.4 2006/06/29 20:10:37 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4NeutronEvaporationChannel.cc,v 1.5 2010/11/17 12:14:59 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Nov. 1999)
    3232//
     33// 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup
    3334
    3435#include "G4NeutronEvaporationChannel.hh"
    3536
     37G4NeutronEvaporationChannel::G4NeutronEvaporationChannel()
     38: G4EvaporationChannel(1,0,"neutron",&theEvaporationProbability,&theCoulombBarrier)
     39{}
    3640
    37 const G4NeutronEvaporationChannel & G4NeutronEvaporationChannel::operator=(const G4NeutronEvaporationChannel & )
    38 {
    39     throw G4HadronicException(__FILE__, __LINE__, "G4NeutronEvaporationChannel::operator= meant to not be accessable");
    40     return *this;
    41 }
    42 
    43 G4NeutronEvaporationChannel::G4NeutronEvaporationChannel(const G4NeutronEvaporationChannel & ) : G4EvaporationChannel()
    44 {
    45     throw G4HadronicException(__FILE__, __LINE__, "G4NeutronEvaporationChannel::CopyConstructor meant to not be accessable");
    46 }
    47 
    48 G4bool G4NeutronEvaporationChannel::operator==(const G4NeutronEvaporationChannel & right) const
    49 {
    50     return (this == (G4NeutronEvaporationChannel *) &right);
    51     //  return false;
    52 }
    53 
    54 G4bool G4NeutronEvaporationChannel::operator!=(const G4NeutronEvaporationChannel & right) const
    55 {
    56     return (this != (G4NeutronEvaporationChannel *) &right);
    57     //  return true;
    58 }
    59 
     41G4NeutronEvaporationChannel::~G4NeutronEvaporationChannel()
     42{}

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationProbability.cc

    r962 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4NeutronEvaporationProbability.cc,v 1.16 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
    3033//
    31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse
    32 // cross section option
     34// Modified:
     35// 03-09-2008 J.M. Quesada for external choice of inverse cross section option
     36// 17-11-2010 V.Ivanchenko integer Z and A
    3337
    3438#include "G4NeutronEvaporationProbability.hh"
    … …  
    3640G4NeutronEvaporationProbability::G4NeutronEvaporationProbability() :
    3741    G4EvaporationProbability(1,0,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
    38 {
    39  
    40 }
    41 
    42 
    43 G4NeutronEvaporationProbability::G4NeutronEvaporationProbability(const G4NeutronEvaporationProbability &) : G4EvaporationProbability()
    44 {
    45     throw G4HadronicException(__FILE__, __LINE__, "G4NeutronEvaporationProbability::copy_constructor meant to not be accessable");
    46 }
    47 
    48 
    49 
    50 
    51 const G4NeutronEvaporationProbability & G4NeutronEvaporationProbability::
    52 operator=(const G4NeutronEvaporationProbability &)
    53 {
    54     throw G4HadronicException(__FILE__, __LINE__, "G4NeutronEvaporationProbability::operator= meant to not be accessable");
    55     return *this;
    56 }
    57 
    58 
    59 G4bool G4NeutronEvaporationProbability::operator==(const G4NeutronEvaporationProbability &) const
    60 {
    61     return false;
    62 }
    63 
    64 G4bool G4NeutronEvaporationProbability::operator!=(const G4NeutronEvaporationProbability &) const
    65 {
    66     return true;
    67 }
    68 
    69  G4double G4NeutronEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
    70   { return 0.76+2.2/std::pow(static_cast<G4double>(fragment.GetA()-GetA()),1.0/3.0);}
    71 
     42{}
     43
     44G4NeutronEvaporationProbability::~G4NeutronEvaporationProbability()
     45{}
     46
     47G4double G4NeutronEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
     48{
     49  return 0.76+2.2/fG4pow->Z13(fragment.GetA_asInt()-GetA());
     50}
    7251       
    73   G4double G4NeutronEvaporationProbability::CalcBetaParam(const G4Fragment &  fragment)
    74   { return (2.12/std::pow(static_cast<G4double>(fragment.GetA()-GetA()),2.0/3.0) - 0.05)*MeV/
    75       CalcAlphaParam(fragment); }
    76 
     52G4double G4NeutronEvaporationProbability::CalcBetaParam(const G4Fragment &  fragment)
     53{
     54  return (2.12/fG4pow->Z23(fragment.GetA_asInt()-GetA()) - 0.05)*MeV/
     55    CalcAlphaParam(fragment);
     56}
    7757
    7858////////////////////////////////////////////////////////////////////////////////////
    … …  
    8262//OPT=3,4 Kalbach's parameterization
    8363//
    84  G4double G4NeutronEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
     64G4double
     65G4NeutronEvaporationProbability::CrossSection(const  G4Fragment & fragment, G4double K)
    8566{
    8667  theA=GetA();
    8768  theZ=GetZ();
    88   ResidualA=fragment.GetA()-theA;
    89   ResidualZ=fragment.GetZ()-theZ;
    90  
    91   ResidualAthrd=std::pow(ResidualA,0.33333);
    92   FragmentA=fragment.GetA();
    93   FragmentAthrd=std::pow(FragmentA,0.33333);
    94 
     69  ResidualA=fragment.GetA_asInt()-theA;
     70  ResidualZ=fragment.GetZ_asInt()-theZ;
     71 
     72  ResidualAthrd=fG4pow->Z13(ResidualA);
     73  FragmentA=fragment.GetA_asInt();
     74  FragmentAthrd=fG4pow->Z13(FragmentA);
    9575
    9676  if (OPTxs==0) {std::ostringstream errOs;
    … …  
    10787  }
    10888}
    109 
    110 
    111 
    112 
    113 
    114 
    115 //********************* OPT=1,2 : Chatterjee's cross section ************************
     89 
     90//********************* OPT=1,2 : Chatterjee's cross section ***************
    11691//(fitting to cross section from Bechetti & Greenles OM potential)
    11792
    118 G4double G4NeutronEvaporationProbability::GetOpt12(const  G4double K)
     93G4double G4NeutronEvaporationProbability::GetOpt12(G4double K)
    11994{
    120 
    12195  G4double Kc=K;
    12296
    123 // JMQ  xsec is set constat above limit of validity
    124   if (K>50) Kc=50;
     97  // Pramana (Bechetti & Greenles) for neutrons is chosen
     98
     99  // JMQ  xsec is set constat above limit of validity
     100  if (K > 50*MeV) { Kc = 50*MeV; }
    125101
    126102  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs;
    … …  
    143119    errOs <<"  xsec("<<Kc<<" MeV) ="<<xs <<G4endl;
    144120    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
    145   }
     121              }
    146122  return xs;
    147123}
    148124
    149 
    150125// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
    151 G4double G4NeutronEvaporationProbability::GetOpt34(const  G4double K)
     126G4double G4NeutronEvaporationProbability::GetOpt34(G4double K)
    152127{
    153 
    154128  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2;
    155129  G4double p, p0, p1, p2;
    … …  
    157131  G4double b,ecut,cut,ecut2,geom,elab;
    158132
    159 //safety initialization
     133  //safety initialization
    160134  landa0=0;
    161135  landa1=0;
    … …  
    169143  p2=0.;
    170144
    171 
    172145  flow = 1.e-18;
    173146  spill= 1.0e+18;
    174147
    175  
    176 
    177 // PRECO xs for neutrons is choosen
    178 
     148  // PRECO xs for neutrons is choosen
    179149  p0 = -312.;
    180150  p1= 0.;
    … …  
    188158  nu2 = 1280.8;
    189159
    190   if (ResidualA < 40.) signor=0.7+ResidualA*0.0075;
    191   if (ResidualA > 210.) signor = 1. + (ResidualA-210.)/250.;
     160  if (ResidualA < 40)  { signor =0.7 + ResidualA*0.0075; }
     161  if (ResidualA > 210) { signor = 1. + (ResidualA-210)/250.; }
    192162  landa = landa0/ResidualAthrd + landa1;
    193163  mu = mu0*ResidualAthrd + mu1*ResidualAthrd*ResidualAthrd;
    … …  
    195165
    196166  // JMQ very low energy behaviour corrected (problem  for A (apprx.)>60)
    197   if (nu < 0.)nu=-nu;
     167  if (nu < 0.) { nu=-nu; }
    198168
    199169  ec = 0.5;
    … …  
    202172  xnulam = 1.;
    203173  etest = 32.;
    204 //          ** etest is the energy above which the rxn cross section is
    205 //          ** compared with the geometrical limit and the max taken.
    206 //          ** xnulam here is a dummy value to be used later.
    207 
    208 
     174  //          ** etest is the energy above which the rxn cross section is
     175  //          ** compared with the geometrical limit and the max taken.
     176  //          ** xnulam here is a dummy value to be used later.
    209177
    210178  a = -2.*p*ec + landa - nu/ecsq;
    … …  
    212180  ecut = 0.;
    213181  cut = a*a - 4.*p*b;
    214   if (cut > 0.) ecut = std::sqrt(cut);
     182  if (cut > 0.) { ecut = std::sqrt(cut); }
    215183  ecut = (ecut-a) / (p+p);
    216184  ecut2 = ecut;
    217   if (cut < 0.) ecut2 = ecut - 2.;
    218   elab = K * FragmentA / ResidualA;
     185  if (cut < 0.) { ecut2 = ecut - 2.; }
     186  elab = K * FragmentA / G4double(ResidualA);
    219187  sig = 0.;
    220188  if (elab <= ec) { //start for E<Ec
    221     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;   
     189    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
    222190  }              //end for E<Ec
    223191  else {           //start for  E>Ec
    224192    sig = (landa*elab+mu+nu/elab) * signor;
    225193    geom = 0.;
    226     if (xnulam < flow || elab < etest) return sig;
     194    if (xnulam < flow || elab < etest) { return sig; }
    227195    geom = std::sqrt(theA*K);
    228196    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
    229197    geom = 31.416 * geom * geom;
    230198    sig = std::max(geom,sig);
     199
    231200  }
    232   return sig;}
    233 
    234 //   ************************** end of cross sections *******************************
    235 
     201  return sig;
     202}
     203

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationChannel.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4ProtonEvaporationChannel.cc,v 1.4 2006/06/29 20:10:41 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4ProtonEvaporationChannel.cc,v 1.5 2010/11/17 12:14:59 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Nov. 1999)
    3232//
     33// 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup
    3334
    3435#include "G4ProtonEvaporationChannel.hh"
    3536
     37G4ProtonEvaporationChannel::G4ProtonEvaporationChannel()
     38: G4EvaporationChannel(1,1,"proton",&theEvaporationProbability,&theCoulombBarrier)
     39{}
    3640
    37 const G4ProtonEvaporationChannel & G4ProtonEvaporationChannel::operator=(const G4ProtonEvaporationChannel & )
    38 {
    39     throw G4HadronicException(__FILE__, __LINE__, "G4ProtonEvaporationChannel::operator= meant to not be accessable");
    40     return *this;
    41 }
    42 
    43 G4ProtonEvaporationChannel::G4ProtonEvaporationChannel(const G4ProtonEvaporationChannel & ) : G4EvaporationChannel()
    44 {
    45     throw G4HadronicException(__FILE__, __LINE__, "G4ProtonEvaporationChannel::CopyConstructor meant to not be accessable");
    46 }
    47 
    48 G4bool G4ProtonEvaporationChannel::operator==(const G4ProtonEvaporationChannel & right) const
    49 {
    50     return (this == (G4ProtonEvaporationChannel *) &right);
    51     //  return false;
    52 }
    53 
    54 G4bool G4ProtonEvaporationChannel::operator!=(const G4ProtonEvaporationChannel & right) const
    55 {
    56     return (this != (G4ProtonEvaporationChannel *) &right);
    57     //  return true;
    58 }
    59 
    60 
     41G4ProtonEvaporationChannel::~G4ProtonEvaporationChannel()
     42{}

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationProbability.cc

    r1315 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4ProtonEvaporationProbability.cc,v 1.17 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
    3033//
    31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse
    32 // cross section option
     34// Modified:
     35// 03-09-2008 J.M. Quesada for external choice of inverse cross section option
     36// 17-11-2010 V.Ivanchenko integer Z and A
    3337
    3438#include "G4ProtonEvaporationProbability.hh"
    … …  
    3640G4ProtonEvaporationProbability::G4ProtonEvaporationProbability() :
    3741    G4EvaporationProbability(1,1,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
    38 {
    39    
    40 }
    41 
    42 G4ProtonEvaporationProbability::G4ProtonEvaporationProbability(const G4ProtonEvaporationProbability &) : G4EvaporationProbability()
    43 {
    44     throw G4HadronicException(__FILE__, __LINE__, "G4ProtonEvaporationProbability::copy_constructor meant to not be accessable");
    45 }
    46 
    47 const G4ProtonEvaporationProbability & G4ProtonEvaporationProbability::
    48 operator=(const G4ProtonEvaporationProbability &)
    49 {
    50     throw G4HadronicException(__FILE__, __LINE__, "G4ProtonEvaporationProbability::operator= meant to not be accessable");
    51     return *this;
    52 }
    53 
    54 
    55 G4bool G4ProtonEvaporationProbability::operator==(const G4ProtonEvaporationProbability &) const
    56 {
    57     return false;
    58 }
    59 
    60 G4bool G4ProtonEvaporationProbability::operator!=(const G4ProtonEvaporationProbability &) const
    61 {
    62     return true;
    63 }
    64 
    65   G4double G4ProtonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
    66   { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
    67        
    68   G4double G4ProtonEvaporationProbability::CalcBetaParam(const G4Fragment & ) 
     42{}
     43
     44G4ProtonEvaporationProbability::~G4ProtonEvaporationProbability()
     45{}
     46
     47G4double G4ProtonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
     48  { return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());}
     49       
     50G4double G4ProtonEvaporationProbability::CalcBetaParam(const G4Fragment & ) 
    6951  { return 0.0; }
    7052
    71   G4double G4ProtonEvaporationProbability::CCoeficient(const G4double aZ)
    72 {
    73     // Data comes from
    74     // Dostrovsky, Fraenkel and Friedlander
    75     // Physical Review, vol 116, num. 3 1959
    76     //
    77     // const G4int size = 5;
    78     // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
    79     // G4double Cp[5] = { 0.50, 0.28, 0.20, 0.15, 0.10};
    80     G4double C = 0.0;
    81        
    82     if (aZ >= 70) {
    83         C = 0.10;
    84     } else {
    85         C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
    86     }
    87        
    88     return C;
     53G4double G4ProtonEvaporationProbability::CCoeficient(G4int aZ)
     54{
     55  // Data comes from
     56  // Dostrovsky, Fraenkel and Friedlander
     57  // Physical Review, vol 116, num. 3 1959
     58  //
     59  // const G4int size = 5;
     60  // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
     61  // G4double Cp[5] = { 0.50, 0.28, 0.20, 0.15, 0.10};
     62  G4double C = 0.0;
     63       
     64  if (aZ >= 70) {
     65    C = 0.10;
     66  } else {
     67    C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
     68  }
     69       
     70  return C;
    8971       
    9072}
    … …  
    9779//OPT=3 Kalbach's parameterization
    9880//
    99 G4double G4ProtonEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
    100 {
    101 //  G4cout<<" In G4ProtonEVaporationProbability OPTxs="<<OPTxs<<G4endl;
    102 //  G4cout<<" In G4ProtonEVaporationProbability useSICB="<<useSICB<<G4endl;
     81G4double
     82G4ProtonEvaporationProbability::CrossSection(const  G4Fragment & fragment, G4double K)
     83{
     84  //  G4cout<<" In G4ProtonEVaporationProbability OPTxs="<<OPTxs<<G4endl;
     85  //  G4cout<<" In G4ProtonEVaporationProbability useSICB="<<useSICB<<G4endl;
    10386
    10487  theA=GetA();
    10588  theZ=GetZ();
    106   ResidualA=fragment.GetA()-theA;
    107   ResidualZ=fragment.GetZ()-theZ;
    108  
    109   ResidualAthrd=std::pow(ResidualA,0.33333);
    110   FragmentA=fragment.GetA();
    111   FragmentAthrd=std::pow(FragmentA,0.33333);
     89  ResidualA=fragment.GetA_asInt()-theA;
     90  ResidualZ=fragment.GetZ_asInt()-theZ;
     91 
     92  ResidualAthrd=fG4pow->Z13(ResidualA);
     93  FragmentA=fragment.GetA_asInt();
     94  FragmentAthrd=fG4pow->Z13(FragmentA);
     95
    11296  U=fragment.GetExcitationEnergy();
    11397
    … …  
    126110  }
    127111}
    128 //********************* OPT=1 : Chatterjee's cross section ************************
     112
     113//********************* OPT=1 : Chatterjee's cross section *********************
    129114//(fitting to cross section from Bechetti & Greenles OM potential)
    130115
    131 G4double G4ProtonEvaporationProbability::GetOpt1(const  G4double K)
     116G4double G4ProtonEvaporationProbability::GetOpt1(G4double K)
    132117{
    133118  G4double Kc=K;
    134119
    135 // JMQ  xsec is set constat above limit of validity
    136   if (K>50)  Kc=50;
     120  // JMQ  xsec is set constat above limit of validity
     121  if (K > 50*MeV) { Kc = 50*MeV; }
    137122
    138123  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs;
    … …  
    154139  p = p0 + p1/Ec + p2/(Ec*Ec);
    155140  landa = landa0*ResidualA + landa1;
    156   mu = mu0*std::pow(ResidualA,mu1);
    157   nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
     141
     142  G4double resmu1 = fG4pow->powZ(ResidualA,mu1);
     143  mu = mu0*resmu1;
     144  nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
    158145  q = landa - nu/(Ec*Ec) - 2*p*Ec;
    159146  r = mu + 2*nu/Ec + p*(Ec*Ec);
    … …  
    162149  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
    163150  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
    164    if (xs <0.0) {xs=0.0;}
    165 
    166    return xs;
    167 
    168 }
    169 
    170 
    171 
    172 //************* OPT=2 : Wellisch's proton reaction cross section ************************
    173 
    174 G4double G4ProtonEvaporationProbability::GetOpt2(const  G4double K)
    175 {
    176 
    177   G4double rnpro,rnneu,eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0);
     151  if (xs <0.0) {xs=0.0;}
     152
     153  return xs;
     154}
     155
     156//************* OPT=2 : Welisch's proton reaction cross section ***************
     157
     158G4double G4ProtonEvaporationProbability::GetOpt2(G4double K)
     159{
     160
     161  G4double eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0);
    178162 
    179 //This is redundant when the Coulomb  barrier is overimposed to all cross sections
    180 //It should be kept when Coulomb barrier only imposed at OPTxs=2, this is why ..
    181 
    182          if(!useSICB && K <= theCoulombBarrier.GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U)) return xine_th=0.0;
     163  // This is redundant when the Coulomb  barrier is overimposed to all
     164  // cross sections
     165  // It should be kept when Coulomb barrier only imposed at OPTxs=2
     166
     167  if(!useSICB && K<=theCoulombBarrier.GetCoulombBarrier(ResidualA,ResidualZ,U))
     168    { return 0.0; }
    183169
    184170  eekin=K;
    185   rnpro=ResidualZ;
    186   rnneu=ResidualA-ResidualZ;
     171  G4int rnneu=ResidualA-ResidualZ;
    187172  ekin=eekin/1000;
    188 
    189173  r0=1.36*1.e-15;
    190174  fac=pi*r0*r0;
    … …  
    192176  fac1=b0*(1.-1./ResidualAthrd);
    193177  fac2=1.;
    194   if(rnneu > 1.5) fac2=std::log(rnneu);
     178  if(rnneu > 1.5) { fac2 = fG4pow->logZ(rnneu); }
    195179  xine_th= 1.e+31*fac*fac2*(1.+ResidualAthrd-fac1);
    196180  xine_th=(1.-0.15*std::exp(-ekin))*xine_th/(1.00-0.0007*ResidualA);   
    197   ff1=0.70-0.0020*ResidualA ;
    198   ff2=1.00+1/ResidualA;
    199   ff3=0.8+18/ResidualA-0.002*ResidualA;
     181  ff1=0.70-0.0020*ResidualA;
     182  ff2=1.00+1/G4double(ResidualA);
     183  ff3=0.8+18/G4double(ResidualA)-0.002*ResidualA;
    200184  fac=1.-(1./(1.+std::exp(-8.*ff1*(std::log10(ekin)+1.37*ff2))));
    201185  xine_th=xine_th*(1.+ff3*fac);
    202   ff1=1.-1/ResidualA-0.001*ResidualA;
    203   ff2=1.17-2.7/ResidualA-0.0014*ResidualA;
     186  ff1=1.-1/G4double(ResidualA)-0.001*ResidualA;
     187  ff2=1.17-2.7/G4double(ResidualA)-0.0014*ResidualA;
    204188  fac=-8.*ff1*(std::log10(ekin)+2.0*ff2);
    205189  fac=1./(1.+std::exp(fac));
    206   xine_th=xine_th*fac;               
     190  xine_th=xine_th*fac;           
    207191  if (xine_th < 0.0){
    208192    std::ostringstream errOs;
    … …  
    212196    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
    213197  }
    214 
    215198  return xine_th;
    216            
    217 }
    218 
     199}
    219200
    220201// *********** OPT=3 : Kalbach's cross sections (from PRECO code)*************
    221202G4double G4ProtonEvaporationProbability::GetOpt3(const  G4double K)
    222203{
    223 //     ** p from  becchetti and greenlees (but modified with sub-barrier
    224 //     ** correction function and xp2 changed from -449)
     204  //     ** p from  becchetti and greenlees (but modified with sub-barrier
     205  //     ** correction function and xp2 changed from -449)
    225206
    226207  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2;
    … …  
    236217  nu1 = -182.4;
    237218  nu2 = -1.872;
    238 
    239 // parameters for  proton cross section refinement
     219 
     220  // parameters for  proton cross section refinement
    240221  G4double afit,bfit,a2,b2;
    241222  afit=-0.0785656;
    … …  
    243224  a2= -0.00089076;
    244225  b2= 0.0231597; 
    245 
     226 
    246227  G4double ec,ecsq,xnulam,etest(0.),ra(0.),a,w,c,signor(1.),signor2,sig;
    247228  G4double b,ecut,cut,ecut2,geom,elab;
    248 
    249 
     229   
    250230  G4double      flow = 1.e-18;
    251231  G4double       spill= 1.e+18;
    252 
    253 
    254   if (ResidualA <= 60.)  signor = 0.92;
    255   else if (ResidualA < 100.) signor = 0.8 + ResidualA*0.002;
    256 
    257 
     232   
     233  if (ResidualA <= 60)      { signor = 0.92; }
     234  else if (ResidualA < 100) { signor = 0.8 + ResidualA*0.002; }
     235 
    258236  ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
    259237  ecsq = ec * ec;
    260238  p = p0 + p1/ec + p2/ecsq;
    261239  landa = landa0*ResidualA + landa1;
    262   a = std::pow(ResidualA,mu1);
     240  a = fG4pow->powZ(ResidualA,mu1);
    263241  mu = mu0 * a;
    264242  nu = a* (nu0+nu1*ec+nu2*ecsq);
    265  
     243  
    266244  c =std::min(3.15,ec*0.5);
    267245  w = 0.7 * c / 3.15;
    268 
     246 
    269247  xnulam = nu / landa;
    270   if (xnulam > spill) xnulam=0.;
    271   if (xnulam >= flow) etest =std::sqrt(xnulam) + 7.;
    272 
     248  if (xnulam > spill) { xnulam=0.; }
     249  if (xnulam >= flow) { etest =std::sqrt(xnulam) + 7.; }
     250 
    273251  a = -2.*p*ec + landa - nu/ecsq;
    274252  b = p*ecsq + mu + 2.*nu/ec;
    275253  ecut = 0.;
    276254  cut = a*a - 4.*p*b;
    277   if (cut > 0.) ecut = std::sqrt(cut);
     255  if (cut > 0.) { ecut = std::sqrt(cut); }
    278256  ecut = (ecut-a) / (p+p);
    279257  ecut2 = ecut;
    280 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
    281 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
    282 //  if (cut < 0.) ecut2 = ecut - 2.;
    283  if (cut < 0.) ecut2 = ecut;
    284   elab = K * FragmentA / ResidualA;
     258  //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     259  // ecut<0 means that there is no cut with energy axis, i.e. xs is set
     260  // to 0 bellow minimum
     261  //  if (cut < 0.) ecut2 = ecut - 2.;
     262  if (cut < 0.) { ecut2 = ecut; }
     263  elab = K * FragmentA /G4double(ResidualA);
    285264  sig = 0.;
    286265  if (elab <= ec) { //start for E<Ec
    287     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     266    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
     267   
    288268    signor2 = (ec-elab-c) / w;
    289269    signor2 = 1. + std::exp(signor2);
    290     sig = sig / signor2;     
    291                        }       //end for E<=Ec
    292   else            {           //start for  E>Ec
     270    sig = sig / signor2;
     271  }              //end for E<=Ec
     272  else{           //start for  E>Ec
    293273    sig = (landa*elab+mu+nu/elab) * signor;
    294274    geom = 0.;
    295 
    296     if (xnulam < flow || elab < etest)
    297      {
     275   
     276    if (xnulam < flow || elab < etest) 
     277      {
    298278        if (sig <0.0) {sig=0.0;}
    299279        return sig;
    … …  
    303283    geom = 31.416 * geom * geom;
    304284    sig = std::max(geom,sig);
    305 
     285   
    306286  }   //end for E>Ec
    307  return sig;}
    308 
    309 
    310 
    311 //   ************************** end of cross sections *******************************
    312 
     287  return sig;
     288}
     289

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationChannel.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4TritonEvaporationChannel.cc,v 1.4 2006/06/29 20:10:45 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4TritonEvaporationChannel.cc,v 1.5 2010/11/17 12:14:59 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    3131// by V. Lara (Nov. 1999)
    3232//
     33// 17-11-2010 V.Ivanchenko moved constructor and destructor to source and cleanup
    3334
    3435#include "G4TritonEvaporationChannel.hh"
    3536
     37G4TritonEvaporationChannel::G4TritonEvaporationChannel()
     38: G4EvaporationChannel(3,1,"triton",&theEvaporationProbability,&theCoulombBarrier)
     39{}
    3640
    37 const G4TritonEvaporationChannel & G4TritonEvaporationChannel::
    38 operator=(const G4TritonEvaporationChannel & )
    39 {
    40     throw G4HadronicException(__FILE__, __LINE__, "G4TritonEvaporationChannel::operator= meant to not be accessable");
    41     return *this;
    42 }
     41G4TritonEvaporationChannel::~G4TritonEvaporationChannel()
     42{}
    4343
    44 G4TritonEvaporationChannel::G4TritonEvaporationChannel(const G4TritonEvaporationChannel & ) : G4EvaporationChannel()
    45 {
    46     throw G4HadronicException(__FILE__, __LINE__, "G4TritonEvaporationChannel::CopyConstructor meant to not be accessable");
    47 }
    48 
    49 G4bool G4TritonEvaporationChannel::operator==(const G4TritonEvaporationChannel & right) const
    50 {
    51     return (this == (G4TritonEvaporationChannel *) &right);
    52 }
    53 
    54 G4bool G4TritonEvaporationChannel::operator!=(const G4TritonEvaporationChannel & right) const
    55 {
    56     return (this != (G4TritonEvaporationChannel *) &right);
    57 }

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationProbability.cc

    r1315 r1347  
    2424// ********************************************************************
    2525//
    26 //J.M. Quesada (August2008). Based on:
     26// $Id: G4TritonEvaporationProbability.cc,v 1.18 2010/11/17 11:06:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
     28//
     29// J.M. Quesada (August2008). Based on:
    2730//
    2831// Hadronic Process: Nuclear De-excitations
    2932// by V. Lara (Oct 1998)
    3033//
    31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse
    32 // cross section option
     34// Modified:
     35// 03-09-2008 J.M. Quesada for external choice of inverse cross section option
     36// 17-11-2010 V.Ivanchenko integer Z and A
    3337
    3438#include "G4TritonEvaporationProbability.hh"
    … …  
    3640G4TritonEvaporationProbability::G4TritonEvaporationProbability() :
    3741    G4EvaporationProbability(3,1,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
    38 {
    39 
    40 }
    41 
    42 G4TritonEvaporationProbability::G4TritonEvaporationProbability(const G4TritonEvaporationProbability &) : G4EvaporationProbability()
    43 {
    44     throw G4HadronicException(__FILE__, __LINE__, "G4TritonEvaporationProbability::copy_constructor meant to not be accessable");
    45 }
    46 
    47 
    48 
    49 
    50 const G4TritonEvaporationProbability & G4TritonEvaporationProbability::
    51 operator=(const G4TritonEvaporationProbability &)
    52 {
    53     throw G4HadronicException(__FILE__, __LINE__, "G4TritonEvaporationProbability::operator= meant to not be accessable");
    54     return *this;
    55 }
    56 
    57 
    58 G4bool G4TritonEvaporationProbability::operator==(const G4TritonEvaporationProbability &) const
    59 {
    60     return false;
    61 }
    62 
    63 G4bool G4TritonEvaporationProbability::operator!=(const G4TritonEvaporationProbability &) const
    64 {
    65     return true;
    66 }
    67 
    68    G4double G4TritonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
    69   { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
     42{}
     43
     44G4TritonEvaporationProbability::~G4TritonEvaporationProbability()
     45{}
     46
     47G4double G4TritonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
     48{
     49  return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());
     50}
    7051       
    71   G4double G4TritonEvaporationProbability::CalcBetaParam(const G4Fragment & )
    72   { return 0.0; }
    73 
    74 G4double G4TritonEvaporationProbability::CCoeficient(const G4double aZ)
    75 {
    76     // Data comes from
    77     // Dostrovsky, Fraenkel and Friedlander
    78     // Physical Review, vol 116, num. 3 1959
    79     //
    80     // const G4int size = 5;
    81     // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
    82     // G4double Cp[5] = { 0.50, 0.28, 0.20, 0.15, 0.10};
    83     // C for triton is equal to C for protons divided by 3
    84     G4double C = 0.0;
     52G4double G4TritonEvaporationProbability::CalcBetaParam(const G4Fragment & )
     53{
     54  return 0.0;
     55}
     56
     57G4double G4TritonEvaporationProbability::CCoeficient(G4int aZ)
     58{
     59  // Data comes from
     60  // Dostrovsky, Fraenkel and Friedlander
     61  // Physical Review, vol 116, num. 3 1959
     62  //
     63  // const G4int size = 5;
     64  // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
     65  // G4double Cp[5] = { 0.50, 0.28, 0.20, 0.15, 0.10};
     66  // C for triton is equal to C for protons divided by 3
     67  G4double C = 0.0;
    8568       
    86     if (aZ >= 70) {
    87         C = 0.10;
    88     } else {
    89         C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
    90     }
     69  if (aZ >= 70) {
     70    C = 0.10;
     71  } else {
     72    C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
     73  }
    9174       
    92     return C/3.0;
    93        
     75  return C/3.0;
    9476}
    9577
    … …  
    10082//OPT=3,4 Kalbach's parameterization
    10183//
    102 G4double G4TritonEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
     84G4double
     85G4TritonEvaporationProbability::CrossSection(const  G4Fragment & fragment, G4double K)
    10386{
    10487  theA=GetA();
    10588  theZ=GetZ();
    106   ResidualA=fragment.GetA()-theA;
    107   ResidualZ=fragment.GetZ()-theZ;
     89  ResidualA=fragment.GetA_asInt()-theA;
     90  ResidualZ=fragment.GetZ_asInt()-theZ;
    10891 
    109   ResidualAthrd=std::pow(ResidualA,0.33333);
    110   FragmentA=fragment.GetA();
    111   FragmentAthrd=std::pow(FragmentA,0.33333);
    112 
    113  
     92  ResidualAthrd=fG4pow->Z13(ResidualA);
     93  FragmentA=fragment.GetA_asInt();
     94  FragmentAthrd=fG4pow->Z13(FragmentA);
     95
    11496  if (OPTxs==0) {std::ostringstream errOs;
    11597    errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (tritons)!!" 
    … …  
    128110
    129111//
    130 //********************* OPT=1,2 : Chatterjee's cross section ************************
     112//********************* OPT=1,2 : Chatterjee's cross section *****************
    131113//(fitting to cross section from Bechetti & Greenles OM potential)
    132114
    133 G4double G4TritonEvaporationProbability::GetOpt12(const  G4double K)
    134 {
    135 
     115G4double G4TritonEvaporationProbability::GetOpt12(G4double K)
     116{
    136117  G4double Kc=K;
    137118
    138 // JMQ xsec is set constat above limit of validity
    139   if (K>50) Kc=50;
    140 
    141  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
    142 //G4double Eo(0),epsilon1(0),epsilon2(0),discri(0);
    143 
     119  // JMQ xsec is set constat above limit of validity
     120  if (K > 50*MeV) { Kc=50*MeV; }
     121
     122  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
    144123 
    145  G4double    p0 = -11.04;
    146  G4double    p1 = 619.1;
    147  G4double    p2 = -2147.;
    148  G4double    landa0 = -0.0426;
    149  G4double    landa1 = -10.33;
    150  G4double    mu0 = 601.9;
    151  G4double    mu1 = 0.37;
    152  G4double    nu0 = 583.0;
    153  G4double    nu1 = -546.2;
    154  G4double    nu2 = 1.718; 
    155  G4double    delta=1.2;           
    156 
    157  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
    158  p = p0 + p1/Ec + p2/(Ec*Ec);
    159  landa = landa0*ResidualA + landa1;
    160  mu = mu0*std::pow(ResidualA,mu1);
    161  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
    162  q = landa - nu/(Ec*Ec) - 2*p*Ec;
    163  r = mu + 2*nu/Ec + p*(Ec*Ec);
    164 
    165  ji=std::max(Kc,Ec);
    166  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
    167  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
    168  
    169  if (xs <0.0) {xs=0.0;}
    170  
    171  return xs;
    172 
     124  G4double    p0 = -11.04;
     125  G4double    p1 = 619.1;
     126  G4double    p2 = -2147.;
     127  G4double    landa0 = -0.0426;
     128  G4double    landa1 = -10.33;
     129  G4double    mu0 = 601.9;
     130  G4double    mu1 = 0.37;
     131  G4double    nu0 = 583.0;
     132  G4double    nu1 = -546.2;
     133  G4double    nu2 = 1.718; 
     134  G4double    delta=1.2;           
     135
     136  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
     137  p = p0 + p1/Ec + p2/(Ec*Ec);
     138  landa = landa0*ResidualA + landa1;
     139
     140  G4double resmu1 = fG4pow->powZ(ResidualA,mu1);
     141  mu = mu0*resmu1;
     142  nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
     143  q = landa - nu/(Ec*Ec) - 2*p*Ec;
     144  r = mu + 2*nu/Ec + p*(Ec*Ec);
     145 
     146  ji=std::max(Kc,Ec);
     147  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
     148  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
     149                 
     150  if (xs <0.0) {xs=0.0;}
     151             
     152  return xs;
    173153}
    174154
    175155// *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
    176 G4double G4TritonEvaporationProbability::GetOpt34(const  G4double K)
     156G4double G4TritonEvaporationProbability::GetOpt34(G4double K)
    177157//     ** t from o.m. of hafele, flynn et al
    178158{
    179  
    180159  G4double landa, mu, nu, p , signor(1.),sig;
    181160  G4double ec,ecsq,xnulam,etest(0.),a;
    182161  G4double b,ecut,cut,ecut2,geom,elab;
    183  
    184  
     162
    185163  G4double     flow = 1.e-18;
    186164  G4double     spill= 1.e+18;
    … …  
    205183  p = p0 + p1/ec + p2/ecsq;
    206184  landa = landa0*ResidualA + landa1;
    207   a = std::pow(ResidualA,mu1);
     185  a = fG4pow->powZ(ResidualA,mu1);
    208186  mu = mu0 * a;
    209187  nu = a* (nu0+nu1*ec+nu2*ecsq); 
    210188  xnulam = nu / landa;
    211   if (xnulam > spill) xnulam=0.;
    212   if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam);
    213   
     189  if (xnulam > spill) { xnulam=0.; }
     190  if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
     191 
    214192  a = -2.*p*ec + landa - nu/ecsq;
    215193  b = p*ecsq + mu + 2.*nu/ec;
    216194  ecut = 0.;
    217195  cut = a*a - 4.*p*b;
    218   if (cut > 0.) ecut = std::sqrt(cut);
     196  if (cut > 0.) { ecut = std::sqrt(cut); }
    219197  ecut = (ecut-a) / (p+p);
    220198  ecut2 = ecut;
    221 //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
    222 //ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
    223 //  if (cut < 0.) ecut2 = ecut - 2.;
    224   if (cut < 0.) ecut2 = ecut;
    225   elab = K * FragmentA / ResidualA;
     199  //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     200  // ecut<0 means that there is no cut with energy axis, i.e. xs is set
     201  // to 0 bellow minimum
     202  //  if (cut < 0.) ecut2 = ecut - 2.;
     203  if (cut < 0.) { ecut2 = ecut; }
     204  elab = K * FragmentA / G4double(ResidualA);
    226205  sig = 0.;
    227 
     206 
    228207  if (elab <= ec) { //start for E<Ec
    229     if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
     208    if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
    230209  }           //end for E<Ec
    231   else {           //start for E>Ec 
     210  else {           //start for E>Ec
    232211    sig = (landa*elab+mu+nu/elab) * signor;
    233212    geom = 0.;
    234     if (xnulam < flow || elab < etest) return sig;
     213    if (xnulam < flow || elab < etest) { return sig; }
    235214    geom = std::sqrt(theA*K);
    236215    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
    … …  
    239218  }           //end for E>Ec
    240219  return sig;
    241  
    242 }
    243 
    244 //   ************************** end of cross sections *******************************
    245 
     220}
     221

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4UnstableFragmentBreakUp.cc

    r1340 r1347  
    2525//
    2626// $Id: G4UnstableFragmentBreakUp.cc,v 1.8 2010/06/25 09:46:09 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2828//
    2929// -------------------------------------------------------------------

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4VEvaporation.cc

    r1340 r1347  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEvaporation.cc,v 1.7 2010/04/27 14:00:40 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     26// $Id: G4VEvaporation.cc,v 1.8 2010/10/29 17:35:03 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2828//
    2929// Hadronic Process: Nuclear De-excitations
    … …  
    3333#include "G4VEvaporation.hh"
    3434
    35 G4VEvaporation::G4VEvaporation() 
     35G4VEvaporation::G4VEvaporation():OPTxs(3),useSICB(false)
    3636{}
    3737
Note: See TracChangeset for help on using the changeset viewer.