00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "globals.hh"
00043 #include "G4ios.hh"
00044 #include "G4KineticTrack.hh"
00045 #include "G4VXResonance.hh"
00046 #include "Randomize.hh"
00047 #include "G4Proton.hh"
00048 #include "G4HadTmpUtil.hh"
00049
00050 G4VXResonance::G4VXResonance()
00051 { }
00052
00053
00054 G4VXResonance::~G4VXResonance()
00055 { }
00056
00057
00058 G4bool G4VXResonance::operator==(const G4VXResonance &right) const
00059 {
00060 return (this == (G4VXResonance *) &right);
00061 }
00062
00063
00064 G4bool G4VXResonance::operator!=(const G4VXResonance &right) const
00065 {
00066 return (this != (G4VXResonance *) &right);
00067 }
00068
00069
00070 G4double G4VXResonance::IsospinCorrection(const G4KineticTrack& trk1,
00071 const G4KineticTrack& trk2,
00072 G4int isoOut1, G4int isoOut2,
00073 G4double , G4double ) const
00074 {
00075 G4double result = 0.;
00076
00077 G4ParticleDefinition* in1 = trk1.GetDefinition();
00078 G4ParticleDefinition* in2 = trk2.GetDefinition();
00079
00080 G4int isoIn1 = in1->GetPDGiIsospin();
00081 G4int iso3In1 = in1->GetPDGiIsospin3();
00082 G4int isoIn2 = in2->GetPDGiIsospin();
00083 G4int iso3In2 = in2->GetPDGiIsospin3();
00084
00085 G4int isoProton = G4Proton::ProtonDefinition()->GetPDGiIsospin();
00086 G4int iso3Proton = G4Proton::ProtonDefinition()->GetPDGiIsospin3();
00087
00088 G4double pWeight = clebsch.Weight(isoProton,iso3Proton, isoProton,iso3Proton, isoOut1,isoOut2);
00089 if (pWeight == 0.) throw G4HadronicException(__FILE__, __LINE__, "G4VXResonance::IsospinCorrection, no resonances - pWeight is zero");
00090
00091 if (in1->IsShortLived() || in2->IsShortLived())
00092 {
00093
00094 G4int iSpinProton = G4Proton::ProtonDefinition()->GetPDGiSpin();
00095 G4double degeneracyFactor = DegeneracyFactor(trk1,trk2,iSpinProton,iSpinProton);
00096
00097 G4double factor = degeneracyFactor * pWeight;
00098 if (factor > DBL_MIN)
00099 {
00100
00101 std::vector<G4double> iso = clebsch.GenerateIso3(isoIn1,iso3In1,
00102 isoIn2,iso3In2,
00103 isoProton,isoProton);
00104 G4int isoA = G4lrint(iso[0]);
00105 G4int isoB = G4lrint(iso[1]);
00106 G4double rWeight = clebsch.Weight(isoProton,isoA,
00107 isoProton,isoB,
00108 isoOut1,isoOut2);
00109 result = rWeight / pWeight;
00110 }
00111 }
00112 else
00113 {
00114 G4double weight = clebsch.Weight(isoIn1,iso3In1, isoIn2,iso3In2, isoOut1,isoOut2);
00115 result = weight / pWeight;
00116 }
00117
00118 return result;
00119 }
00120
00121
00122 #include "G4DetailedBalancePhaseSpaceIntegral.hh"
00123
00124 G4double G4VXResonance::DetailedBalance(const G4KineticTrack& trk1,
00125 const G4KineticTrack& trk2,
00126 G4int isoOut1, G4int isoOut2,
00127 G4double iSpinOut1, G4double iSpinOut2,
00128 G4double mOut1, G4double mOut2) const
00129 {
00130
00131
00132
00133
00134
00135 G4ParticleDefinition* in1 = trk1.GetDefinition();
00136 G4ParticleDefinition* in2 = trk2.GetDefinition();
00137 if(in1->IsShortLived() && in2->IsShortLived())
00138 {
00139 throw G4HadronicException(__FILE__, __LINE__, "Detailed balance for resonance scattering still on the schedule.");
00140 }
00141
00142 G4double result = 0.;
00143
00144 G4int isoIn1 = in1->GetPDGiIsospin();
00145 G4int iso3In1 = in1->GetPDGiIsospin3();
00146 G4int isoIn2 = in2->GetPDGiIsospin();
00147 G4int iso3In2 = in2->GetPDGiIsospin3();
00148 G4double weight = clebsch.Weight(isoIn1, iso3In1, isoIn2, iso3In2, isoOut1, isoOut2);
00149
00150 if (weight > 00001)
00151 {
00152
00153
00154 G4double degeneracy = DegeneracyFactor(trk1,trk2,iSpinOut1,iSpinOut2);
00155 G4double factor = degeneracy * weight;
00156
00157
00158 G4double S = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag2();
00159 G4double m_1 = in1->GetPDGMass();
00160 G4double m_2 = in2->GetPDGMass();
00161
00162
00163 G4double pinitial2 = (S - (m_1+m_2) * (m_1+m_2)) * (S - (m_1-m_2) * (m_1-m_2)) / (4.0*S);
00164 G4double pfinal2 = (S - (mOut1+mOut2) * (mOut1+mOut2)) * (S - ( mOut1-mOut2) * (mOut1-mOut2)) / (4.0*S);
00165 G4double relativeMomsquared = pfinal2/pinitial2;
00166
00167
00168 if(in1->IsShortLived())
00169 {
00170 G4DetailedBalancePhaseSpaceIntegral theI(in1);
00171 relativeMomsquared = 1./theI.GetPhaseSpaceIntegral(std::sqrt(S));
00172 }
00173 else if(in2->IsShortLived())
00174 {
00175 G4DetailedBalancePhaseSpaceIntegral theI(in2);
00176 relativeMomsquared = 1./theI.GetPhaseSpaceIntegral(std::sqrt(S));
00177 }
00178
00179 result = factor * relativeMomsquared;
00180 }
00181
00182 return result;
00183 }
00184
00185
00186 G4double G4VXResonance::DegeneracyFactor(const G4KineticTrack& trk1,
00187 const G4KineticTrack& trk2,
00188 G4double iSpinOut1, G4double iSpinOut2) const
00189 {
00190 G4double value = 0.;
00191
00192 G4ParticleDefinition* in1 = trk1.GetDefinition();
00193 G4ParticleDefinition* in2 = trk2.GetDefinition();
00194
00195 G4double sIn1 = in1->GetPDGiSpin() + 1.;
00196 G4double sIn2 = in2->GetPDGiSpin() + 1.;
00197
00198 G4double denom = sIn1 * sIn2;
00199 if (denom > 0.)
00200 {
00201 value = (iSpinOut1+1) * (iSpinOut2+1) / denom;
00202 }
00203 return value;
00204 }
00205
00206