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 #include "globals.hh"
00031 #include "G4ios.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "G4XPDGTotal.hh"
00034 #include "G4KineticTrack.hh"
00035 #include "G4ParticleDefinition.hh"
00036 #include "G4DataVector.hh"
00037 #include "G4AntiProton.hh"
00038 #include "G4AntiNeutron.hh"
00039 #include "G4Proton.hh"
00040 #include "G4Neutron.hh"
00041 #include "G4PionPlus.hh"
00042 #include "G4PionMinus.hh"
00043 #include "G4Gamma.hh"
00044 #include "G4KaonMinus.hh"
00045 #include "G4KaonPlus.hh"
00046
00047 const G4double G4XPDGTotal::_lowLimit = 3. * GeV;
00048 const G4double G4XPDGTotal::_highLimit = DBL_MAX;
00049
00050
00051
00052 const G4int G4XPDGTotal::nFit = 5;
00053
00054 const G4double G4XPDGTotal::ppPDGFit[5] = { 3., 40000., 18.256, 60.19, 33.43 };
00055
00056 const G4double G4XPDGTotal::npPDGFit[5] = { 3., 40., 18.256, 61.14, 29.80 };
00057
00058 const G4double G4XPDGTotal::pipPDGFit[5] = { 3., 40., 11.568, 27.55, 5.62 };
00059
00060 const G4double G4XPDGTotal::KpPDGFit[5] = { 3., 40., 10.376, 15.57, 13.19 };
00061
00062 const G4double G4XPDGTotal::KnPDGFit[5] = { 3., 40., 10.376, 14.29, 7.38 };
00063
00064 const G4double G4XPDGTotal::gammapPDGFit[5] = { 3., 300., 0.0577, 0.1171, 0. };
00065
00066 const G4double G4XPDGTotal::gammagammaPDGFit[5] = { 3., 300., 0.000156, 0.00032, 0. };
00067
00068
00069 G4XPDGTotal::G4XPDGTotal()
00070 {
00071 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> pp(G4Proton::ProtonDefinition(),
00072 G4Proton::ProtonDefinition());
00073 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> pn(G4Proton::ProtonDefinition(),
00074 G4Neutron::NeutronDefinition());
00075 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> piPlusp(G4PionPlus::PionPlusDefinition(),
00076 G4Proton::ProtonDefinition());
00077 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> piMinusp(G4PionMinus::PionMinusDefinition(),
00078 G4Proton::ProtonDefinition());
00079 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> KPlusp(G4KaonPlus::KaonPlusDefinition(),
00080 G4Proton::ProtonDefinition());
00081 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> KPlusn(G4KaonPlus::KaonPlusDefinition(),
00082 G4Neutron::NeutronDefinition());
00083 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> KMinusp(G4KaonMinus::KaonMinusDefinition(),
00084 G4Proton::ProtonDefinition());
00085 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> KMinusn(G4KaonMinus::KaonMinusDefinition(),
00086 G4Neutron::NeutronDefinition());
00087 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> gp(G4Gamma::GammaDefinition(),
00088 G4Proton::ProtonDefinition());
00089 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> gg(G4Gamma::GammaDefinition(),
00090 G4Gamma::GammaDefinition());
00091 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> nn(G4Neutron::NeutronDefinition(),
00092 G4Neutron::NeutronDefinition());
00093
00094 std::vector<G4double> nnData;
00095 std::vector<G4double> ppData;
00096 std::vector<G4double> pnData;
00097 std::vector<G4double> pipData;
00098 std::vector<G4double> KpData;
00099 std::vector<G4double> KnData;
00100 std::vector<G4double> gpData;
00101 std::vector<G4double> ggData;
00102
00103 G4int i;
00104 for (i=0; i<2; i++)
00105 {
00106 nnData.push_back(ppPDGFit[i] * GeV);
00107 ppData.push_back(ppPDGFit[i] * GeV);
00108 pnData.push_back(npPDGFit[i] * GeV);
00109 pipData.push_back(pipPDGFit[i] * GeV);
00110 KpData.push_back(KpPDGFit[i] * GeV);
00111 KnData.push_back(KnPDGFit[i] * GeV);
00112 gpData.push_back(gammapPDGFit[i] * GeV);
00113 ggData.push_back(gammagammaPDGFit[i] * GeV);
00114 }
00115 for (i=2; i<nFit; i++)
00116 {
00117 nnData.push_back(ppPDGFit[i]);
00118 ppData.push_back(ppPDGFit[i]);
00119 pnData.push_back(npPDGFit[i]);
00120 pipData.push_back(pipPDGFit[i]);
00121 KpData.push_back(KpPDGFit[i]);
00122 KnData.push_back(KnPDGFit[i]);
00123 gpData.push_back(gammapPDGFit[i]);
00124 ggData.push_back(gammagammaPDGFit[i]);
00125 }
00126
00127 xMap[pp] = ppData;
00128 xMap[pn] = pnData;
00129 xMap[piPlusp] = pipData;
00130 xMap[piMinusp] = pipData;
00131 xMap[KPlusp] = KpData;
00132 xMap[KPlusn] = KnData;
00133 xMap[KMinusp] = KpData;
00134 xMap[KMinusn] = KnData;
00135 xMap[gp] = gpData;
00136 xMap[gg] = ggData;
00137 xMap[nn] = nnData;
00138 }
00139
00140
00141 G4XPDGTotal::~G4XPDGTotal()
00142 { }
00143
00144
00145 G4bool G4XPDGTotal::operator==(const G4XPDGTotal &right) const
00146 {
00147 return (this == (G4XPDGTotal *) &right);
00148 }
00149
00150
00151 G4bool G4XPDGTotal::operator!=(const G4XPDGTotal &right) const
00152 {
00153 return (this != (G4XPDGTotal *) &right);
00154 }
00155
00156
00157 G4double G4XPDGTotal::CrossSection(const G4KineticTrack& trk1,
00158 const G4KineticTrack& trk2) const
00159 {
00160 G4double sigma = 0.;
00161
00162 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
00163
00164 G4ParticleDefinition* def1 = trk1.GetDefinition();
00165 G4ParticleDefinition* def2 = trk2.GetDefinition();
00166
00167 G4double enc1 = def1->GetPDGEncoding();
00168 G4double enc2 = def2->GetPDGEncoding();
00169 G4double coeff = -1.;
00170 if ( (enc1 < 0 && enc2 >0) || (enc2 < 0 && enc1 >0) ) coeff = 1.;
00171
00172
00173 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> trkPair(def1,def2);
00174
00175 if (def1->GetPDGMass() > def2->GetPDGMass())
00176 trkPair = std::pair<G4ParticleDefinition *,G4ParticleDefinition *>(def2,def1);
00177
00178 std::vector<G4double> data;
00179
00180 if (xMap.find(trkPair) != xMap.end())
00181 {
00182
00183 PairDoubleMap::const_iterator iter;
00184 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
00185 {
00186 std::pair<G4ParticleDefinition *,G4ParticleDefinition *> thePair = (*iter).first;
00187 if (thePair == trkPair)
00188 {
00189 data = (*iter).second;
00190
00191 G4double eMinFit = data[0];
00192 G4double eMaxFit = data[1];
00193 G4double xFit = data[2];
00194 G4double y1Fit = data[3];
00195 G4double y2Fit = data[4];
00196
00197
00198
00199
00200 const G4double epsilon = 0.095;
00201 const G4double eta1 = -0.34;
00202 const G4double eta2 = -0.55;
00203
00204 if (sqrtS < eMinFit || sqrtS > eMaxFit)
00205 {
00206 G4cout << "WARNING! G4XPDGTotal::PDGTotal extrapolating cross section at "
00207 << sqrtS / GeV
00208 << " GeV outside the PDG fit range "
00209 << eMinFit / GeV << " - " << eMaxFit / GeV << " GeV " << G4endl;
00210 }
00211
00212 G4double S = (sqrtS * sqrtS) / (GeV*GeV);
00213
00214 sigma = ( (xFit * std::pow(S,epsilon)) +
00215 (y1Fit * std::pow(S,eta1)) +
00216 (coeff * y2Fit * std::pow(S,eta2)) ) * millibarn;
00217
00218 if (sigma < 0.)
00219 {
00220 G4String name1 = def1->GetParticleName();
00221 G4String name2 = def2->GetParticleName();
00222 G4cout << "WARNING! G4XPDGTotal::PDGTotal "
00223 << name1 << "-" << name2
00224 << " total cross section: Ecm "
00225 << sqrtS / GeV << " GeV, negative cross section "
00226 << sigma / millibarn << " mb set to 0" << G4endl;
00227 sigma = 0.;
00228 }
00229 }
00230 }
00231 }
00232 return sigma;
00233 }
00234
00235
00236 G4String G4XPDGTotal::Name() const
00237 {
00238 G4String name = "PDGTotal ";
00239 return name;
00240 }
00241
00242
00243 G4bool G4XPDGTotal::IsValid(G4double e) const
00244 {
00245 G4bool answer = InLimits(e,_lowLimit,_highLimit);
00246
00247 return answer;
00248 }
00249
00250
00251 G4double G4XPDGTotal::PDGTotal(G4double ,G4double ) const
00252 {
00253 return 0.;
00254 }
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265