Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // >> 27 // $Id: G4XPDGElastic.cc,v 1.3 2006/06/29 20:42:42 gunter Exp $ // 27 // ------------------------------------------- 28 // ------------------------------------------------------------------- 28 // 29 // 29 // PDG Elastic cross section 30 // PDG Elastic cross section 30 // PDG fits, Phys.Rev. D50 (1994), 1335 31 // PDG fits, Phys.Rev. D50 (1994), 1335 31 // 32 // 32 // 33 // 33 // ------------------------------------------- 34 // ------------------------------------------------------------------- 34 35 35 #include "globals.hh" 36 #include "globals.hh" 36 #include "G4ios.hh" 37 #include "G4ios.hh" 37 #include "G4Log.hh" << 38 #include "G4Pow.hh" << 39 #include "G4SystemOfUnits.hh" << 40 #include "G4XPDGElastic.hh" 38 #include "G4XPDGElastic.hh" 41 #include "G4KineticTrack.hh" 39 #include "G4KineticTrack.hh" 42 #include "G4ParticleDefinition.hh" 40 #include "G4ParticleDefinition.hh" 43 #include "G4DataVector.hh" 41 #include "G4DataVector.hh" 44 42 45 #include "G4AntiProton.hh" 43 #include "G4AntiProton.hh" 46 #include "G4AntiNeutron.hh" 44 #include "G4AntiNeutron.hh" 47 #include "G4Proton.hh" 45 #include "G4Proton.hh" 48 #include "G4Neutron.hh" 46 #include "G4Neutron.hh" 49 #include "G4PionPlus.hh" 47 #include "G4PionPlus.hh" 50 #include "G4PionMinus.hh" 48 #include "G4PionMinus.hh" 51 #include "G4KaonMinus.hh" 49 #include "G4KaonMinus.hh" 52 #include "G4KaonPlus.hh" 50 #include "G4KaonPlus.hh" 53 51 54 const G4double G4XPDGElastic::_lowLimit = 5. * 52 const G4double G4XPDGElastic::_lowLimit = 5. * GeV; 55 const G4double G4XPDGElastic::_highLimit = DBL 53 const G4double G4XPDGElastic::_highLimit = DBL_MAX; 56 54 57 // Parameters of the PDG Elastic cross-section 55 // Parameters of the PDG Elastic cross-section fit (Rev. Particle Properties, 1998) 58 // Columns are: lower and higher fit range, X, 56 // Columns are: lower and higher fit range, X, Y1, Y2 59 const G4int G4XPDGElastic::nPar = 7; 57 const G4int G4XPDGElastic::nPar = 7; 60 // p pi+ 58 // p pi+ 61 const G4double G4XPDGElastic::pPiPlusPDGFit[7] 59 const G4double G4XPDGElastic::pPiPlusPDGFit[7] = { 2., 200., 0., 11.4, -0.4, 0.079, 0. }; 62 // p pi- 60 // p pi- 63 const G4double G4XPDGElastic::pPiMinusPDGFit[7 61 const G4double G4XPDGElastic::pPiMinusPDGFit[7] = { 2., 360., 1.76, 11.2, -0.64, 0.043, 0. }; 64 // p K+ 62 // p K+ 65 const G4double G4XPDGElastic::pKPlusPDGFit[7] 63 const G4double G4XPDGElastic::pKPlusPDGFit[7] = { 2., 175., 5.0, 8.1, -1.8, 0.16, -1.3 }; 66 // p K- 64 // p K- 67 const G4double G4XPDGElastic::pKMinusPDGFit[7] 65 const G4double G4XPDGElastic::pKMinusPDGFit[7] = { 2., 175., 7.3, 0., 0., 0.29, -2.40 }; 68 // p p 66 // p p 69 const G4double G4XPDGElastic::ppPDGFit[7] = 67 const G4double G4XPDGElastic::ppPDGFit[7] = { 2., 2100., 11.9, 26.9, -1.21, 0.169, -1.85 }; 70 // p pbar 68 // p pbar 71 const G4double G4XPDGElastic::ppbarPDGFit[7] = 69 const G4double G4XPDGElastic::ppbarPDGFit[7] = { 5., 1730000., 10.2, 52.7, -1.16, 0.125, -1.28 }; 72 // n pbar 70 // n pbar 73 const G4double G4XPDGElastic::npbarPDGFit[7] = 71 const G4double G4XPDGElastic::npbarPDGFit[7] = { 1.1, 5.55, 36.5, 0., 0., 0., -11.9 }; 74 72 75 73 76 G4XPDGElastic::G4XPDGElastic() 74 G4XPDGElastic::G4XPDGElastic() 77 { 75 { 78 const G4ParticleDefinition * proton = G4Prot << 76 G4String proton = G4Proton::ProtonDefinition()->GetParticleName(); 79 const G4ParticleDefinition * neutron = G4Neu << 77 G4String neutron = G4Neutron::NeutronDefinition()->GetParticleName(); 80 const G4ParticleDefinition * piPlus = G4Pion << 78 G4String piPlus = G4PionPlus::PionPlusDefinition()->GetParticleName(); 81 const G4ParticleDefinition * piMinus = G4Pio << 79 G4String piMinus = G4PionMinus::PionMinusDefinition()->GetParticleName(); 82 const G4ParticleDefinition * KPlus = G4KaonP << 80 G4String KPlus = G4KaonPlus::KaonPlusDefinition()->GetParticleName(); 83 const G4ParticleDefinition * KMinus = G4Kaon << 81 G4String KMinus = G4KaonMinus::KaonMinusDefinition()->GetParticleName(); 84 const G4ParticleDefinition * antiproton = G4 << 82 G4String antiproton = G4AntiProton::AntiProtonDefinition()->GetParticleName(); 85 83 86 std::pair<const G4ParticleDefinition *,const << 84 std::pair<G4String,G4String> pp(proton,proton); 87 std::pair<const G4ParticleDefinition *,const << 85 std::pair<G4String,G4String> pn(proton,neutron); 88 std::pair<const G4ParticleDefinition *,const << 86 std::pair<G4String,G4String> piPlusp(piPlus,proton); 89 std::pair<const G4ParticleDefinition *,const << 87 std::pair<G4String,G4String> piMinusp(piMinus,proton); 90 std::pair<const G4ParticleDefinition *,const << 88 std::pair<G4String,G4String> KPlusp(KPlus,proton); 91 std::pair<const G4ParticleDefinition *,const << 89 std::pair<G4String,G4String> KMinusp(KMinus,proton); 92 std::pair<const G4ParticleDefinition *,const << 90 std::pair<G4String,G4String> nn(neutron,neutron); 93 std::pair<const G4ParticleDefinition *,const << 91 std::pair<G4String,G4String> ppbar(proton,antiproton); 94 std::pair<const G4ParticleDefinition *,const << 92 std::pair<G4String,G4String> npbar(antiproton,neutron); 95 93 96 std::vector<G4double> ppData; 94 std::vector<G4double> ppData; 97 std::vector<G4double> pPiPlusData; 95 std::vector<G4double> pPiPlusData; 98 std::vector<G4double> pPiMinusData; 96 std::vector<G4double> pPiMinusData; 99 std::vector<G4double> pKPlusData; 97 std::vector<G4double> pKPlusData; 100 std::vector<G4double> pKMinusData; 98 std::vector<G4double> pKMinusData; 101 std::vector<G4double> ppbarData; 99 std::vector<G4double> ppbarData; 102 std::vector<G4double> npbarData; 100 std::vector<G4double> npbarData; 103 101 104 G4int i; 102 G4int i; 105 for (i=0; i<2; i++) 103 for (i=0; i<2; i++) 106 { 104 { 107 ppData.push_back(ppPDGFit[i] * GeV); 105 ppData.push_back(ppPDGFit[i] * GeV); 108 pPiPlusData.push_back(pPiPlusPDGFit[i] * 106 pPiPlusData.push_back(pPiPlusPDGFit[i] * GeV); 109 pPiMinusData.push_back(pPiMinusPDGFit[i] 107 pPiMinusData.push_back(pPiMinusPDGFit[i] * GeV); 110 pKPlusData.push_back(pKPlusPDGFit[i] * G 108 pKPlusData.push_back(pKPlusPDGFit[i] * GeV); 111 pKMinusData.push_back(pKMinusPDGFit[i] * 109 pKMinusData.push_back(pKMinusPDGFit[i] * GeV); 112 ppbarData.push_back(ppbarPDGFit[i] * GeV 110 ppbarData.push_back(ppbarPDGFit[i] * GeV); 113 npbarData.push_back(npbarPDGFit[i] * GeV 111 npbarData.push_back(npbarPDGFit[i] * GeV); 114 } 112 } 115 113 116 for (i=2; i<nPar; i++) 114 for (i=2; i<nPar; i++) 117 { 115 { 118 ppData.push_back(ppPDGFit[i]); 116 ppData.push_back(ppPDGFit[i]); 119 pPiPlusData.push_back(pPiPlusPDGFit[i]); 117 pPiPlusData.push_back(pPiPlusPDGFit[i]); 120 pPiMinusData.push_back(pPiMinusPDGFit[i] 118 pPiMinusData.push_back(pPiMinusPDGFit[i]); 121 pKPlusData.push_back(pKPlusPDGFit[i]); 119 pKPlusData.push_back(pKPlusPDGFit[i]); 122 pKMinusData.push_back(pKMinusPDGFit[i]); 120 pKMinusData.push_back(pKMinusPDGFit[i]); 123 ppbarData.push_back(ppbarPDGFit[i]); 121 ppbarData.push_back(ppbarPDGFit[i]); 124 npbarData.push_back(npbarPDGFit[i]); 122 npbarData.push_back(npbarPDGFit[i]); 125 } 123 } 126 124 127 xMap[nn] = ppData; 125 xMap[nn] = ppData; 128 xMap[pp] = ppData; 126 xMap[pp] = ppData; 129 xMap[pn] = std::move(ppData); << 127 xMap[pn] = ppData; 130 xMap[piPlusp] = std::move(pPiPlusData); << 128 xMap[piPlusp] = pPiPlusData; 131 xMap[piMinusp] = std::move(pPiMinusData); << 129 xMap[piMinusp] = pPiMinusData; 132 xMap[KPlusp] = std::move(pKPlusData); << 130 xMap[KPlusp] = pKPlusData; 133 xMap[KMinusp] = std::move(pKMinusData); << 131 xMap[KMinusp] = pKMinusData; 134 xMap[ppbar] = std::move(ppbarData); << 132 xMap[ppbar] = ppbarData; 135 xMap[npbar] = std::move(npbarData); << 133 xMap[npbar] = npbarData; 136 } 134 } 137 135 138 136 139 G4XPDGElastic::~G4XPDGElastic() 137 G4XPDGElastic::~G4XPDGElastic() 140 { } 138 { } 141 139 142 140 143 G4bool G4XPDGElastic::operator==(const G4XPDGE 141 G4bool G4XPDGElastic::operator==(const G4XPDGElastic &right) const 144 { 142 { 145 return (this == (G4XPDGElastic *) &right); 143 return (this == (G4XPDGElastic *) &right); 146 } 144 } 147 145 148 146 149 G4bool G4XPDGElastic::operator!=(const G4XPDGE 147 G4bool G4XPDGElastic::operator!=(const G4XPDGElastic &right) const 150 { 148 { 151 return (this != (G4XPDGElastic *) &right); 149 return (this != (G4XPDGElastic *) &right); 152 } 150 } 153 151 154 152 155 G4double G4XPDGElastic::CrossSection(const G4K 153 G4double G4XPDGElastic::CrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const 156 { 154 { 157 // Elastic Cross-section fit, 1994 Review of 155 // Elastic Cross-section fit, 1994 Review of Particle Properties, (1994), 1 158 156 159 G4double sigma = 0.; 157 G4double sigma = 0.; 160 158 161 const G4ParticleDefinition* def1 = trk1.GetD << 159 G4ParticleDefinition* def1 = trk1.GetDefinition(); 162 const G4ParticleDefinition* def2 = trk2.GetD << 160 G4ParticleDefinition* def2 = trk2.GetDefinition(); >> 161 G4String name1 = def1->GetParticleName(); >> 162 G4String name2 = def2->GetParticleName(); 163 163 164 G4double sqrtS = (trk1.Get4Momentum() + trk2 164 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag(); 165 G4double m_1 = def1->GetPDGMass(); << 165 G4double m1 = def1->GetPDGMass(); 166 G4double m_2 = def2->GetPDGMass(); << 166 G4double m2 = def2->GetPDGMass(); 167 G4double m_max = std::max(m_1,m_2); << 167 G4double m = std::max(m1,m2); 168 // if (m1 > m) m = m1; 168 // if (m1 > m) m = m1; 169 169 170 G4double pLab = 0.; 170 G4double pLab = 0.; 171 171 172 if (m_max > 0. && sqrtS > (m_1 + m_2)) << 172 if (m > 0. && sqrtS > (m1 + m2)) 173 { 173 { 174 pLab = std::sqrt( (sqrtS*sqrtS - (m_1+m_ << 174 pLab = std::sqrt( (sqrtS*sqrtS - (m1+m2)*(m1+m2) ) * (sqrtS*sqrtS - (m1-m2)*(m1-m2)) ) / (2*m); 175 175 176 // The PDG fit formula requires p in GeV 176 // The PDG fit formula requires p in GeV/c 177 177 >> 178 G4double enc1 = def1->GetPDGEncoding(); >> 179 G4double enc2 = def2->GetPDGEncoding(); >> 180 G4double coeff = -1.; >> 181 if ( (enc1 < 0 && enc2 >0) || (enc2 < 0 && enc1 >0) ) coeff = 1.; >> 182 178 // Order the pair: first is the lower ma 183 // Order the pair: first is the lower mass particle, second is the higher mass one 179 std::pair<const G4ParticleDefinition *,c << 184 std::pair<G4String,G4String> trkPair(def1->GetParticleName(),def2->GetParticleName()); 180 if (def1->GetPDGMass() > def2->GetPDGMas 185 if (def1->GetPDGMass() > def2->GetPDGMass()) 181 trkPair = std::pair<const G4ParticleDefiniti << 186 trkPair = std::pair<G4String,G4String>(def2->GetParticleName(),def1->GetParticleName()); 182 187 183 std::vector<G4double> data; 188 std::vector<G4double> data; 184 G4double pMinFit = 0.; 189 G4double pMinFit = 0.; 185 G4double pMaxFit = 0.; 190 G4double pMaxFit = 0.; 186 G4double aFit = 0.; 191 G4double aFit = 0.; 187 G4double bFit = 0.; 192 G4double bFit = 0.; 188 G4double cFit = 0.; 193 G4double cFit = 0.; 189 G4double dFit = 0.; 194 G4double dFit = 0.; 190 G4double nFit = 0.; 195 G4double nFit = 0.; 191 196 192 // Debug 197 // Debug 193 // G4cout << "Map has " << xMap.size() << 198 // G4cout << "Map has " << xMap.size() << " elements" << G4endl; 194 // Debug end 199 // Debug end 195 200 196 if (xMap.find(trkPair) != xMap.end()) 201 if (xMap.find(trkPair) != xMap.end()) 197 { 202 { 198 PairDoubleMap::const_iterator iter; 203 PairDoubleMap::const_iterator iter; 199 for (iter = xMap.begin(); iter != xMap.end 204 for (iter = xMap.begin(); iter != xMap.end(); ++iter) 200 { 205 { 201 std::pair<const G4ParticleDefinition * << 206 std::pair<G4String,G4String> thePair = (*iter).first; 202 if (thePair == trkPair) 207 if (thePair == trkPair) 203 { 208 { 204 data = (*iter).second; 209 data = (*iter).second; 205 pMinFit = data[0]; 210 pMinFit = data[0]; 206 pMaxFit = data[1]; 211 pMaxFit = data[1]; 207 aFit = data[2]; 212 aFit = data[2]; 208 bFit = data[3]; 213 bFit = data[3]; 209 cFit = data[5]; 214 cFit = data[5]; 210 dFit = data[6]; 215 dFit = data[6]; 211 nFit = data[4]; 216 nFit = data[4]; 212 217 213 if (pLab < pMinFit) return 0.0; 218 if (pLab < pMinFit) return 0.0; 214 if (pLab > pMaxFit ) 219 if (pLab > pMaxFit ) 215 G4cout << "WARNING! G4XPDGElastic::PDG 220 G4cout << "WARNING! G4XPDGElastic::PDGElastic " 216 << trk1.GetDefinition()->GetParticleN 221 << trk1.GetDefinition()->GetParticleName() << "-" 217 << trk2.GetDefinition()->GetParticleN 222 << trk2.GetDefinition()->GetParticleName() 218 << " elastic cross section: momentum 223 << " elastic cross section: momentum " 219 << pLab / GeV << " GeV outside valid 224 << pLab / GeV << " GeV outside valid fit range " 220 << pMinFit /GeV << "-" << pMaxFit / G 225 << pMinFit /GeV << "-" << pMaxFit / GeV 221 << G4endl; 226 << G4endl; 222 227 223 pLab /= GeV; 228 pLab /= GeV; 224 if (pLab > 0.) 229 if (pLab > 0.) 225 { 230 { 226 G4double logP = G4Log(pLab); << 231 G4double logP = std::log(pLab); 227 sigma = aFit + bFit * G4Pow::GetInst << 232 sigma = aFit + bFit * std::pow(pLab, nFit) + cFit * logP*logP + dFit * logP; 228 sigma = sigma * millibarn; 233 sigma = sigma * millibarn; 229 } 234 } 230 } 235 } 231 } 236 } 232 } 237 } 233 else 238 else 234 { 239 { 235 G4cout << "G4XPDGElastic::CrossSection - D 240 G4cout << "G4XPDGElastic::CrossSection - Data not found in Map" << G4endl; 236 } 241 } 237 } 242 } 238 243 239 if (sigma < 0.) 244 if (sigma < 0.) 240 { 245 { 241 G4cout << "WARNING! G4XPDGElastic::PDGEl 246 G4cout << "WARNING! G4XPDGElastic::PDGElastic " 242 << def1->GetParticleName() << "-" << de << 247 << name1 << "-" << name2 243 << " elastic cross section: momentum " 248 << " elastic cross section: momentum " 244 << pLab << " GeV, negative cross sectio 249 << pLab << " GeV, negative cross section " 245 << sigma / millibarn << " mb set to 0" 250 << sigma / millibarn << " mb set to 0" << G4endl; 246 sigma = 0.; 251 sigma = 0.; 247 } 252 } 248 253 249 return sigma; 254 return sigma; 250 } 255 } 251 256 252 257 253 G4String G4XPDGElastic::Name() const 258 G4String G4XPDGElastic::Name() const 254 { 259 { 255 G4String name = "PDGElastic "; 260 G4String name = "PDGElastic "; 256 return name; 261 return name; 257 } 262 } 258 263 259 264 260 G4bool G4XPDGElastic::IsValid(G4double e) cons 265 G4bool G4XPDGElastic::IsValid(G4double e) const 261 { 266 { 262 G4bool answer = InLimits(e,_lowLimit,_highLi 267 G4bool answer = InLimits(e,_lowLimit,_highLimit); 263 268 264 return answer; 269 return answer; 265 } 270 } 266 271 267 272 268 273 269 274 270 275 271 276 272 277 273 278 274 279 275 280 276 281 277 282