- Timestamp:
- Apr 6, 2009, 12:30:29 PM (15 years ago)
- Location:
- trunk/source/processes/hadronic/models/de_excitation/evaporation
- Files:
-
- 34 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4AlphaEvaporationChannel.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4AlphaEvaporationChannel.hh,v 1. 3 2006/06/29 20:09:43 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4AlphaEvaporationChannel.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Nov. 1999) 32 32 // 33 34 33 35 34 #ifndef G4AlphaEvaporationChannel_h … … 62 61 private: 63 62 64 G4AlphaCoulombBarrier theCoulombBarrier; 65 63 G4AlphaCoulombBarrier theCoulombBarrier; 64 66 65 G4AlphaEvaporationProbability theEvaporationProbability; 67 66 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4AlphaEvaporationProbability.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4AlphaEvaporationProbability.hh,v 1.3 2006/06/29 20:09:45 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara ( Nov 1999)29 // by V. Lara (Oct 1998) 32 30 // 33 34 35 36 31 #ifndef G4AlphaEvaporationProbability_h 37 32 #define G4AlphaEvaporationProbability_h 1 … … 39 34 40 35 #include "G4EvaporationProbability.hh" 41 36 #include "G4AlphaCoulombBarrier.hh" 42 37 43 38 class G4AlphaEvaporationProbability : public G4EvaporationProbability … … 55 50 G4bool operator==(const G4AlphaEvaporationProbability &right) const; 56 51 G4bool operator!=(const G4AlphaEvaporationProbability &right) const; 57 52 58 53 59 54 private: 60 55 61 virtual G4double CalcAlphaParam(const G4Fragment & fragment) const 62 { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));} 63 64 virtual G4double CalcBetaParam(const G4Fragment & ) const 65 { return 0.0; } 56 virtual G4double CrossSection(const G4Fragment & fragment, const G4double K); 66 57 67 68 G4double CCoeficient(const G4double ) const; 58 G4double GetOpt0(const G4double K); 59 G4double GetOpt12(const G4double K); 60 G4double GetOpt34(const G4double K); 69 61 70 // Excitation energy levels 71 std::vector<G4double> ExcitEnergies; 72 // Spin of excitation energy levels 73 std::vector<G4int> ExcitSpins; 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 81 74 82 }; 75 83 76 84 77 85 #endif 86 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4DeuteronEvaporationChannel.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4DeuteronEvaporationChannel.hh,v 1. 3 2006/06/29 20:09:47 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4DeuteronEvaporationChannel.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Nov. 1999) 32 32 // 33 34 33 35 34 #ifndef G4DeuteronEvaporationChannel_h -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4DeuteronEvaporationProbability.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4DeuteronEvaporationProbability.hh,v 1.3 2006/06/29 20:09:49 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara ( Nov 1999)29 // by V. Lara (Oct 1998) 32 30 // 33 34 35 31 36 32 #ifndef G4DeuteronEvaporationProbability_h … … 39 35 40 36 #include "G4EvaporationProbability.hh" 41 37 #include "G4DeuteronCoulombBarrier.hh" 42 38 43 39 class G4DeuteronEvaporationProbability : public G4EvaporationProbability … … 55 51 G4bool operator==(const G4DeuteronEvaporationProbability &right) const; 56 52 G4bool operator!=(const G4DeuteronEvaporationProbability &right) const; 57 53 58 54 59 55 private: 60 56 61 virtual G4double CalcAlphaParam(const G4Fragment & fragment) const 62 { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));} 63 64 virtual G4double CalcBetaParam(const G4Fragment & ) const 65 { return 0.0; } 57 virtual G4double CrossSection(const G4Fragment & fragment, const G4double K); 66 58 67 68 G4double CCoeficient(const G4double aZ) const; 59 G4double GetOpt0(const G4double K); 60 G4double GetOpt12(const G4double K); 61 G4double GetOpt34(const G4double K); 69 62 70 // Excitation energy levels 71 std::vector<G4double> ExcitEnergies; 72 // Spin of excitation energy levels 73 std::vector<G4int> ExcitSpins; 63 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 71 72 G4DeuteronCoulombBarrier theCoulombBarrier; 73 74 G4double ResidualA; 75 G4double ResidualZ; 76 G4double theA; 77 G4double theZ; 78 G4double ResidualAthrd; 79 G4double FragmentA; 80 G4double FragmentAthrd; 81 74 82 75 83 }; 84 85 76 86 #endif 87 88 89 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4Evaporation.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4Evaporation.hh,v 1. 4 2007/02/14 13:37:49ahoward Exp $28 // GEANT4 tag $Name: $27 // $Id: G4Evaporation.hh,v 1.5 2008/09/19 13:32:54 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations … … 69 69 void SetDefaultChannel(); 70 70 void SetGEMChannel(); 71 72 private:73 71 74 72 #ifdef debug -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4EvaporationChannel.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4EvaporationChannel.hh,v 1.3 2006/06/29 20:09:53 gunter Exp $ 28 // GEANT4 tag $Name: $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations … … 39 37 #include "G4VEmissionProbability.hh" 40 38 #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" 41 45 #include "G4VLevelDensityParameter.hh" 42 46 #include "G4VCoulombBarrier.hh" … … 51 55 { 52 56 public: 53 // Available constructors 54 G4EvaporationChannel(const G4int theA, const G4int theZ, 55 G4VEmissionProbability * aEmissionStrategy, 56 G4VCoulombBarrier * aCoulombBarrier); 57 // constructor 58 57 59 58 60 G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String & aName, 59 61 G4VEmissionProbability * aEmissionStrategy, 60 G4VCoulombBarrier * aCoulombBarrier); 61 62 G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String * aName, 63 G4VEmissionProbability * aEmissionStrategy, 64 G4VCoulombBarrier * aCoulombBarrier); 65 62 G4VCoulombBarrier * aCoulombBarrier); 66 63 public: 67 64 // destructor … … 70 67 void SetEmissionStrategy(G4VEmissionProbability * aEmissionStrategy) 71 68 {theEvaporationProbabilityPtr = aEmissionStrategy;} 72 69 73 70 void SetCoulombBarrierStrategy(G4VCoulombBarrier * aCoulombBarrier) 74 71 {theCoulombBarrierPtr = aCoulombBarrier;} 72 73 75 74 76 75 protected: … … 93 92 94 93 G4FragmentVector * BreakUp(const G4Fragment & theNucleus); 95 96 // inline void SetEmissionStrategy(G4VEmissionProbability * aStrategy)97 // {98 // if (MyOwnEvaporationProbability) delete theEvaporationProbabilityPtr;99 // theEvaporationProbabilityPtr = aStrategy;100 // MyOwnEvaporationProbability = false;101 // }102 94 103 104 inline void SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity) 105 { 106 if (MyOwnLevelDensity) delete theLevelDensityPtr; 107 theLevelDensityPtr = aLevelDensity; 108 MyOwnLevelDensity = false; 109 } 110 95 // void SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity); 96 111 97 public: 112 98 … … 119 105 { return MaximalKineticEnergy; } 120 106 121 // ----------------------122 123 107 private: 124 108 … … 130 114 131 115 // Samples fragment kinetic energy. 132 G4double CalcKineticEnergy(void);116 G4double GetKineticEnergy(const G4Fragment & aFragment); 133 117 134 118 // This has to be removed and put in Random Generator … … 142 126 // They are intializated at object creation (constructor) time. 143 127 144 // Atomic Number 145 G4int A;128 // Atomic Number of ejectile 129 G4int theA; 146 130 147 // Charge 148 G4int Z;131 // Charge of ejectile 132 G4int theZ; 149 133 150 134 151 // For evaporation probability calcual tion135 // For evaporation probability calcualation 152 136 G4VEmissionProbability * theEvaporationProbabilityPtr; 153 137 154 138 // For Level Density calculation 155 G4bool MyOwnLevelDensity;139 // G4bool MyOwnLevelDensity; 156 140 G4VLevelDensityParameter * theLevelDensityPtr; 157 141 142 158 143 // For Coulomb Barrier calculation 159 144 G4VCoulombBarrier * theCoulombBarrierPtr; 160 145 G4double CoulombBarrier; 161 146 147 162 148 //--------------------------------------------------- 163 149 … … 166 152 // the atomic number, charge and excitation energy of nucleus. 167 153 168 // Residual AtomicNumber169 G4int AResidual;154 // Residual Mass Number 155 G4int ResidualA; 170 156 171 157 // Residual Charge 172 G4int ZResidual; 173 174 // // Binding Energy 175 // G4double BindingEnergy; 176 177 // // Level Density Parameter 178 // G4double LevelDensityParameter; 158 G4int ResidualZ; 179 159 180 160 // Emission Probability … … 184 164 // Maximal Kinetic Energy that can be carried by fragment 185 165 G4double MaximalKineticEnergy; 166 167 186 168 }; 187 169 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4EvaporationFactory.hh
r819 r962 26 26 // 27 27 // $Id: G4EvaporationFactory.hh,v 1.3 2006/06/29 20:09:55 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4EvaporationProbability.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4EvaporationProbability.hh,v 1.3 2006/06/29 20:09:57 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara (Oct 1998) 29 // by V. Lara (Oct 1998) 32 30 // 33 34 35 36 31 #ifndef G4EvaporationProbability_h 37 32 #define G4EvaporationProbability_h 1 … … 41 36 #include "G4VLevelDensityParameter.hh" 42 37 #include "G4EvaporationLevelDensityParameter.hh" 38 #include "G4VCoulombBarrier.hh" 39 #include "G4CoulombBarrier.hh" 43 40 44 41 … … 47 44 public: 48 45 // Only available constructor 49 G4EvaporationProbability(const G4int anA, const G4int aZ, const G4double aGamma ) :46 G4EvaporationProbability(const G4int anA, const G4int aZ, const G4double aGamma,G4VCoulombBarrier * aCoulombBarrier) : 50 47 theA(anA), 51 48 theZ(aZ), 52 Gamma(aGamma) 49 Gamma(aGamma) 50 , theCoulombBarrierptr(aCoulombBarrier) 53 51 { 54 52 theEvapLDPptr = new G4EvaporationLevelDensityParameter; 53 54 55 55 } 56 56 … … 58 58 { 59 59 if (theEvapLDPptr != 0) delete theEvapLDPptr; 60 60 61 } 61 62 … … 67 68 68 69 protected: 69 70 void SetExcitationEnergiesPtr(std::vector<G4double> * anExcitationEnergiesPtr)71 {ExcitationEnergies = anExcitationEnergiesPtr;}72 73 void SetExcitationSpinsPtr(std::vector<G4int> * anExcitationSpinsPtr)74 {ExcitationSpins = anExcitationSpinsPtr;}75 76 70 77 71 // Default constructor 78 72 G4EvaporationProbability() {} 73 79 74 private: 80 75 // Copy constructor … … 86 81 87 82 public: 88 G4double EmissionProbability(const G4Fragment & fragment, const G4double anEnergy); 83 84 G4double ProbabilityDistributionFunction( const G4Fragment & aFragment, const G4double K); 85 86 G4double EmissionProbability(const G4Fragment & fragment, const G4double anEnergy); 89 87 90 88 private: 91 89 92 G4double CalcProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy); 93 virtual G4double CCoeficient(const G4double ) const {return 0.0;}; 90 G4double CalculateProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy ); 94 91 95 virtual G4double CalcAlphaParam(const G4Fragment & ) const {return 1.0;} 96 virtual G4double CalcBetaParam(const G4Fragment & ) const {return 1.0;} 92 G4double IntegrateEmissionProbability(const G4Fragment & aFragment, const G4double & Low, const G4double & Up ); 93 94 protected: 95 96 virtual G4double CrossSection( const G4Fragment & fragment, const G4double K )= 0; 97 98 virtual G4double CalcAlphaParam(const G4Fragment & fragment)=0 ; 99 100 virtual G4double CalcBetaParam(const G4Fragment & fragment)=0 ; 101 102 private: 97 103 98 104 // Data Members … … 107 113 G4double Gamma; 108 114 109 // Discrete Excitation Energies 110 std::vector<G4double> * ExcitationEnergies;115 //The Coulomb Barrier 116 G4VCoulombBarrier * theCoulombBarrierptr; 111 117 112 //113 std::vector<G4int> * ExcitationSpins;114 118 115 119 }; 116 120 117 121 122 123 118 124 #endif -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4He3EvaporationChannel.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4He3EvaporationChannel.hh,v 1. 3 2006/06/29 20:09:59 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4He3EvaporationChannel.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Nov. 1999) 32 32 // 33 34 33 35 34 #ifndef G4He3EvaporationChannel_h … … 61 60 private: 62 61 63 G4He3CoulombBarrier theCoulombBarrier;62 G4He3CoulombBarrier theCoulombBarrier; 64 63 65 64 G4He3EvaporationProbability theEvaporationProbability; -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4He3EvaporationProbability.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4He3EvaporationProbability.hh,v 1.3 2006/06/29 20:10:01 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara ( Nov 1999)29 // by V. Lara (Oct 1998) 32 30 // 33 34 35 36 31 #ifndef G4He3EvaporationProbability_h 37 32 #define G4He3EvaporationProbability_h 1 … … 39 34 40 35 #include "G4EvaporationProbability.hh" 41 36 #include "G4He3CoulombBarrier.hh" 42 37 43 38 class G4He3EvaporationProbability : public G4EvaporationProbability … … 55 50 G4bool operator==(const G4He3EvaporationProbability &right) const; 56 51 G4bool operator!=(const G4He3EvaporationProbability &right) const; 57 52 58 53 59 54 private: 60 55 61 virtual G4double CalcAlphaParam(const G4Fragment & fragment) const 62 { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));} 63 64 virtual G4double CalcBetaParam(const G4Fragment & ) const 65 { return 0.0; } 56 virtual G4double CrossSection(const G4Fragment & fragment, const G4double K); 66 57 67 68 G4double CCoeficient(const G4double aZ) const; 58 G4double GetOpt0(const G4double K); 59 G4double GetOpt12(const G4double K); 60 G4double GetOpt34(const G4double K); 69 61 70 // Excitation energy levels 71 std::vector<G4double> ExcitEnergies; 72 // Spin of excitation energy levels 73 std::vector<G4int> ExcitSpins; 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 G4He3CoulombBarrier theCoulombBarrier; 72 73 G4double ResidualA; 74 G4double ResidualZ; 75 G4double theA; 76 G4double theZ; 77 G4double ResidualAthrd; 78 G4double FragmentA; 79 G4double FragmentAthrd; 80 81 74 82 }; 75 83 76 84 77 85 #endif 86 87 88 89 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4NeutronEvaporationChannel.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4NeutronEvaporationChannel.hh,v 1. 3 2006/06/29 20:10:03 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4NeutronEvaporationChannel.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Nov. 1999) 32 32 // 33 34 33 35 34 #ifndef G4NeutronEvaporationChannel_h … … 61 60 private: 62 61 63 62 G4NeutronCoulombBarrier theCoulombBarrier; 64 63 65 64 G4NeutronEvaporationProbability theEvaporationProbability; -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4NeutronEvaporationProbability.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4NeutronEvaporationProbability.hh,v 1.4 2006/06/29 20:10:05 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara ( Nov 1999)29 // by V. Lara (Oct 1998) 32 30 // 33 34 35 31 36 32 #ifndef G4NeutronEvaporationProbability_h … … 39 35 40 36 #include "G4EvaporationProbability.hh" 41 37 #include "G4NeutronCoulombBarrier.hh" 42 38 43 39 class G4NeutronEvaporationProbability : public G4EvaporationProbability 44 40 { 45 41 public: 46 // Only available constructor42 47 43 G4NeutronEvaporationProbability(); 48 44 49 45 ~G4NeutronEvaporationProbability() {} 50 46 private: 51 // Copy constructor47 52 48 G4NeutronEvaporationProbability(const G4NeutronEvaporationProbability &right); 53 49 … … 55 51 G4bool operator==(const G4NeutronEvaporationProbability &right) const; 56 52 G4bool operator!=(const G4NeutronEvaporationProbability &right) const; 57 58 53 59 54 private: 60 55 61 virtual G4double CalcAlphaParam(const G4Fragment & fragment) const 62 { return 0.76+2.2/std::pow(static_cast<G4double>(fragment.GetA()-GetA()),1.0/3.0);} 63 64 virtual G4double CalcBetaParam(const G4Fragment & fragment) const 65 { return (2.12/std::pow(static_cast<G4double>(fragment.GetA()-GetA()),2.0/3.0) - 0.05)*MeV/ 66 CalcAlphaParam(fragment); } 56 virtual G4double CrossSection(const G4Fragment & fragment, const G4double K); 67 57 68 // Excitation energy levels 69 std::vector<G4double> ExcitEnergies; 70 // Spin of excitation energy levels 71 std::vector<G4int> ExcitSpins; 58 G4double GetOpt12(const G4double K); 59 G4double GetOpt34(const G4double K); 60 61 virtual G4double CalcAlphaParam(const G4Fragment & fragment); 62 63 virtual G4double CalcBetaParam(const G4Fragment & fragment); 64 65 66 //data members 67 68 G4NeutronCoulombBarrier theCoulombBarrier; 69 70 G4double ResidualA; 71 G4double ResidualZ; 72 G4double theA; 73 G4double theZ; 74 G4double ResidualAthrd; 75 G4double FragmentA; 76 G4double FragmentAthrd; 77 78 79 72 80 73 81 }; -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4ProtonEvaporationChannel.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4ProtonEvaporationChannel.hh,v 1. 3 2006/06/29 20:10:07 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4ProtonEvaporationChannel.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Nov. 1999) 32 32 // 33 34 33 35 34 #ifndef G4ProtonEvaporationChannel_h … … 62 61 private: 63 62 63 G4ProtonCoulombBarrier theCoulombBarrier; 64 64 65 G4ProtonEvaporationProbability theEvaporationProbability; 65 66 66 G4ProtonCoulombBarrier theCoulombBarrier; 67 67 68 }; 68 69 #endif -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4ProtonEvaporationProbability.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4ProtonEvaporationProbability.hh,v 1.3 2006/06/29 20:10:09 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara ( Nov 1999)29 // by V. Lara (Oct 1998) 32 30 // 33 34 35 31 36 32 #ifndef G4ProtonEvaporationProbability_h … … 39 35 40 36 #include "G4EvaporationProbability.hh" 41 37 #include "G4ProtonCoulombBarrier.hh" 42 38 43 39 class G4ProtonEvaporationProbability : public G4EvaporationProbability … … 55 51 G4bool operator==(const G4ProtonEvaporationProbability &right) const; 56 52 G4bool operator!=(const G4ProtonEvaporationProbability &right) const; 57 53 58 54 59 55 private: 60 56 61 virtual G4double CalcAlphaParam(const G4Fragment & fragment) const 62 { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));} 63 64 virtual G4double CalcBetaParam(const G4Fragment & ) const 65 { return 0.0; } 57 virtual G4double CrossSection(const G4Fragment & fragment, const G4double K); 66 58 67 68 G4double CCoeficient(const G4double aZ) const; 59 G4double GetOpt0(const G4double K); 60 G4double GetOpt1(const G4double K); 61 G4double GetOpt2(const G4double K); 62 G4double GetOpt3(const G4double K); 69 63 70 // Excitation energy levels 71 std::vector<G4double> ExcitEnergies; 72 // Spin of excitation energy levels 73 std::vector<G4int> ExcitSpins; 64 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 72 73 G4ProtonCoulombBarrier theCoulombBarrier; 74 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 74 86 75 87 }; 88 89 76 90 #endif -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4TritonEvaporationChannel.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4TritonEvaporationChannel.hh,v 1. 3 2006/06/29 20:10:11 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4TritonEvaporationChannel.hh,v 1.8 2008/09/19 13:32:54 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Nov. 1999) 32 32 // 33 34 33 35 34 #ifndef G4TritonEvaporationChannel_h -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4TritonEvaporationProbability.hh
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4TritonEvaporationProbability.hh,v 1.3 2006/06/29 20:10:13 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara ( Nov 1999)29 // by V. Lara (Oct 1998) 32 30 // 33 34 35 36 31 #ifndef G4TritonEvaporationProbability_h 37 32 #define G4TritonEvaporationProbability_h 1 … … 39 34 40 35 #include "G4EvaporationProbability.hh" 41 36 #include "G4TritonCoulombBarrier.hh" 42 37 43 38 class G4TritonEvaporationProbability : public G4EvaporationProbability … … 55 50 G4bool operator==(const G4TritonEvaporationProbability &right) const; 56 51 G4bool operator!=(const G4TritonEvaporationProbability &right) const; 57 52 58 53 59 54 private: 60 55 61 virtual G4double CalcAlphaParam(const G4Fragment & fragment) const 62 { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));} 63 64 virtual G4double CalcBetaParam(const G4Fragment & ) const 65 { return 0.0; } 56 virtual G4double CrossSection(const G4Fragment & fragment, const G4double K); 66 57 67 G4double CCoeficient(const G4double aZ) const; 58 G4double GetOpt0(const G4double K); 59 G4double GetOpt12(const G4double K); 60 G4double GetOpt34(const G4double K); 68 61 69 // Excitation energy levels 70 std::vector<G4double> ExcitEnergies; 71 // Spin of excitation energy levels 72 std::vector<G4int> ExcitSpins; 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 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 73 82 74 83 }; … … 76 85 77 86 #endif 87 88 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/include/G4VEvaporation.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4VEvaporation.hh,v 1. 3 2006/06/29 20:10:15 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4VEvaporation.hh,v 1.4 2008/09/19 13:32:54 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations 31 31 // by V. Lara (Oct 1998) written from G4Evaporation.hh (May 1998) 32 32 // 33 33 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 34 // cross section option 35 // JMQ (06 September 2008) Also external choices have been added for 36 // superimposed Coulomb barrier (if useSICBis set true, by default is false) 34 37 35 38 … … 56 59 virtual G4FragmentVector * BreakItUp(const G4Fragment &theNucleus) = 0; 57 60 61 // for inverse cross section choice 62 inline void SetOPTxs(G4int opt) { OPTxs = opt;} 63 // for superimposed Coulomb Barrier for inverse cross sections 64 inline void UseSICB(G4bool use) { useSICB = use; } 65 protected: 66 G4int OPTxs; 67 G4bool useSICB; 58 68 59 69 }; -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationChannel.cc
r819 r962 26 26 // 27 27 // $Id: G4AlphaEvaporationChannel.cc,v 1.4 2006/06/29 20:10:17 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4AlphaEvaporationProbability.cc,v 1.4 2006/06/29 20:10:19 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara (Nov 1999) 32 // 33 29 // by V. Lara (Oct 1998) 30 // 31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 32 // cross section option 34 33 35 34 #include "G4AlphaEvaporationProbability.hh" 36 35 37 36 G4AlphaEvaporationProbability::G4AlphaEvaporationProbability() : 38 G4EvaporationProbability(4,2,4) // A,Z,Gamma 39 { 40 // const G4int NumExcitedStates = 31+1; 41 std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1; 42 std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1; 43 ExcitEnergies.reserve(NumExcitedStatesEnergy); 44 ExcitSpins.reserve(NumExcitedStatesSpin); 45 ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0); 46 ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0); 47 // for (G4int i = 0; i < NumExcitedStates; i++) { 48 // ExcitEnergies(i) = 0.0; 49 // ExcitSpins(i) = 0; 50 // } 51 52 ExcitEnergies[18] = 7.98*MeV; 53 ExcitEnergies[20] = 6.90*MeV; 54 ExcitEnergies[25] = 5.83*MeV; 55 ExcitEnergies[26] = 8.57*MeV; 56 ExcitEnergies[31] = 5.33*MeV; 57 58 ExcitSpins[18] = 4; 59 ExcitSpins[20] = 6; 60 ExcitSpins[25] = 7; 61 ExcitSpins[26] = 4; 62 ExcitSpins[31] = 13; 63 64 SetExcitationEnergiesPtr(&ExcitEnergies); 65 SetExcitationSpinsPtr(&ExcitSpins); 37 G4EvaporationProbability(4,2,1,&theCoulombBarrier) // A,Z,Gamma,&theCoumlombBarrier 38 { 39 66 40 } 67 41 … … 93 67 } 94 68 95 G4double G4AlphaEvaporationProbability::CCoeficient(const G4double aZ) const 69 70 G4double G4AlphaEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 71 { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));} 72 73 G4double G4AlphaEvaporationProbability::CalcBetaParam(const G4Fragment &) 74 { return 0.0; } 75 76 G4double G4AlphaEvaporationProbability::CCoeficient(const G4double aZ) 96 77 { 97 78 // Data comes from … … 116 97 return C; 117 98 } 99 100 /////////////////////////////////////////////////////////////////////////////////// 101 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 102 //OPT=0 Dostrovski's parameterization 103 //OPT=1,2 Chatterjee's paramaterization 104 //OPT=3,4 Kalbach's parameterization 105 // 106 G4double G4AlphaEvaporationProbability::CrossSection(const G4Fragment & fragment, 107 const G4double K) 108 { 109 theA=GetA(); 110 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 118 119 if (OPTxs==0) {std::ostringstream errOs; 120 errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (Alpha's)!!" 121 <<G4endl; 122 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 123 return 0.;} 124 125 if( OPTxs==1 || OPTxs==2) return G4AlphaEvaporationProbability::GetOpt12( K); 126 else if (OPTxs==3 || OPTxs==4) return G4AlphaEvaporationProbability::GetOpt34( K); 127 else{ 128 std::ostringstream errOs; 129 errOs << "BAD Alpha CROSS SECTION OPTION AT EVAPORATION!!" <<G4endl; 130 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 131 return 0.; 132 } 133 } 134 135 136 // 137 //********************* OPT=1,2 : Chatterjee's cross section ************************ 138 //(fitting to cross section from Bechetti & Greenles OM potential) 139 140 G4double G4AlphaEvaporationProbability::GetOpt12(const G4double K) 141 { 142 // c ** alpha from huizenga and igo 143 G4double Kc=K; 144 145 // JMQ xsec is set constat above limit of validity 146 if (K>50) Kc=50; 147 148 G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs; 149 150 G4double p0 = 10.95; 151 G4double p1 = -85.2; 152 G4double p2 = 1146.; 153 G4double landa0 = 0.0643; 154 G4double landa1 = -13.96; 155 G4double mu0 = 781.2; 156 G4double mu1 = 0.29; 157 G4double nu0 = -304.7; 158 G4double nu1 = -470.0; 159 G4double nu2 = -8.580; 160 G4double delta=1.2; 161 162 163 Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta); 164 p = p0 + p1/Ec + p2/(Ec*Ec); 165 landa = landa0*ResidualA + landa1; 166 mu = mu0*std::pow(ResidualA,mu1); 167 nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec)); 168 q = landa - nu/(Ec*Ec) - 2*p*Ec; 169 r = mu + 2*nu/Ec + p*(Ec*Ec); 170 171 ji=std::max(Kc,Ec); 172 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;} 173 else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;} 174 175 if (xs <0.0) {xs=0.0;} 176 177 return xs; 178 179 } 180 181 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)************* 182 G4double G4AlphaEvaporationProbability::GetOpt34(const G4double K) 183 // c ** alpha from huizenga and igo 184 { 185 186 G4double landa, mu, nu, p , signor(1.),sig; 187 G4double ec,ecsq,xnulam,etest(0.),a; 188 G4double b,ecut,cut,ecut2,geom,elab; 189 190 191 G4double flow = 1.e-18; 192 G4double spill= 1.e+18; 193 194 G4double p0 = 10.95; 195 G4double p1 = -85.2; 196 G4double p2 = 1146.; 197 G4double landa0 = 0.0643; 198 G4double landa1 = -13.96; 199 G4double mu0 = 781.2; 200 G4double mu1 = 0.29; 201 G4double nu0 = -304.7; 202 G4double nu1 = -470.0; 203 G4double nu2 = -8.580; 204 205 G4double ra=1.20; 206 207 //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 ecsq = ec * ec; 211 p = p0 + p1/ec + p2/ecsq; 212 landa = landa0*ResidualA + landa1; 213 a = std::pow(ResidualA,mu1); 214 mu = mu0 * a; 215 nu = a* (nu0+nu1*ec+nu2*ecsq); 216 xnulam = nu / landa; 217 if (xnulam > spill) xnulam=0.; 218 if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam); 219 220 a = -2.*p*ec + landa - nu/ecsq; 221 b = p*ecsq + mu + 2.*nu/ec; 222 ecut = 0.; 223 cut = a*a - 4.*p*b; 224 if (cut > 0.) ecut = std::sqrt(cut); 225 ecut = (ecut-a) / (p+p); 226 ecut2 = ecut; 227 if (cut < 0.) ecut2 = ecut - 2.; 228 elab = K * FragmentA / ResidualA; 229 sig = 0.; 230 231 if (elab <= ec) { //start for E<Ec 232 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 233 } //end for E<Ec 234 else { //start for E>Ec 235 sig = (landa*elab+mu+nu/elab) * signor; 236 geom = 0.; 237 if (xnulam < flow || elab < etest) return sig; 238 geom = std::sqrt(theA*K); 239 geom = 1.23*ResidualAthrd + ra + 4.573/geom; 240 geom = 31.416 * geom * geom; 241 sig = std::max(geom,sig); 242 } //end for E>Ec 243 return sig; 244 245 } 246 247 // ************************** end of cross sections ******************************* -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationChannel.cc
r819 r962 26 26 // 27 27 // $Id: G4DeuteronEvaporationChannel.cc,v 1.4 2006/06/29 20:10:21 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationProbability.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4DeuteronEvaporationProbability.cc,v 1.4 2006/06/29 20:10:23 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara (Nov 1999) 32 // 29 // by V. Lara (Oct 1998) 30 // 31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 32 // cross section option 33 33 34 34 35 35 #include "G4DeuteronEvaporationProbability.hh" 36 36 37 37 38 G4DeuteronEvaporationProbability::G4DeuteronEvaporationProbability() : 38 G4EvaporationProbability(2,1,6) // A,Z,Gamma 39 { 40 std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1; 41 std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1; 42 ExcitEnergies.reserve(NumExcitedStatesEnergy); 43 ExcitSpins.reserve(NumExcitedStatesSpin); 44 ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0); 45 ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0); 46 47 48 ExcitEnergies[15] = 6.18*MeV; 49 ExcitEnergies[17] = 2.15*MeV; 50 ExcitEnergies[18] = 5.02*MeV; 51 ExcitEnergies[19] = 2.65*MeV; 52 ExcitEnergies[20] = 4.80*MeV; 53 ExcitEnergies[22] = 3.85*MeV; 54 ExcitEnergies[23] = 6.96*MeV; 55 ExcitEnergies[25] = 4.92*MeV; 56 ExcitEnergies[26] = 7.22*MeV; 57 ExcitEnergies[27] = 0.40*MeV; 58 ExcitEnergies[28] = 6.83*MeV; 59 ExcitEnergies[29] = 7.12*MeV; 60 ExcitEnergies[30] = 3.84*MeV; 61 ExcitEnergies[31] = 3.92*MeV; 62 63 ExcitSpins[15] = 1; 64 ExcitSpins[17] = 3; 65 ExcitSpins[18] = 4; 66 ExcitSpins[19] = 4; 67 ExcitSpins[20] = 4; 68 ExcitSpins[22] = 6; 69 ExcitSpins[23] = 6; 70 ExcitSpins[25] = 1; 71 ExcitSpins[26] = 10; 72 ExcitSpins[27] = 3; 73 ExcitSpins[28] = 10; 74 ExcitSpins[29] = 3; 75 ExcitSpins[30] = 6; 76 ExcitSpins[31] = 5; 77 78 SetExcitationEnergiesPtr(&ExcitEnergies); 79 SetExcitationSpinsPtr(&ExcitSpins); 80 81 } 39 G4EvaporationProbability(2,1,3,&theCoulombBarrier) // A,Z,Gamma (fixed JMQ) 40 { } 82 41 G4DeuteronEvaporationProbability::G4DeuteronEvaporationProbability(const G4DeuteronEvaporationProbability &) : G4EvaporationProbability() 83 42 { 84 43 throw G4HadronicException(__FILE__, __LINE__, "G4DeuteronEvaporationProbability::copy_constructor meant to not be accessable"); 85 44 } 86 87 88 45 89 46 … … 101 58 } 102 59 60 103 61 G4bool G4DeuteronEvaporationProbability::operator!=(const G4DeuteronEvaporationProbability &) const 104 62 { … … 106 64 } 107 65 108 G4double G4DeuteronEvaporationProbability::CCoeficient(const G4double aZ) const 66 67 G4double G4DeuteronEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 68 { 69 return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ())); 70 } 71 72 73 G4double G4DeuteronEvaporationProbability::CalcBetaParam(const G4Fragment & ) 74 { 75 return 0.0; 76 } 77 78 79 G4double G4DeuteronEvaporationProbability::CCoeficient(const G4double aZ) 109 80 { 110 81 // Data comes from … … 127 98 128 99 } 100 101 102 /////////////////////////////////////////////////////////////////////////////////// 103 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 104 //OPT=0 Dostrovski's parameterization 105 //OPT=1,2 Chatterjee's paramaterization 106 //OPT=3,4 Kalbach's parameterization 107 // 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; 114 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 ;} 172 173 if (xs <0.0) {xs=0.0;} 174 175 return xs; 176 177 } 178 179 180 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)************* 181 G4double G4DeuteronEvaporationProbability::GetOpt34(const G4double K) 182 // ** d from o.m. of perey and perey 183 { 184 185 G4double landa, mu, nu, p , signor(1.),sig; 186 G4double ec,ecsq,xnulam,etest(0.),a; 187 G4double b,ecut,cut,ecut2,geom,elab; 188 189 190 G4double flow = 1.e-18; 191 G4double spill= 1.e+18; 192 193 194 G4double p0 = 0.798; 195 G4double p1 = 420.3; 196 G4double p2 = -1651.; 197 G4double landa0 = 0.00619; 198 G4double landa1 = -7.54; 199 G4double mu0 = 583.5; 200 G4double mu1 = 0.337; 201 G4double nu0 = 421.8; 202 G4double nu1 = -474.5; 203 G4double nu2 = -3.592; 204 205 G4double ra=0.80; 206 207 //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 ecsq = ec * ec; 211 p = p0 + p1/ec + p2/ecsq; 212 landa = landa0*ResidualA + landa1; 213 a = std::pow(ResidualA,mu1); 214 mu = mu0 * a; 215 nu = a* (nu0+nu1*ec+nu2*ecsq); 216 xnulam = nu / landa; 217 if (xnulam > spill) xnulam=0.; 218 if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam); 219 220 a = -2.*p*ec + landa - nu/ecsq; 221 b = p*ecsq + mu + 2.*nu/ec; 222 ecut = 0.; 223 cut = a*a - 4.*p*b; 224 if (cut > 0.) ecut = std::sqrt(cut); 225 ecut = (ecut-a) / (p+p); 226 ecut2 = ecut; 227 if (cut < 0.) ecut2 = ecut - 2.; 228 elab = K * FragmentA / ResidualA; 229 sig = 0.; 230 231 if (elab <= ec) { //start for E<Ec 232 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 233 } //end for E<Ec 234 else { //start for E>Ec 235 sig = (landa*elab+mu+nu/elab) * signor; 236 geom = 0.; 237 if (xnulam < flow || elab < etest) return sig; 238 geom = std::sqrt(theA*K); 239 geom = 1.23*ResidualAthrd + ra + 4.573/geom; 240 geom = 31.416 * geom * geom; 241 sig = std::max(geom,sig); 242 } //end for E>Ec 243 return sig; 244 245 } 246 247 // ************************** end of cross sections ******************************* -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4Evaporation.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4Evaporation.cc,v 1. 7 2007/02/14 13:37:49ahoward Exp $28 // GEANT4 tag $Name: $27 // $Id: G4Evaporation.cc,v 1.12 2008/12/09 17:57:36 ahoward Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations … … 33 33 // Alex Howard - added protection for negative probabilities in the sum, 14/2/07 34 34 // 35 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 36 // cross section option 37 // JMQ (06 September 2008) Also external choices have been added for 38 // superimposed Coulomb barrier (if useSICBis set true, by default is false) 39 35 40 #include "G4Evaporation.hh" 36 41 #include "G4EvaporationFactory.hh" … … 89 94 } 90 95 91 92 96 G4FragmentVector * G4Evaporation::BreakItUp(const G4Fragment &theNucleus) 93 97 { … … 105 109 // Number of channels 106 110 G4int TotNumberOfChannels = theChannels->size(); 107 108 111 109 112 // Starts loop over evaporated particles 110 113 for (;;) 111 { 114 115 { 112 116 // loop over evaporation channels 113 117 std::vector<G4VEvaporationChannel*>::iterator i; 114 118 for (i=theChannels->begin(); i != theChannels->end(); i++) 115 119 { 120 // for inverse cross section choice 121 (*i)->SetOPTxs(OPTxs); 122 // for superimposed Coulomb Barrier for inverse cross sections 123 (*i)->UseSICB(useSICB); 124 116 125 (*i)->Initialize(theResidualNucleus); 117 126 } … … 172 181 if( j >= TotNumberOfChannels ) 173 182 { 183 G4cerr << " Residual A: " << theResidualNucleus.GetA() << " Residual Z: " << theResidualNucleus.GetZ() << " Excitation Energy: " << theResidualNucleus.GetExcitationEnergy() << G4endl; 184 G4cerr << " j has not chosen a channel, j = " << j << " TotNumberOfChannels " << TotNumberOfChannels << " Total Probability: " << TotalProbability << G4endl; 185 for (j=0; j < TotNumberOfChannels; j++) 186 { 187 G4cerr << " j: " << j << " EmissionProbChannel: " << EmissionProbChannel[j] << " and shoot: " << shoot << " (<ProbChannel?) " << G4endl; 188 } 174 189 throw G4HadronicException(__FILE__, __LINE__, "G4Evaporation::BreakItUp: Can't define emission probability of the channels!" ); 175 190 } … … 216 231 theResidualNucleus.SetCreatorModel(G4String("ResidualNucleus")); 217 232 #endif 233 218 234 } 219 235 } -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationChannel.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4EvaporationChannel.cc,v 1.5 2006/06/29 20:10:27 gunter Exp $ 28 // GEANT4 tag $Name: $ 27 //J.M. Quesada (August2008). Based on: 29 28 // 30 29 // Hadronic Process: Nuclear De-excitations 31 30 // by V. Lara (Oct 1998) 32 31 // 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 33 37 34 38 #include "G4EvaporationChannel.hh" … … 36 40 37 41 38 G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ, 42 43 G4EvaporationChannel::G4EvaporationChannel(const G4int anA, const G4int aZ, const G4String & aName, 39 44 G4VEmissionProbability * aEmissionStrategy, 40 G4VCoulombBarrier * aCoulombBarrier): 41 A(theA), 42 Z(theZ), 45 G4VCoulombBarrier * aCoulombBarrier): 46 G4VEvaporationChannel(aName), 47 theA(anA), 48 theZ(aZ), 43 49 theEvaporationProbabilityPtr(aEmissionStrategy), 44 50 theCoulombBarrierPtr(aCoulombBarrier), … … 47 53 { 48 54 theLevelDensityPtr = new G4EvaporationLevelDensityParameter; 49 MyOwnLevelDensity = true; 50 } 51 52 G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String & aName, 53 G4VEmissionProbability * aEmissionStrategy, 54 G4VCoulombBarrier * aCoulombBarrier): 55 G4VEvaporationChannel(aName), 56 A(theA), 57 Z(theZ), 58 theEvaporationProbabilityPtr(aEmissionStrategy), 59 theCoulombBarrierPtr(aCoulombBarrier), 60 EmissionProbability(0.0), 61 MaximalKineticEnergy(-1000.0) 62 { 63 theLevelDensityPtr = new G4EvaporationLevelDensityParameter; 64 MyOwnLevelDensity = true; 65 } 66 67 G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String * aName, 68 G4VEmissionProbability * aEmissionStrategy, 69 G4VCoulombBarrier * aCoulombBarrier): 70 G4VEvaporationChannel(aName), 71 A(theA), 72 Z(theZ), 73 theEvaporationProbabilityPtr(aEmissionStrategy), 74 theCoulombBarrierPtr(aCoulombBarrier), 75 EmissionProbability(0.0), 76 MaximalKineticEnergy(-1000.0) 77 { 78 theLevelDensityPtr = new G4EvaporationLevelDensityParameter; 79 MyOwnLevelDensity = true; 80 } 81 55 // MyOwnLevelDensity = true; 56 } 82 57 83 58 G4EvaporationChannel::~G4EvaporationChannel() 84 59 { 85 if (MyOwnLevelDensity) delete theLevelDensityPtr; 86 } 87 88 89 60 delete theLevelDensityPtr; 61 } 90 62 91 63 G4EvaporationChannel::G4EvaporationChannel(const G4EvaporationChannel & ) : G4VEvaporationChannel() … … 112 84 } 113 85 114 86 //void G4EvaporationChannel::SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity) 87 // { 88 // if (MyOwnLevelDensity) delete theLevelDensityPtr; 89 // theLevelDensityPtr = aLevelDensity; 90 // MyOwnLevelDensity = false; 91 // } 115 92 116 93 void G4EvaporationChannel::Initialize(const G4Fragment & fragment) 117 94 { 118 119 G4int anA = static_cast<G4int>(fragment.GetA()); 120 G4int aZ = static_cast<G4int>(fragment.GetZ()); 121 AResidual = anA - A; 122 ZResidual = aZ - Z; 123 124 // Effective excitation energy 125 G4double ExEnergy = fragment.GetExcitationEnergy() - 126 G4PairingCorrection::GetInstance()->GetPairingCorrection(anA,aZ); 127 128 // We only take into account channels which are physically allowed 129 if (AResidual <= 0 || ZResidual <= 0 || AResidual < ZResidual || 130 (AResidual == ZResidual && AResidual > 1) || ExEnergy <= 0.0) { 131 // LevelDensityParameter = 0.0; 132 CoulombBarrier = 0.0; 133 // BindingEnergy = 0.0; 134 MaximalKineticEnergy = -1000.0*MeV; 135 EmissionProbability = 0.0; 136 } else { 137 // // Get Level Density 138 // LevelDensityParameter = theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy); 139 140 // Coulomb Barrier calculation 141 CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(AResidual,ZResidual,ExEnergy); 142 143 // // Binding Enegy (for separate fragment from nucleus) 144 // BindingEnergy = CalcBindingEnergy(anA,aZ); 145 146 // Maximal Kinetic Energy 147 MaximalKineticEnergy = CalcMaximalKineticEnergy(G4ParticleTable::GetParticleTable()-> 148 GetIonTable()->GetNucleusMass(aZ,anA)+ExEnergy); 149 150 // Emission probability 151 if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0; 152 else { 153 // Total emission probability for this channel 154 EmissionProbability = theEvaporationProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy); 155 156 } 95 //for inverse cross section choice 96 theEvaporationProbabilityPtr->SetOPTxs(OPTxs); 97 // for superimposed Coulomb Barrier for inverse cross sections 98 theEvaporationProbabilityPtr->UseSICB(useSICB); 99 100 101 G4int FragmentA = static_cast<G4int>(fragment.GetA()); 102 G4int FragmentZ = static_cast<G4int>(fragment.GetZ()); 103 ResidualA = FragmentA - theA; 104 ResidualZ = FragmentZ - theZ; 105 106 //Effective excitation energy 107 G4double ExEnergy = fragment.GetExcitationEnergy() - 108 G4PairingCorrection::GetInstance()->GetPairingCorrection(FragmentA,FragmentZ); 109 110 111 // Only channels which are physically allowed are taken into account 112 if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ || 113 (ResidualA == ResidualZ && ResidualA > 1) || ExEnergy <= 0.0) { 114 CoulombBarrier=0.0; 115 MaximalKineticEnergy = -1000.0*MeV; 116 EmissionProbability = 0.0; 117 } else { 118 CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy); 119 // Maximal Kinetic Energy 120 MaximalKineticEnergy = CalcMaximalKineticEnergy 121 (G4ParticleTable::GetParticleTable()-> 122 GetIonTable()->GetNucleusMass(FragmentZ,FragmentA)+ExEnergy); 123 124 // Emission probability 125 // Protection for the case Tmax<V. If not set in this way we could end up in an 126 // infinite loop in the method GetKineticEnergy if OPTxs!=0 && useSICB=true. 127 // Of course for OPTxs=0 we have the Coulomb barrier 128 129 G4double limit; 130 if (OPTxs==0 || (OPTxs!=0 && useSICB)) 131 limit= CoulombBarrier; 132 else limit=0.; 133 134 // The threshold for charged particle emission must be set to 0 if Coulomb 135 //cutoff is included in the cross sections 136 //if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0; 137 if (MaximalKineticEnergy <= limit) EmissionProbability = 0.0; 138 else { 139 // Total emission probability for this channel 140 EmissionProbability = theEvaporationProbabilityPtr-> 141 EmissionProbability(fragment, MaximalKineticEnergy); 157 142 } 158 159 return;160 } 161 143 } 144 145 return; 146 } 162 147 163 148 G4FragmentVector * G4EvaporationChannel::BreakUp(const G4Fragment & theNucleus) 164 149 { 165 166 G4double EvaporatedKineticEnergy = CalcKineticEnergy(); 167 G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A); 168 G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass; 169 170 G4ThreeVector momentum(IsotropicVector(std::sqrt(EvaporatedKineticEnergy* 171 (EvaporatedKineticEnergy+2.0*EvaporatedMass)))); 172 173 momentum.rotateUz(theNucleus.GetMomentum().vect().unit()); 174 175 G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy); 176 EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector()); 177 178 G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum); 150 G4double EvaporatedKineticEnergy=GetKineticEnergy(theNucleus); 151 152 G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()-> 153 GetIonMass(theZ,theA); 154 G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass; 155 156 G4ThreeVector momentum(IsotropicVector 157 (std::sqrt(EvaporatedKineticEnergy* 158 (EvaporatedKineticEnergy+2.0*EvaporatedMass)))); 159 160 momentum.rotateUz(theNucleus.GetMomentum().vect().unit()); 161 162 G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy); 163 EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector()); 164 165 G4Fragment * EvaporatedFragment = new G4Fragment(theA,theZ,EvaporatedMomentum); 179 166 #ifdef PRECOMPOUND_TEST 180 167 EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation")); 181 168 #endif 182 // ** And now the residual nucleus ** 183 G4double theExEnergy = theNucleus.GetExcitationEnergy(); 184 G4double theMass = G4ParticleTable::GetParticleTable()->GetIonTable()-> 185 GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()), 186 static_cast<G4int>(theNucleus.GetA())); 187 G4double ResidualEnergy = theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass; 188 189 G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy); 190 ResidualMomentum.boost(theNucleus.GetMomentum().boostVector()); 191 192 G4Fragment * ResidualFragment = new G4Fragment( AResidual, ZResidual, ResidualMomentum ); 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 193 182 #ifdef PRECOMPOUND_TEST 194 183 ResidualFragment->SetCreatorModel(G4String("ResidualNucleus")); 195 184 #endif 196 197 185 G4FragmentVector * theResult = new G4FragmentVector; 186 198 187 #ifdef debug 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 188 G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e(); 189 G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect(); 190 if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 1.0*keV) { 191 G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl; 192 G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV Final :" 193 <<Efinal/MeV << " MeV Delta: " << (Efinal-theNucleus.GetMomentum().e())/MeV 194 << " MeV" << G4endl; 195 } 196 if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 1.0*keV || 197 std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 1.0*keV || 198 std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 1.0*keV ) { 199 G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl; 200 G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV Final :" 201 <<Pfinal/MeV << " MeV Delta: " << Pfinal-theNucleus.GetMomentum().vect() 202 << " MeV" << G4endl; 203 204 } 216 205 #endif 217 218 219 206 theResult->push_back(EvaporatedFragment); 207 theResult->push_back(ResidualFragment); 208 return theResult; 220 209 } 221 210 222 223 224 // G4double G4EvaporationChannel::CalcBindingEnergy(const G4int anA, const G4int aZ) 225 // // Calculate Binding Energy for separate fragment from nucleus 226 // { 227 // // Mass Excess for residual nucleus 228 // G4double ResNucMassExcess = G4NucleiProperties::GetNuclearMass(AResidual,ZResidual); 229 // // Mass Excess for fragment 230 // G4double FragmentMassExcess = G4NucleiProperties::GetNuclearMass(A,Z); 231 // // Mass Excess for Compound Nucleus 232 // G4double NucleusMassExcess = G4NucleiProperties::GetNuclearMass(anA,aZ); 233 // 234 // return ResNucMassExcess + FragmentMassExcess - NucleusMassExcess; 235 // } 236 237 211 ///////////////////////////////////////// 212 // Calculates the maximal kinetic energy that can be carried by fragment. 238 213 G4double G4EvaporationChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE) 239 // Calculate maximal kinetic energy that can be carried by fragment. 240 { 241 G4double ResidualMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass( ZResidual, AResidual ); 242 G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass( Z, A ); 243 244 G4double T = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - ResidualMass*ResidualMass)/ 245 (2.0*NucleusTotalE) - 246 EvaporatedMass - CoulombBarrier; 247 248 return T; 249 } 250 251 252 253 254 G4double G4EvaporationChannel::CalcKineticEnergy(void) 255 // Samples fragment kinetic energy. 214 { 215 G4double ResidualMass = G4ParticleTable::GetParticleTable()-> 216 GetIonTable()->GetNucleusMass( ResidualZ, ResidualA ); 217 G4double EvaporatedMass = G4ParticleTable::GetParticleTable()-> 218 GetIonTable()->GetNucleusMass( theZ, theA ); 219 220 // This is the "true" assimptotic kinetic energy (from energy conservation) 221 G4double Tmax = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - 222 ResidualMass*ResidualMass)/(2.0*NucleusTotalE) - EvaporatedMass; 223 224 //JMQ (13-09-08) bug fixed: in the original version the Tmax is calculated 225 //at the Coulomb barrier 226 //IMPORTANT: meaning of Tmax differs in OPTxs=0 and OPTxs!=0 227 //When OPTxs!=0 Tmax is the TRUE (assimptotic) maximal kinetic energy 228 229 if(OPTxs==0) 230 Tmax=Tmax- CoulombBarrier; 231 232 return Tmax; 233 } 234 235 /////////////////////////////////////////// 236 //JMQ: New method for MC sampling of kinetic energy. Substitutes old CalcKineticEnergy 237 G4double G4EvaporationChannel::GetKineticEnergy(const G4Fragment & aFragment) 238 { 239 240 if (OPTxs==0) { 256 241 // It uses Dostrovsky's approximation for the inverse reaction cross 257 // in the probability for fragment emisson 258 { 259 if (MaximalKineticEnergy < 0.0) 260 throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic energy is less than 0"); 261 262 G4double Rb = 4.0*theLevelDensityPtr->LevelDensityParameter(AResidual+A,ZResidual+Z,MaximalKineticEnergy)* 263 MaximalKineticEnergy; 242 // in the probability for fragment emission 243 // MaximalKineticEnergy energy in the original version (V.Lara) was calculated at 244 //the Coulomb barrier. 245 246 247 if (MaximalKineticEnergy < 0.0) 248 throw G4HadronicException(__FILE__, __LINE__, 249 "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic at the Coulomb barrier is less than 0"); 250 251 G4double Rb = 4.0*theLevelDensityPtr-> 252 LevelDensityParameter(ResidualA+theA,ResidualZ+theZ,MaximalKineticEnergy)* 253 MaximalKineticEnergy; 264 254 G4double RbSqrt = std::sqrt(Rb); 265 255 G4double PEX1 = 0.0; … … 268 258 G4double FRk = 0.0; 269 259 do { 270 271 272 273 274 if (Z == 0) { // for emitted neutron275 G4double Beta = (2.12/std::pow(G4double(AResidual),2./3.) - 0.05)*MeV/276 (0.76 + 2.2/std::pow(G4double(AResidual),1./3.));277 278 279 280 281 282 260 G4double RandNumber = G4UniformRand(); 261 Rk = 1.0 + (1./RbSqrt)*std::log(RandNumber + (1.0-RandNumber)*PEX1); 262 G4double Q1 = 1.0; 263 G4double Q2 = 1.0; 264 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.)); 267 Q1 = 1.0 + Beta/(MaximalKineticEnergy); 268 Q2 = Q1*std::sqrt(Q1); 269 } 270 271 FRk = (3.0*std::sqrt(3.0)/2.0)/Q2 * Rk * (Q1 - Rk*Rk); 272 283 273 } while (FRk < G4UniformRand()); 284 274 285 275 G4double result = MaximalKineticEnergy * (1.0-Rk*Rk) + CoulombBarrier; 286 287 276 return result; 288 } 289 277 } else if (OPTxs==1 || OPTxs==2 || OPTxs==3 || OPTxs==4) { 278 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.; 284 285 G4double Tmax=MaximalKineticEnergy; 286 G4double T(0.0); 287 G4double NormalizedProbability(1.0); 288 289 // A pointer is created in order to access the distribution function. 290 G4EvaporationProbability * G4EPtemp; 291 292 if (theA==1 && theZ==0) G4EPtemp=new G4NeutronEvaporationProbability(); 293 else if (theA==1 && theZ==1) G4EPtemp=new G4ProtonEvaporationProbability(); 294 else if (theA==2 && theZ==1 ) G4EPtemp=new G4DeuteronEvaporationProbability(); 295 else if (theA==3 && theZ==1 ) G4EPtemp=new G4TritonEvaporationProbability(); 296 else if (theA==3 && theZ==2 ) G4EPtemp=new G4He3EvaporationProbability(); 297 else if (theA==4 && theZ==2) G4EPtemp=new G4AlphaEvaporationProbability(); 298 else { 299 std::ostringstream errOs; 300 errOs << "ejected particle out of range in G4EvaporationChannel" << G4endl; 301 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 302 } 303 304 //for cross section selection and superimposed Coulom Barrier for xs 305 G4EPtemp->SetOPTxs(OPTxs); 306 G4EPtemp->UseSICB(useSICB); 307 308 do 309 { 310 T=V+G4UniformRand()*(Tmax-V); 311 NormalizedProbability=G4EPtemp->ProbabilityDistributionFunction(aFragment,T)/ 312 (this->GetEmissionProbability()); 313 314 } 315 while (G4UniformRand() > NormalizedProbability); 316 delete G4EPtemp; 317 return T; 318 } else{ 319 std::ostringstream errOs; 320 errOs << "Bad option for energy sampling in evaporation" <<G4endl; 321 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 322 } 323 } 290 324 291 325 G4ThreeVector G4EvaporationChannel::IsotropicVector(const G4double Magnitude) … … 300 334 Magnitude*CosTheta); 301 335 return Vector; 302 } 303 304 305 336 } 337 338 339 340 341 342 343 344 345 346 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationFactory.cc
r819 r962 26 26 // 27 27 // $Id: G4EvaporationFactory.cc,v 1.4 2006/06/29 20:10:29 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4EvaporationProbability.cc,v 1.5 2006/06/29 20:10:31 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 29 // by V. Lara (Oct 1998) 32 30 // 33 31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 32 // cross section option 33 // JMQ (06 September 2008) Also external choices have been added for 34 // superimposed Coulomb barrier (if useSICB is set true, by default is false) 35 // 36 // JMQ (14 february 2009) bug fixed in emission width: hbarc instead of hbar_Planck in the denominator 37 // 38 #include <iostream> 39 using namespace std; 34 40 35 41 #include "G4EvaporationProbability.hh" … … 63 69 return true; 64 70 } 65 71 66 72 G4double G4EvaporationProbability::EmissionProbability(const G4Fragment & fragment, const G4double anEnergy) 67 73 { 68 74 G4double EmissionProbability = 0.0; 69 75 G4double MaximalKineticEnergy = anEnergy; 76 70 77 if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) { 71 EmissionProbability = CalcProbability(fragment,MaximalKineticEnergy); 72 // // Next there is a loop over excited states for this channel summing probabilities 73 // G4double SavedGamma = Gamma; 74 // for (G4int i = 0; i < ExcitationEnergies->length(); i++) { 75 // if (ExcitationSpins->operator()(i) < 0.1) continue; 76 // Gamma = ExcitationSpins->operator()(i)*A; // A is the channel mass number 77 // // substract excitation energies 78 // MaximalKineticEnergy -= ExcitationEnergies->operator()(i); 79 // // update probability 80 // EmissionProbability += CalcProbability(fragment,MaximalKineticEnergy); 81 // EmissionProbability += tmp; 82 // } 83 // // restore Gamma and MaximalKineticEnergy 84 // MaximalKineticEnergy = SavedMaximalKineticEnergy; 85 // Gamma = SavedGamma; 86 // } 78 EmissionProbability = CalculateProbability(fragment, MaximalKineticEnergy); 79 87 80 } 88 81 return EmissionProbability; 89 82 } 90 83 91 G4double G4EvaporationProbability::CalcProbability(const G4Fragment & fragment, 92 const G4double MaximalKineticEnergy) 93 // Calculate integrated probability (width) for rvaporation channel 94 { 95 G4double ResidualA = static_cast<G4double>(fragment.GetA() - theA); 96 G4double ResidualZ = static_cast<G4double>(fragment.GetZ() - theZ); 84 //////////////////////////////////// 85 86 // Computes the integrated probability of evaporation channel 87 G4double G4EvaporationProbability::CalculateProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy) 88 { 89 G4double ResidualA = fragment.GetA() - theA; 90 G4double ResidualZ = fragment.GetZ() - theZ; 97 91 G4double U = fragment.GetExcitationEnergy(); 92 93 if (OPTxs==0) { 94 98 95 99 96 G4double NuclearMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(theZ,theA); … … 107 104 (U-delta0)); 108 105 109 // compute the integrated probability of evaporation channel 106 110 107 G4double RN = 1.5*fermi; 111 108 … … 133 130 134 131 return Width; 135 } 136 137 138 139 132 133 } else if (OPTxs==1 || OPTxs==2 ||OPTxs==3 || OPTxs==4) { 134 135 G4double limit; 136 if (useSICB) limit=theCoulombBarrierptr->GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U); 137 else limit=0.; 138 139 if (MaximalKineticEnergy <= limit) return 0.0; 140 141 142 // if Coulomb barrier cutoff is superimposed for all cross sections the limit is the Coulomb Barrier 143 G4double LowerLimit= limit; 144 145 // MaximalKineticEnergy: asimptotic value (already accounted for in G4EvaporationChannel) 146 147 G4double UpperLimit = MaximalKineticEnergy; 148 149 150 G4double Width = IntegrateEmissionProbability(fragment,LowerLimit,UpperLimit); 151 152 return Width; 153 } else{ 154 std::ostringstream errOs; 155 errOs << "Bad option for cross sections at evaporation" <<G4endl; 156 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 157 } 158 159 } 160 161 ///////////////////////////////////////////////////////////////////// 162 163 G4double G4EvaporationProbability:: 164 IntegrateEmissionProbability(const G4Fragment & fragment, const G4double & Low, const G4double & Up ) 165 { 166 167 static const G4int N = 10; 168 // 10-Points Gauss-Legendre abcisas and weights 169 static const G4double w[N] = { 170 0.0666713443086881, 171 0.149451349150581, 172 0.219086362515982, 173 0.269266719309996, 174 0.295524224714753, 175 0.295524224714753, 176 0.269266719309996, 177 0.219086362515982, 178 0.149451349150581, 179 0.0666713443086881 180 }; 181 static const G4double x[N] = { 182 -0.973906528517172, 183 -0.865063366688985, 184 -0.679409568299024, 185 -0.433395394129247, 186 -0.148874338981631, 187 0.148874338981631, 188 0.433395394129247, 189 0.679409568299024, 190 0.865063366688985, 191 0.973906528517172 192 }; 193 194 G4double Total = 0.0; 195 196 197 for (G4int i = 0; i < N; i++) 198 { 199 200 G4double KineticE = ((Up-Low)*x[i]+(Up+Low))/2.0; 201 202 Total += w[i]*ProbabilityDistributionFunction(fragment, KineticE); 203 204 } 205 Total *= (Up-Low)/2.0; 206 return Total; 207 } 208 209 210 ///////////////////////////////////////////////////////// 211 //New method (OPT=1,2,3,4) 212 213 G4double G4EvaporationProbability::ProbabilityDistributionFunction( const G4Fragment & fragment, const G4double K) 214 { 215 216 217 218 219 G4double ResidualA = fragment.GetA() - theA; 220 G4double ResidualZ = fragment.GetZ() - theZ; 221 G4double U = fragment.GetExcitationEnergy(); 222 223 // if(K <= theCoulombBarrierptr->GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U)) return 0.0; 224 225 G4double delta1 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ)); 226 227 228 G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ())); 229 230 231 G4double ParticleMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(theZ,theA); 232 233 G4double theSeparationEnergy= G4NucleiProperties::GetMassExcess(static_cast<G4int>(theA),static_cast<G4int>(theZ)) + 234 G4NucleiProperties::GetMassExcess(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ)) - 235 G4NucleiProperties::GetMassExcess(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ())); 236 237 G4double a0 = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()),U - delta0); 238 239 G4double a1 = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA),static_cast<G4int>(ResidualZ),U - theSeparationEnergy - delta1); 240 241 242 G4double E0=U-delta0; 243 244 G4double E1=U-theSeparationEnergy-delta1-K; 245 246 if (E1<0.) return 0.; 247 248 //JMQ 14/02/09 BUG fixed: hbarc should be in the denominator instead of hbar_Planck 249 //Without 1/hbar_Panck remains as a width 250 // G4double Prob=Gamma*ParticleMass/((pi*hbar_Planck)*(pi*hbar_Planck)* 251 //std::exp(2*std::sqrt(a0*E0)))*K*CrossSection(fragment,K)*std::exp(2*std::sqrt(a1*E1))*millibarn; 252 253 G4double Prob=Gamma*ParticleMass/((pi*hbarc)*(pi*hbarc)*std::exp(2*std::sqrt(a0*E0))) 254 *K*CrossSection(fragment,K)*std::exp(2*std::sqrt(a1*E1))*millibarn; 255 256 return Prob; 257 } 258 259 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationChannel.cc
r819 r962 26 26 // 27 27 // $Id: G4He3EvaporationChannel.cc,v 1.4 2006/06/29 20:10:33 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationProbability.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4He3EvaporationProbability.cc,v 1.4 2006/06/29 20:10:35 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara (Nov 1999) 32 // 33 29 // by V. Lara (Oct 1998) 30 // 31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 32 // cross section option 34 33 35 34 #include "G4He3EvaporationProbability.hh" 36 35 37 36 G4He3EvaporationProbability::G4He3EvaporationProbability() : 38 G4EvaporationProbability(3,2,6) // A,Z,Gamma 39 { 40 std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1; 41 std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1; 42 ExcitEnergies.reserve(NumExcitedStatesEnergy); 43 ExcitSpins.reserve(NumExcitedStatesSpin); 44 ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0); 45 ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0); 46 47 48 49 ExcitEnergies[18] = 7.29*MeV; 50 ExcitEnergies[20] = 6.48*MeV; 51 ExcitEnergies[25] = 5.69*MeV; 52 ExcitEnergies[26] = 8.31*MeV; 53 ExcitEnergies[31] = 5.10*MeV; 54 55 ExcitSpins[18] = 6; 56 ExcitSpins[20] = 8; 57 ExcitSpins[25] = 3; 58 ExcitSpins[26] = 2; 59 ExcitSpins[31] = 7; 60 61 SetExcitationEnergiesPtr(&ExcitEnergies); 62 SetExcitationSpinsPtr(&ExcitSpins); 37 G4EvaporationProbability(3,2,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier 38 { 39 63 40 } 64 41 … … 89 66 } 90 67 91 G4double G4He3EvaporationProbability::CCoeficient(const G4double aZ) const 68 G4double G4He3EvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 69 { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));} 70 71 G4double G4He3EvaporationProbability::CalcBetaParam(const G4Fragment & ) 72 { return 0.0; } 73 74 75 G4double G4He3EvaporationProbability::CCoeficient(const G4double aZ) 92 76 { 93 77 // Data comes from … … 113 97 return C*(4.0/3.0); 114 98 } 99 100 /////////////////////////////////////////////////////////////////////////////////// 101 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 102 //OPT=0 Dostrovski's parameterization 103 //OPT=1,2 Chatterjee's paramaterization 104 //OPT=3,4 Kalbach's parameterization 105 // 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 ************************ 134 //(fitting to cross section from Bechetti & Greenles OM potential) 135 136 G4double G4He3EvaporationProbability::GetOpt12(const G4double K) 137 { 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;} 171 172 return xs; 173 174 } 175 176 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)************* 177 G4double G4He3EvaporationProbability::GetOpt34(const G4double K) 178 //c ** 3he from o.m. of gibson et al 179 { 180 181 G4double landa, mu, nu, p , signor(1.),sig; 182 G4double ec,ecsq,xnulam,etest(0.),a; 183 G4double b,ecut,cut,ecut2,geom,elab; 184 185 186 G4double flow = 1.e-18; 187 G4double spill= 1.e+18; 188 189 190 G4double p0 = -2.88; 191 G4double p1 = 205.6; 192 G4double p2 = -1487.; 193 G4double landa0 = 0.00459; 194 G4double landa1 = -8.93; 195 G4double mu0 = 611.2; 196 G4double mu1 = 0.35; 197 G4double nu0 = 473.8; 198 G4double nu1 = -468.2; 199 G4double nu2 = -2.225; 200 201 G4double ra=0.80; 202 203 //JMQ 13/02/09 increase of reduced radius to lower the barrier 204 // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); 205 ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra); 206 ecsq = ec * ec; 207 p = p0 + p1/ec + p2/ecsq; 208 landa = landa0*ResidualA + landa1; 209 a = std::pow(ResidualA,mu1); 210 mu = mu0 * a; 211 nu = a* (nu0+nu1*ec+nu2*ecsq); 212 xnulam = nu / landa; 213 if (xnulam > spill) xnulam=0.; 214 if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam); 215 216 a = -2.*p*ec + landa - nu/ecsq; 217 b = p*ecsq + mu + 2.*nu/ec; 218 ecut = 0.; 219 cut = a*a - 4.*p*b; 220 if (cut > 0.) ecut = std::sqrt(cut); 221 ecut = (ecut-a) / (p+p); 222 ecut2 = ecut; 223 if (cut < 0.) ecut2 = ecut - 2.; 224 elab = K * FragmentA / ResidualA; 225 sig = 0.; 226 227 if (elab <= ec) { //start for E<Ec 228 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 229 } //end for E<Ec 230 else { //start for E>Ec 231 sig = (landa*elab+mu+nu/elab) * signor; 232 geom = 0.; 233 if (xnulam < flow || elab < etest) return sig; 234 geom = std::sqrt(theA*K); 235 geom = 1.23*ResidualAthrd + ra + 4.573/geom; 236 geom = 31.416 * geom * geom; 237 sig = std::max(geom,sig); 238 } //end for E>Ec 239 return sig; 240 241 } 242 243 // ************************** end of cross sections ******************************* 244 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationChannel.cc
r819 r962 26 26 // 27 27 // $Id: G4NeutronEvaporationChannel.cc,v 1.4 2006/06/29 20:10:37 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4NeutronEvaporationProbability.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4NeutronEvaporationProbability.cc,v 1.4 2006/06/29 20:10:39 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara (Nov 1999) 32 // 33 29 // by V. Lara (Oct 1998) 30 // 31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 32 // cross section option 34 33 35 34 #include "G4NeutronEvaporationProbability.hh" 36 35 37 36 G4NeutronEvaporationProbability::G4NeutronEvaporationProbability() : 38 G4EvaporationProbability(1,0,2) // A,Z,Gamma 39 { 40 std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1; 41 std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1; 42 ExcitEnergies.reserve(NumExcitedStatesEnergy); 43 ExcitSpins.reserve(NumExcitedStatesSpin); 44 ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0); 45 ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0); 46 47 48 ExcitEnergies[ 9] = 3.56*MeV; 49 ExcitEnergies[10] = 0.48*MeV; 50 ExcitEnergies[11] = 0.98*MeV; 51 ExcitEnergies[12] = 0.43*MeV; 52 ExcitEnergies[15] = 3.37*MeV; 53 ExcitEnergies[17] = 0.72*MeV; 54 ExcitEnergies[18] = 2.13*MeV; 55 ExcitEnergies[19] = 0.95*MeV; 56 ExcitEnergies[20] = 2.00*MeV; 57 ExcitEnergies[21] = 4.44*MeV; 58 ExcitEnergies[22] = 3.09*MeV; 59 ExcitEnergies[23] = 6.09*MeV; 60 ExcitEnergies[25] = 2.31*MeV; 61 ExcitEnergies[26] = 5.28*MeV; 62 ExcitEnergies[27] = 0.12*MeV; 63 ExcitEnergies[28] = 5.22*MeV; 64 ExcitEnergies[29] = 6.10*MeV; 65 ExcitEnergies[30] = 0.87*MeV; 66 ExcitEnergies[31] = 1.98*MeV; 67 68 69 ExcitSpins[ 9] = 1; 70 ExcitSpins[10] = 2; 71 ExcitSpins[11] = 3; 72 ExcitSpins[12] = 2; 73 ExcitSpins[15] = 5; 74 ExcitSpins[17] = 3; 75 ExcitSpins[18] = 2; 76 ExcitSpins[19] = 5; 77 ExcitSpins[20] = 2; 78 ExcitSpins[21] = 5; 79 ExcitSpins[22] = 2; 80 ExcitSpins[23] = 3; 81 ExcitSpins[25] = 1; 82 ExcitSpins[26] = 8; 83 ExcitSpins[27] = 1; 84 ExcitSpins[28] = 8; 85 ExcitSpins[29] = 8; 86 ExcitSpins[30] = 2; 87 ExcitSpins[31] = 5; 88 89 SetExcitationEnergiesPtr(&ExcitEnergies); 90 SetExcitationSpinsPtr(&ExcitSpins); 37 G4EvaporationProbability(1,0,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier 38 { 39 91 40 } 92 41 … … 117 66 return true; 118 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 72 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 77 78 //////////////////////////////////////////////////////////////////////////////////// 79 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 80 //OPT=0 Dostrovski's parameterization 81 //OPT=1,2 Chatterjee's paramaterization 82 //OPT=3,4 Kalbach's parameterization 83 // 84 G4double G4NeutronEvaporationProbability::CrossSection(const G4Fragment & fragment, const G4double K) 85 { 86 theA=GetA(); 87 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 95 96 if (OPTxs==0) {std::ostringstream errOs; 97 errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (neutrons)!!" <<G4endl; 98 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 99 return 0.;} 100 else if( OPTxs==1 ||OPTxs==2) return GetOpt12( K); 101 else if (OPTxs==3 || OPTxs==4) return GetOpt34( K); 102 else{ 103 std::ostringstream errOs; 104 errOs << "BAD NEUTRON CROSS SECTION OPTION AT EVAPORATION!!" <<G4endl; 105 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 106 return 0.; 107 } 108 } 109 110 111 112 113 114 115 //********************* OPT=1,2 : Chatterjee's cross section ************************ 116 //(fitting to cross section from Bechetti & Greenles OM potential) 117 118 G4double G4NeutronEvaporationProbability::GetOpt12(const G4double K) 119 { 120 121 G4double Kc=K; 122 123 // JMQ xsec is set constat above limit of validity 124 if (K>50) Kc=50; 125 126 G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs; 127 128 landa0 = 18.57; 129 landa1 = -22.93; 130 mu0 = 381.7; 131 mu1 = 24.31; 132 nu0 = 0.172; 133 nu1 = -15.39; 134 nu2 = 804.8; 135 landa = landa0/ResidualAthrd + landa1; 136 mu = mu0*ResidualAthrd + mu1*ResidualAthrd*ResidualAthrd; 137 nu = nu0*ResidualAthrd*ResidualA + nu1*ResidualAthrd*ResidualAthrd + nu2 ; 138 xs=landa*Kc + mu + nu/Kc; 139 if (xs <= 0.0 ){ 140 std::ostringstream errOs; 141 G4cout<<"WARNING: NEGATIVE OPT=1 neutron cross section "<<G4endl; 142 errOs << "RESIDUAL: Ar=" << ResidualA << " Zr=" << ResidualZ <<G4endl; 143 errOs <<" xsec("<<Kc<<" MeV) ="<<xs <<G4endl; 144 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 145 } 146 return xs; 147 } 148 149 150 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)************* 151 G4double G4NeutronEvaporationProbability::GetOpt34(const G4double K) 152 { 153 154 G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2; 155 G4double p, p0, p1, p2; 156 G4double flow,spill,ec,ecsq,xnulam,etest(0.),ra(0.),a,signor(1.),sig; 157 G4double b,ecut,cut,ecut2,geom,elab; 158 159 //safety initialization 160 landa0=0; 161 landa1=0; 162 mu0=0.; 163 mu1=0.; 164 nu0=0.; 165 nu1=0.; 166 nu2=0.; 167 p0=0.; 168 p1=0.; 169 p2=0.; 170 171 172 flow = 1.e-18; 173 spill= 1.0e+18; 174 175 176 177 // PRECO xs for neutrons is choosen 178 179 p0 = -312.; 180 p1= 0.; 181 p2 = 0.; 182 landa0 = 12.10; 183 landa1= -11.27; 184 mu0 = 234.1; 185 mu1 = 38.26; 186 nu0 = 1.55; 187 nu1 = -106.1; 188 nu2 = 1280.8; 189 190 if (ResidualA < 40.) signor=0.7+ResidualA*0.0075; 191 if (ResidualA > 210.) signor = 1. + (ResidualA-210.)/250.; 192 landa = landa0/ResidualAthrd + landa1; 193 mu = mu0*ResidualAthrd + mu1*ResidualAthrd*ResidualAthrd; 194 nu = nu0*ResidualAthrd*ResidualA + nu1*ResidualAthrd*ResidualAthrd + nu2; 195 196 // JMQ very low energy behaviour corrected (problem for A (apprx.)>60) 197 if (nu < 0.)nu=-nu; 198 199 ec = 0.5; 200 ecsq = 0.25; 201 p = p0; 202 xnulam = 1.; 203 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 209 210 a = -2.*p*ec + landa - nu/ecsq; 211 b = p*ecsq + mu + 2.*nu/ec; 212 ecut = 0.; 213 cut = a*a - 4.*p*b; 214 if (cut > 0.) ecut = std::sqrt(cut); 215 ecut = (ecut-a) / (p+p); 216 ecut2 = ecut; 217 if (cut < 0.) ecut2 = ecut - 2.; 218 elab = K * FragmentA / ResidualA; 219 sig = 0.; 220 if (elab <= ec) { //start for E<Ec 221 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 222 } //end for E<Ec 223 else { //start for E>Ec 224 sig = (landa*elab+mu+nu/elab) * signor; 225 geom = 0.; 226 if (xnulam < flow || elab < etest) return sig; 227 geom = std::sqrt(theA*K); 228 geom = 1.23*ResidualAthrd + ra + 4.573/geom; 229 geom = 31.416 * geom * geom; 230 sig = std::max(geom,sig); 231 } 232 return sig;} 233 234 // ************************** end of cross sections ******************************* 235 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationChannel.cc
r819 r962 26 26 // 27 27 // $Id: G4ProtonEvaporationChannel.cc,v 1.4 2006/06/29 20:10:41 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationProbability.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4ProtonEvaporationProbability.cc,v 1.4 2006/06/29 20:10:43 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara (Nov 1999) 32 // 33 29 // by V. Lara (Oct 1998) 30 // 31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 32 // cross section option 34 33 35 34 #include "G4ProtonEvaporationProbability.hh" 36 35 37 36 G4ProtonEvaporationProbability::G4ProtonEvaporationProbability() : 38 G4EvaporationProbability(1,1,2) // A,Z,Gamma 39 { 40 std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1; 41 std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1; 42 ExcitEnergies.reserve(NumExcitedStatesEnergy); 43 ExcitSpins.reserve(NumExcitedStatesSpin); 44 ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0); 45 ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0); 46 47 48 49 ExcitEnergies[15] = 5.96*MeV; 50 ExcitEnergies[17] = 1.74*MeV; 51 ExcitEnergies[18] = 4.44*MeV; 52 ExcitEnergies[19] = 1.67*MeV; 53 ExcitEnergies[20] = 4.32*MeV; 54 ExcitEnergies[22] = 3.68*MeV; 55 ExcitEnergies[23] = 6.69*MeV; 56 ExcitEnergies[25] = 3.95*MeV; 57 ExcitEnergies[26] = 6.32*MeV; 58 ExcitEnergies[27] = 0.30*MeV; 59 ExcitEnergies[28] = 6.18*MeV; 60 ExcitEnergies[29] = 6.92*MeV; 61 ExcitEnergies[30] = 3.06*MeV; 62 ExcitEnergies[31] = 3.57*MeV; 63 64 65 ExcitSpins[15] = 8; 66 ExcitSpins[17] = 1; 67 ExcitSpins[18] = 6; 68 ExcitSpins[19] = 5; 69 ExcitSpins[20] = 6; 70 ExcitSpins[22] = 4; 71 ExcitSpins[23] = 8; 72 ExcitSpins[25] = 3; 73 ExcitSpins[26] = 4; 74 ExcitSpins[27] = 7; 75 ExcitSpins[28] = 4; 76 ExcitSpins[29] = 5; 77 ExcitSpins[30] = 2; 78 ExcitSpins[31] = 10; 79 80 SetExcitationEnergiesPtr(&ExcitEnergies); 81 SetExcitationSpinsPtr(&ExcitSpins); 82 37 G4EvaporationProbability(1,1,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier 38 { 39 83 40 } 84 41 … … 106 63 } 107 64 108 G4double G4ProtonEvaporationProbability::CCoeficient(const G4double aZ) const 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 & ) 69 { return 0.0; } 70 71 G4double G4ProtonEvaporationProbability::CCoeficient(const G4double aZ) 109 72 { 110 73 // Data comes from … … 126 89 127 90 } 91 92 /////////////////////////////////////////////////////////////////////////////////// 93 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 94 //OPT=0 Dostrovski's parameterization 95 //OPT=1,2 Chatterjee's paramaterization 96 //OPT=3,4 Kalbach's parameterization 97 // 98 G4double G4ProtonEvaporationProbability::CrossSection(const G4Fragment & fragment, const G4double K) 99 { 100 // G4cout<<" In G4ProtonEVaporationProbability OPTxs="<<OPTxs<<G4endl; 101 // G4cout<<" In G4ProtonEVaporationProbability useSICB="<<useSICB<<G4endl; 102 103 theA=GetA(); 104 theZ=GetZ(); 105 ResidualA=fragment.GetA()-theA; 106 ResidualZ=fragment.GetZ()-theZ; 107 108 ResidualAthrd=std::pow(ResidualA,0.33333); 109 FragmentA=fragment.GetA(); 110 FragmentAthrd=std::pow(FragmentA,0.33333); 111 U=fragment.GetExcitationEnergy(); 112 113 114 if (OPTxs==0) {std::ostringstream errOs; 115 errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (protons)!!" <<G4endl; 116 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 117 return 0.;} 118 else if( OPTxs==1 ) return GetOpt1( K); 119 else if( OPTxs==2 ||OPTxs==4) return GetOpt2( K); 120 else if (OPTxs==3 ) return GetOpt3( K); 121 else{ 122 std::ostringstream errOs; 123 errOs << "BAD PROTON CROSS SECTION OPTION AT EVAPORATION!!" <<G4endl; 124 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 125 return 0.; 126 } 127 } 128 //********************* OPT=1 : Chatterjee's cross section ************************ 129 //(fitting to cross section from Bechetti & Greenles OM potential) 130 131 G4double G4ProtonEvaporationProbability::GetOpt1(const G4double K) 132 { 133 G4double Kc=K; 134 135 // JMQ xsec is set constat above limit of validity 136 if (K>50) Kc=50; 137 138 G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs; 139 G4double p, p0, p1, p2,Ec,delta,q,r,ji; 140 141 p0 = 15.72; 142 p1 = 9.65; 143 p2 = -449.0; 144 landa0 = 0.00437; 145 landa1 = -16.58; 146 mu0 = 244.7; 147 mu1 = 0.503; 148 nu0 = 273.1; 149 nu1 = -182.4; 150 nu2 = -1.872; 151 delta=0.; 152 153 Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta); 154 p = p0 + p1/Ec + p2/(Ec*Ec); 155 landa = landa0*ResidualA + landa1; 156 mu = mu0*std::pow(ResidualA,mu1); 157 nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec)); 158 q = landa - nu/(Ec*Ec) - 2*p*Ec; 159 r = mu + 2*nu/Ec + p*(Ec*Ec); 160 161 ji=std::max(Kc,Ec); 162 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;} 163 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); 178 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; 183 184 eekin=K; 185 rnpro=ResidualZ; 186 rnneu=ResidualA-ResidualZ; 187 ekin=eekin/1000; 188 189 r0=1.36*1.e-15; 190 fac=pi*r0*r0; 191 b0=2.247-0.915*(1.-1./ResidualAthrd); 192 fac1=b0*(1.-1./ResidualAthrd); 193 fac2=1.; 194 if(rnneu > 1.5) fac2=std::log(rnneu); 195 xine_th= 1.e+31*fac*fac2*(1.+ResidualAthrd-fac1); 196 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; 200 fac=1.-(1./(1.+std::exp(-8.*ff1*(std::log10(ekin)+1.37*ff2)))); 201 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; 204 fac=-8.*ff1*(std::log10(ekin)+2.0*ff2); 205 fac=1./(1.+std::exp(fac)); 206 xine_th=xine_th*fac; 207 if (xine_th < 0.0){ 208 std::ostringstream errOs; 209 G4cout<<"WARNING: negative Wellisch cross section "<<G4endl; 210 errOs << "RESIDUAL: A=" << ResidualA << " Z=" << ResidualZ <<G4endl; 211 errOs <<" xsec("<<ekin<<" MeV) ="<<xine_th <<G4endl; 212 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 213 } 214 215 return xine_th; 216 217 } 218 219 220 // *********** OPT=3 : Kalbach's cross sections (from PRECO code)************* 221 G4double G4ProtonEvaporationProbability::GetOpt3(const G4double K) 222 { 223 // ** p from becchetti and greenlees (but modified with sub-barrier 224 // ** correction function and xp2 changed from -449) 225 // JMQ (june 2008) : refinement of proton cross section for light systems 226 // 227 G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2; 228 G4double p, p0, p1, p2; 229 p0 = 15.72; 230 p1 = 9.65; 231 p2 = -300.; 232 landa0 = 0.00437; 233 landa1 = -16.58; 234 mu0 = 244.7; 235 mu1 = 0.503; 236 nu0 = 273.1; 237 nu1 = -182.4; 238 nu2 = -1.872; 239 240 // parameters for proton cross section refinement 241 G4double afit,bfit,a2,b2; 242 afit=-0.0785656; 243 bfit=5.10789; 244 a2= -0.00089076; 245 b2= 0.0231597; 246 247 G4double ec,ecsq,xnulam,etest(0.),ra(0.),a,w,c,signor(1.),signor2,sig; 248 G4double b,ecut,cut,ecut2,geom,elab; 249 250 251 G4double flow = 1.e-18; 252 G4double spill= 1.e+18; 253 254 255 if (ResidualA <= 60.) signor = 0.92; 256 else if (ResidualA < 100.) signor = 0.8 + ResidualA*0.002; 257 258 259 ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); 260 ecsq = ec * ec; 261 p = p0 + p1/ec + p2/ecsq; 262 landa = landa0*ResidualA + landa1; 263 a = std::pow(ResidualA,mu1); 264 mu = mu0 * a; 265 nu = a* (nu0+nu1*ec+nu2*ecsq); 266 267 c =std::min(3.15,ec*0.5); 268 w = 0.7 * c / 3.15; 269 270 xnulam = nu / landa; 271 if (xnulam > spill) xnulam=0.; 272 if (xnulam >= flow) etest =std::sqrt(xnulam) + 7.; 273 274 a = -2.*p*ec + landa - nu/ecsq; 275 b = p*ecsq + mu + 2.*nu/ec; 276 ecut = 0.; 277 cut = a*a - 4.*p*b; 278 if (cut > 0.) ecut = std::sqrt(cut); 279 ecut = (ecut-a) / (p+p); 280 ecut2 = ecut; 281 if (cut < 0.) ecut2 = ecut - 2.; 282 elab = K * FragmentA / ResidualA; 283 sig = 0.; 284 if (elab <= ec) { //start for E<Ec 285 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 286 signor2 = (ec-elab-c) / w; 287 signor2 = 1. + std::exp(signor2); 288 sig = sig / signor2; 289 // signor2 is empirical global corr'n at low elab for protons in PRECO, not enough for light targets 290 // refinement for proton cross section 291 if (ResidualZ<=26) 292 sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec)); 293 } //end for E<Ec 294 else { //start for E>Ec 295 sig = (landa*elab+mu+nu/elab) * signor; 296 297 // refinement for proton cross section 298 if ( ResidualZ<=26 && elab <=(afit*ResidualZ+bfit)*ec) 299 sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec)); 300 301 // 302 geom = 0.; 303 304 if (xnulam < flow || elab < etest) 305 { 306 if (sig <0.0) {sig=0.0;} 307 return sig; 308 } 309 geom = std::sqrt(theA*K); 310 geom = 1.23*ResidualAthrd + ra + 4.573/geom; 311 geom = 31.416 * geom * geom; 312 sig = std::max(geom,sig); 313 314 } //end for E>Ec 315 return sig;} 316 317 318 319 // ************************** end of cross sections ******************************* 320 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationChannel.cc
r819 r962 26 26 // 27 27 // $Id: G4TritonEvaporationChannel.cc,v 1.4 2006/06/29 20:10:45 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationProbability.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4TritonEvaporationProbability.cc,v 1.4 2006/06/29 20:10:47 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ 26 //J.M. Quesada (August2008). Based on: 29 27 // 30 28 // Hadronic Process: Nuclear De-excitations 31 // by V. Lara (Nov 1999) 32 // 33 29 // by V. Lara (Oct 1998) 30 // 31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 32 // cross section option 34 33 35 34 #include "G4TritonEvaporationProbability.hh" 36 35 37 36 G4TritonEvaporationProbability::G4TritonEvaporationProbability() : 38 G4EvaporationProbability(3,1,6) // A,Z,Gamma 39 { 40 std::vector<G4double>::size_type NumExcitedStatesEnergy = 31+1; 41 std::vector<G4int>::size_type NumExcitedStatesSpin = 31+1; 42 ExcitEnergies.reserve(NumExcitedStatesEnergy); 43 ExcitSpins.reserve(NumExcitedStatesSpin); 44 ExcitEnergies.insert(ExcitEnergies.begin(),NumExcitedStatesEnergy,0.0); 45 ExcitSpins.insert(ExcitSpins.begin(),NumExcitedStatesSpin,0); 46 47 48 ExcitEnergies[15] = 6.26*MeV; 49 ExcitEnergies[17] = 3.59*MeV; 50 ExcitEnergies[18] = 6.76*MeV; 51 ExcitEnergies[20] = 6.34*MeV; 52 ExcitEnergies[23] = 7.34*MeV; 53 ExcitEnergies[25] = 5.11*MeV; 54 ExcitEnergies[26] = 7.57*MeV; 55 ExcitEnergies[28] = 7.28*MeV; 56 ExcitEnergies[31] = 4.46*MeV; 57 58 ExcitSpins[15] = 5; 59 ExcitSpins[17] = 5; 60 ExcitSpins[18] = 10; 61 ExcitSpins[20] = 2; 62 ExcitSpins[23] = 5; 63 ExcitSpins[25] = 5; 64 ExcitSpins[26] = 8; 65 ExcitSpins[28] = 8; 66 ExcitSpins[31] = 3; 67 68 SetExcitationEnergiesPtr(&ExcitEnergies); 69 SetExcitationSpinsPtr(&ExcitSpins); 37 G4EvaporationProbability(3,1,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier 38 { 39 70 40 } 71 41 … … 96 66 } 97 67 98 G4double G4TritonEvaporationProbability::CCoeficient(const G4double aZ) const 68 G4double G4TritonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 69 { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));} 70 71 G4double G4TritonEvaporationProbability::CalcBetaParam(const G4Fragment & ) 72 { return 0.0; } 73 74 G4double G4TritonEvaporationProbability::CCoeficient(const G4double aZ) 99 75 { 100 76 // Data comes from … … 117 93 118 94 } 95 96 /////////////////////////////////////////////////////////////////////////////////// 97 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections 98 //OPT=0 Dostrovski's parameterization 99 //OPT=1,2 Chatterjee's paramaterization 100 //OPT=3,4 Kalbach's parameterization 101 // 102 G4double G4TritonEvaporationProbability::CrossSection(const G4Fragment & fragment, const G4double K) 103 { 104 theA=GetA(); 105 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); 112 113 114 if (OPTxs==0) {std::ostringstream errOs; 115 errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (tritons)!!" 116 <<G4endl; 117 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 118 return 0.;} 119 if( OPTxs==1 || OPTxs==2) return G4TritonEvaporationProbability::GetOpt12( K); 120 else if (OPTxs==3 || OPTxs==4) return G4TritonEvaporationProbability::GetOpt34( K); 121 else{ 122 std::ostringstream errOs; 123 errOs << "BAD Triton CROSS SECTION OPTION AT EVAPORATION!!" <<G4endl; 124 throw G4HadronicException(__FILE__, __LINE__, errOs.str()); 125 return 0.; 126 } 127 } 128 129 // 130 //********************* OPT=1,2 : Chatterjee's cross section ************************ 131 //(fitting to cross section from Bechetti & Greenles OM potential) 132 133 G4double G4TritonEvaporationProbability::GetOpt12(const G4double K) 134 { 135 136 G4double Kc=K; 137 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 144 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 173 } 174 175 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)************* 176 G4double G4TritonEvaporationProbability::GetOpt34(const G4double K) 177 // ** t from o.m. of hafele, flynn et al 178 { 179 180 G4double landa, mu, nu, p , signor(1.),sig; 181 G4double ec,ecsq,xnulam,etest(0.),a; 182 G4double b,ecut,cut,ecut2,geom,elab; 183 184 185 G4double flow = 1.e-18; 186 G4double spill= 1.e+18; 187 188 G4double p0 = -21.45; 189 G4double p1 = 484.7; 190 G4double p2 = -1608.; 191 G4double landa0 = 0.0186; 192 G4double landa1 = -8.90; 193 G4double mu0 = 686.3; 194 G4double mu1 = 0.325; 195 G4double nu0 = 368.9; 196 G4double nu1 = -522.2; 197 G4double nu2 = -4.998; 198 199 G4double ra=0.80; 200 201 //JMQ 13/02/09 increase of reduced radius to lower the barrier 202 // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra); 203 ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra); 204 ecsq = ec * ec; 205 p = p0 + p1/ec + p2/ecsq; 206 landa = landa0*ResidualA + landa1; 207 a = std::pow(ResidualA,mu1); 208 mu = mu0 * a; 209 nu = a* (nu0+nu1*ec+nu2*ecsq); 210 xnulam = nu / landa; 211 if (xnulam > spill) xnulam=0.; 212 if (xnulam >= flow) etest = 1.2 *std::sqrt(xnulam); 213 214 a = -2.*p*ec + landa - nu/ecsq; 215 b = p*ecsq + mu + 2.*nu/ec; 216 ecut = 0.; 217 cut = a*a - 4.*p*b; 218 if (cut > 0.) ecut = std::sqrt(cut); 219 ecut = (ecut-a) / (p+p); 220 ecut2 = ecut; 221 if (cut < 0.) ecut2 = ecut - 2.; 222 elab = K * FragmentA / ResidualA; 223 sig = 0.; 224 225 if (elab <= ec) { //start for E<Ec 226 if (elab > ecut2) sig = (p*elab*elab+a*elab+b) * signor; 227 } //end for E<Ec 228 else { //start for E>Ec 229 sig = (landa*elab+mu+nu/elab) * signor; 230 geom = 0.; 231 if (xnulam < flow || elab < etest) return sig; 232 geom = std::sqrt(theA*K); 233 geom = 1.23*ResidualAthrd + ra + 4.573/geom; 234 geom = 31.416 * geom * geom; 235 sig = std::max(geom,sig); 236 } //end for E>Ec 237 return sig; 238 239 } 240 241 // ************************** end of cross sections ******************************* 242 -
trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4VEvaporation.cc
r819 r962 26 26 // 27 27 // $Id: G4VEvaporation.cc,v 1.5 2006/06/29 20:10:49 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Hadronic Process: Nuclear De-excitations
Note: See TracChangeset
for help on using the changeset viewer.