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 #include "globals.hh"
00037 #include "G4ios.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "G4XPDGElastic.hh"
00040 #include "G4KineticTrack.hh"
00041 #include "G4ParticleDefinition.hh"
00042 #include "G4DataVector.hh"
00043
00044 #include "G4AntiProton.hh"
00045 #include "G4AntiNeutron.hh"
00046 #include "G4Proton.hh"
00047 #include "G4Neutron.hh"
00048 #include "G4PionPlus.hh"
00049 #include "G4PionMinus.hh"
00050 #include "G4KaonMinus.hh"
00051 #include "G4KaonPlus.hh"
00052
00053 const G4double G4XPDGElastic::_lowLimit = 5. * GeV;
00054 const G4double G4XPDGElastic::_highLimit = DBL_MAX;
00055
00056
00057
00058 const G4int G4XPDGElastic::nPar = 7;
00059
00060 const G4double G4XPDGElastic::pPiPlusPDGFit[7] = { 2., 200., 0., 11.4, -0.4, 0.079, 0. };
00061
00062 const G4double G4XPDGElastic::pPiMinusPDGFit[7] = { 2., 360., 1.76, 11.2, -0.64, 0.043, 0. };
00063
00064 const G4double G4XPDGElastic::pKPlusPDGFit[7] = { 2., 175., 5.0, 8.1, -1.8, 0.16, -1.3 };
00065
00066 const G4double G4XPDGElastic::pKMinusPDGFit[7] = { 2., 175., 7.3, 0., 0., 0.29, -2.40 };
00067
00068 const G4double G4XPDGElastic::ppPDGFit[7] = { 2., 2100., 11.9, 26.9, -1.21, 0.169, -1.85 };
00069
00070 const G4double G4XPDGElastic::ppbarPDGFit[7] = { 5., 1730000., 10.2, 52.7, -1.16, 0.125, -1.28 };
00071
00072 const G4double G4XPDGElastic::npbarPDGFit[7] = { 1.1, 5.55, 36.5, 0., 0., 0., -11.9 };
00073
00074
00075 G4XPDGElastic::G4XPDGElastic()
00076 {
00077 G4ParticleDefinition * proton = G4Proton::ProtonDefinition();
00078 G4ParticleDefinition * neutron = G4Neutron::NeutronDefinition();
00079 G4ParticleDefinition * piPlus = G4PionPlus::PionPlusDefinition();
00080 G4ParticleDefinition * piMinus = G4PionMinus::PionMinusDefinition();
00081 G4ParticleDefinition * KPlus = G4KaonPlus::KaonPlusDefinition();
00082 G4ParticleDefinition * KMinus = G4KaonMinus::KaonMinusDefinition();
00083 G4ParticleDefinition * antiproton = G4AntiProton::AntiProtonDefinition();
00084
00085 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> pp(proton,proton);
00086 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> pn(proton,neutron);
00087 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> piPlusp(piPlus,proton);
00088 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> piMinusp(piMinus,proton);
00089 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> KPlusp(KPlus,proton);
00090 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> KMinusp(KMinus,proton);
00091 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> nn(neutron,neutron);
00092 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> ppbar(proton,antiproton);
00093 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> npbar(antiproton,neutron);
00094
00095 std::vector<G4double> ppData;
00096 std::vector<G4double> pPiPlusData;
00097 std::vector<G4double> pPiMinusData;
00098 std::vector<G4double> pKPlusData;
00099 std::vector<G4double> pKMinusData;
00100 std::vector<G4double> ppbarData;
00101 std::vector<G4double> npbarData;
00102
00103 G4int i;
00104 for (i=0; i<2; i++)
00105 {
00106 ppData.push_back(ppPDGFit[i] * GeV);
00107 pPiPlusData.push_back(pPiPlusPDGFit[i] * GeV);
00108 pPiMinusData.push_back(pPiMinusPDGFit[i] * GeV);
00109 pKPlusData.push_back(pKPlusPDGFit[i] * GeV);
00110 pKMinusData.push_back(pKMinusPDGFit[i] * GeV);
00111 ppbarData.push_back(ppbarPDGFit[i] * GeV);
00112 npbarData.push_back(npbarPDGFit[i] * GeV);
00113 }
00114
00115 for (i=2; i<nPar; i++)
00116 {
00117 ppData.push_back(ppPDGFit[i]);
00118 pPiPlusData.push_back(pPiPlusPDGFit[i]);
00119 pPiMinusData.push_back(pPiMinusPDGFit[i]);
00120 pKPlusData.push_back(pKPlusPDGFit[i]);
00121 pKMinusData.push_back(pKMinusPDGFit[i]);
00122 ppbarData.push_back(ppbarPDGFit[i]);
00123 npbarData.push_back(npbarPDGFit[i]);
00124 }
00125
00126 xMap[nn] = ppData;
00127 xMap[pp] = ppData;
00128 xMap[pn] = ppData;
00129 xMap[piPlusp] = pPiPlusData;
00130 xMap[piMinusp] = pPiMinusData;
00131 xMap[KPlusp] = pKPlusData;
00132 xMap[KMinusp] = pKMinusData;
00133 xMap[ppbar] = ppbarData;
00134 xMap[npbar] = npbarData;
00135 }
00136
00137
00138 G4XPDGElastic::~G4XPDGElastic()
00139 { }
00140
00141
00142 G4bool G4XPDGElastic::operator==(const G4XPDGElastic &right) const
00143 {
00144 return (this == (G4XPDGElastic *) &right);
00145 }
00146
00147
00148 G4bool G4XPDGElastic::operator!=(const G4XPDGElastic &right) const
00149 {
00150 return (this != (G4XPDGElastic *) &right);
00151 }
00152
00153
00154 G4double G4XPDGElastic::CrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const
00155 {
00156
00157
00158 G4double sigma = 0.;
00159
00160 G4ParticleDefinition* def1 = trk1.GetDefinition();
00161 G4ParticleDefinition* def2 = trk2.GetDefinition();
00162
00163 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
00164 G4double m_1 = def1->GetPDGMass();
00165 G4double m_2 = def2->GetPDGMass();
00166 G4double m_max = std::max(m_1,m_2);
00167
00168
00169 G4double pLab = 0.;
00170
00171 if (m_max > 0. && sqrtS > (m_1 + m_2))
00172 {
00173 pLab = std::sqrt( (sqrtS*sqrtS - (m_1+m_2)*(m_1+m_2) ) * (sqrtS*sqrtS - (m_1-m_2)*(m_1-m_2)) ) / (2*m_max);
00174
00175
00176
00177
00178 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> trkPair(def1,def2);
00179 if (def1->GetPDGMass() > def2->GetPDGMass())
00180 trkPair = std::pair<G4ParticleDefinition *,G4ParticleDefinition *>(def2,def1);
00181
00182 std::vector<G4double> data;
00183 G4double pMinFit = 0.;
00184 G4double pMaxFit = 0.;
00185 G4double aFit = 0.;
00186 G4double bFit = 0.;
00187 G4double cFit = 0.;
00188 G4double dFit = 0.;
00189 G4double nFit = 0.;
00190
00191
00192
00193
00194
00195 if (xMap.find(trkPair) != xMap.end())
00196 {
00197 PairDoubleMap::const_iterator iter;
00198 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
00199 {
00200 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> thePair = (*iter).first;
00201 if (thePair == trkPair)
00202 {
00203 data = (*iter).second;
00204 pMinFit = data[0];
00205 pMaxFit = data[1];
00206 aFit = data[2];
00207 bFit = data[3];
00208 cFit = data[5];
00209 dFit = data[6];
00210 nFit = data[4];
00211
00212 if (pLab < pMinFit) return 0.0;
00213 if (pLab > pMaxFit )
00214 G4cout << "WARNING! G4XPDGElastic::PDGElastic "
00215 << trk1.GetDefinition()->GetParticleName() << "-"
00216 << trk2.GetDefinition()->GetParticleName()
00217 << " elastic cross section: momentum "
00218 << pLab / GeV << " GeV outside valid fit range "
00219 << pMinFit /GeV << "-" << pMaxFit / GeV
00220 << G4endl;
00221
00222 pLab /= GeV;
00223 if (pLab > 0.)
00224 {
00225 G4double logP = std::log(pLab);
00226 sigma = aFit + bFit * std::pow(pLab, nFit) + cFit * logP*logP + dFit * logP;
00227 sigma = sigma * millibarn;
00228 }
00229 }
00230 }
00231 }
00232 else
00233 {
00234 G4cout << "G4XPDGElastic::CrossSection - Data not found in Map" << G4endl;
00235 }
00236 }
00237
00238 if (sigma < 0.)
00239 {
00240 G4cout << "WARNING! G4XPDGElastic::PDGElastic "
00241 << def1->GetParticleName() << "-" << def2->GetParticleName()
00242 << " elastic cross section: momentum "
00243 << pLab << " GeV, negative cross section "
00244 << sigma / millibarn << " mb set to 0" << G4endl;
00245 sigma = 0.;
00246 }
00247
00248 return sigma;
00249 }
00250
00251
00252 G4String G4XPDGElastic::Name() const
00253 {
00254 G4String name = "PDGElastic ";
00255 return name;
00256 }
00257
00258
00259 G4bool G4XPDGElastic::IsValid(G4double e) const
00260 {
00261 G4bool answer = InLimits(e,_lowLimit,_highLimit);
00262
00263 return answer;
00264 }
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275