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 // $Id: G4AntiNuclElastic.cc - A.Galoyan 02.05.2011 >> 27 // GEANT4 tag $Name: not supported by cvs2svn $ 26 // 28 // 27 // Geant4 Header : G4AntiNuclElastic 29 // Geant4 Header : G4AntiNuclElastic 28 // 30 // 29 // 31 // 30 32 31 #include "G4AntiNuclElastic.hh" << 32 << 33 #include "G4PhysicalConstants.hh" << 34 #include "G4SystemOfUnits.hh" << 35 #include "G4ParticleTable.hh" 33 #include "G4ParticleTable.hh" 36 #include "G4ParticleDefinition.hh" 34 #include "G4ParticleDefinition.hh" 37 #include "G4IonTable.hh" 35 #include "G4IonTable.hh" 38 #include "Randomize.hh" 36 #include "Randomize.hh" 39 #include "G4AntiProton.hh" 37 #include "G4AntiProton.hh" 40 #include "G4AntiNeutron.hh" 38 #include "G4AntiNeutron.hh" 41 #include "G4AntiDeuteron.hh" 39 #include "G4AntiDeuteron.hh" 42 #include "G4AntiAlpha.hh" 40 #include "G4AntiAlpha.hh" 43 #include "G4AntiTriton.hh" 41 #include "G4AntiTriton.hh" 44 #include "G4AntiHe3.hh" 42 #include "G4AntiHe3.hh" 45 #include "G4Proton.hh" 43 #include "G4Proton.hh" 46 #include "G4Neutron.hh" 44 #include "G4Neutron.hh" 47 #include "G4Deuteron.hh" 45 #include "G4Deuteron.hh" 48 #include "G4Alpha.hh" 46 #include "G4Alpha.hh" 49 #include "G4Pow.hh" 47 #include "G4Pow.hh" 50 #include "G4Exp.hh" << 48 #include "G4AntiNuclElastic.hh" 51 #include "G4Log.hh" << 52 49 53 #include "G4NucleiProperties.hh" << 50 #include "G4NucleiProperties.hh" 54 #include "G4CrossSectionDataSetRegistry.hh" << 55 51 56 G4AntiNuclElastic::G4AntiNuclElastic() 52 G4AntiNuclElastic::G4AntiNuclElastic() 57 : G4HadronElastic("AntiAElastic") 53 : G4HadronElastic("AntiAElastic") 58 { 54 { 59 //V.Ivanchenko commented out 55 //V.Ivanchenko commented out 60 //SetMinEnergy( 0.1*GeV ); 56 //SetMinEnergy( 0.1*GeV ); 61 //SetMaxEnergy( 10.*TeV ); 57 //SetMaxEnergy( 10.*TeV ); 62 58 >> 59 63 theAProton = G4AntiProton::AntiProton( 60 theAProton = G4AntiProton::AntiProton(); 64 theANeutron = G4AntiNeutron::AntiNeutro 61 theANeutron = G4AntiNeutron::AntiNeutron(); 65 theADeuteron = G4AntiDeuteron::AntiDeute 62 theADeuteron = G4AntiDeuteron::AntiDeuteron(); 66 theATriton = G4AntiTriton::AntiTriton( 63 theATriton = G4AntiTriton::AntiTriton(); 67 theAAlpha = G4AntiAlpha::AntiAlpha(); 64 theAAlpha = G4AntiAlpha::AntiAlpha(); 68 theAHe3 = G4AntiHe3::AntiHe3(); 65 theAHe3 = G4AntiHe3::AntiHe3(); 69 66 70 theProton = G4Proton::Proton(); 67 theProton = G4Proton::Proton(); 71 theNeutron = G4Neutron::Neutron(); 68 theNeutron = G4Neutron::Neutron(); 72 theDeuteron = G4Deuteron::Deuteron(); 69 theDeuteron = G4Deuteron::Deuteron(); 73 theAlpha = G4Alpha::Alpha(); 70 theAlpha = G4Alpha::Alpha(); 74 71 75 G4CrossSectionDataSetRegistry* reg = G4Cross << 76 cs = static_cast<G4ComponentAntiNuclNuclearX << 77 if(!cs) { cs = new G4ComponentAntiNuclNuclea << 78 72 >> 73 cs = new G4ComponentAntiNuclNuclearXS(); 79 fParticle = 0; 74 fParticle = 0; 80 fWaveVector = 0.; 75 fWaveVector = 0.; 81 fBeta = 0.; 76 fBeta = 0.; 82 fZommerfeld = 0.; 77 fZommerfeld = 0.; 83 fAm = 0.; 78 fAm = 0.; 84 fTetaCMS = 0.; 79 fTetaCMS = 0.; 85 fRa = 0.; 80 fRa = 0.; 86 fRef = 0.; 81 fRef = 0.; 87 fceff = 0.; 82 fceff = 0.; 88 fptot = 0.; 83 fptot = 0.; 89 fTmax = 0.; 84 fTmax = 0.; 90 fThetaLab = 0.; 85 fThetaLab = 0.; 91 } 86 } 92 87 93 ////////////////////////////////////////////// 88 ///////////////////////////////////////////////////////////////////////// 94 G4AntiNuclElastic::~G4AntiNuclElastic() 89 G4AntiNuclElastic::~G4AntiNuclElastic() 95 {} << 90 { >> 91 delete cs; >> 92 } 96 93 97 ////////////////////////////////////////////// 94 //////////////////////////////////////////////////////////////////////// 98 // sample momentum transfer in the CMS system 95 // sample momentum transfer in the CMS system 99 G4double G4AntiNuclElastic::SampleInvariantT(c 96 G4double G4AntiNuclElastic::SampleInvariantT(const G4ParticleDefinition* particle, 100 G4double Plab, G4int Z, G4int 97 G4double Plab, G4int Z, G4int A) 101 { 98 { 102 G4double T; 99 G4double T; 103 G4double Mproj = particle->GetPDGMass(); 100 G4double Mproj = particle->GetPDGMass(); 104 G4LorentzVector Pproj(0.,0.,Plab,std::sqrt(P 101 G4LorentzVector Pproj(0.,0.,Plab,std::sqrt(Plab*Plab+Mproj*Mproj)); 105 G4double ctet1 = GetcosTeta1(Plab, A); 102 G4double ctet1 = GetcosTeta1(Plab, A); 106 103 107 G4double energy=Pproj.e()-Mproj; 104 G4double energy=Pproj.e()-Mproj; 108 105 109 const G4ParticleDefinition* theParticle = pa 106 const G4ParticleDefinition* theParticle = particle; 110 107 111 G4ParticleDefinition * theTargetDef = 0; << 108 G4ParticleDefinition * theDef = 0; 112 109 113 if (Z == 1 && A == 1) theTargetDef = th << 110 if(Z == 1 && A == 1) theDef = theProton; 114 else if (Z == 1 && A == 2) theTargetDef = th << 111 else if (Z == 1 && A == 2) theDef = theDeuteron; 115 else if (Z == 1 && A == 3) theTargetDef = G4 << 112 else if (Z == 1 && A == 3) theDef = G4Triton::Triton(); 116 else if (Z == 2 && A == 3) theTargetDef = G4 << 113 else if (Z == 2 && A == 3) theDef = G4He3::He3(); 117 else if (Z == 2 && A == 4) theTargetDef = th << 114 else if (Z == 2 && A == 4) theDef = theAlpha; 118 115 119 116 120 G4double TargMass =G4NucleiProperties::GetNu 117 G4double TargMass =G4NucleiProperties::GetNuclearMass(A,Z); 121 118 122 //transform to CMS 119 //transform to CMS 123 120 124 G4LorentzVector lv(0.0,0.0,0.0,TargMass); 121 G4LorentzVector lv(0.0,0.0,0.0,TargMass); 125 lv += Pproj; 122 lv += Pproj; 126 G4double S = lv.mag2()/(GeV*GeV); << 123 G4double S = lv.mag2()/GeV/GeV; 127 124 128 G4ThreeVector bst = lv.boostVector(); 125 G4ThreeVector bst = lv.boostVector(); 129 Pproj.boost(-bst); 126 Pproj.boost(-bst); 130 127 131 G4ThreeVector p1 = Pproj.vect(); 128 G4ThreeVector p1 = Pproj.vect(); 132 G4double ptot = p1.mag(); 129 G4double ptot = p1.mag(); 133 130 134 fbst = bst; 131 fbst = bst; 135 fptot= ptot; 132 fptot= ptot; 136 fTmax = 4.0*ptot*ptot; // In (MeV/c)^2 << 133 fTmax = 4.0*ptot*ptot; 137 134 138 if(Plab < (std::abs(particle->GetBaryonNumbe << 135 if(Plab/std::abs(particle->GetBaryonNumber()) < 100.*MeV) // Uzhi 24 Nov. 2011 139 {return fTmax*G4UniformRand();} << 136 {return fTmax*G4UniformRand();} // Uzhi 24 Nov. 2011 140 137 141 // Calculation of NN collision properties << 142 G4double PlabPerN = Plab/std::abs(theParticl << 143 G4double NucleonMass = 0.5*( theProton->GetP << 144 G4double PrNucleonMass(0.); // Projectile a << 145 if( std::abs(theParticle->GetBaryonNumber()) << 146 else << 147 G4double energyPerN = std::sqrt( sqr(PlabPer << 148 energyPerN -= PrNucleonMass; << 149 //--- << 150 << 151 G4double Z1 = particle->GetPDGCharge(); 138 G4double Z1 = particle->GetPDGCharge(); 152 G4double Z2 = Z; 139 G4double Z2 = Z; 153 140 154 G4double beta = CalculateParticleBeta(partic 141 G4double beta = CalculateParticleBeta(particle, ptot); 155 G4double n = CalculateZommerfeld( beta, Z1 142 G4double n = CalculateZommerfeld( beta, Z1, Z2 ); 156 G4double Am = CalculateAm( ptot, n, Z2 ); 143 G4double Am = CalculateAm( ptot, n, Z2 ); 157 fWaveVector = ptot; // /hbarc; 144 fWaveVector = ptot; // /hbarc; 158 145 159 G4LorentzVector Fproj(0.,0.,0.,0.); 146 G4LorentzVector Fproj(0.,0.,0.,0.); 160 const G4double mevToBarn = 0.38938e+6; << 147 G4double XsCoulomb = sqr(n/fWaveVector)*pi*(1+ctet1)/(1.+Am)/(1.+2.*Am-ctet1); 161 G4double XsCoulomb = mevToBarn*sqr(n/fWaveVe << 148 XsCoulomb=XsCoulomb*0.38938e+6; >> 149 162 150 163 G4double XsElastHadronic =cs->GetElasticElem << 151 G4double XsElastHad =cs->GetElasticElementCrossSection(particle, energy, Z, (G4double)A); 164 G4double XsTotalHadronic =cs->GetTotalElemen << 152 G4double XstotalHad =cs->GetTotalElementCrossSection(particle, energy, Z, (G4double)A); 165 153 166 XsElastHadronic/=millibarn; XsTotalHadronic/ << 167 154 168 G4double CoulombProb = XsCoulomb/(XsCoulomb << 155 XsElastHad/=millibarn; XstotalHad/=millibarn; >> 156 >> 157 >> 158 G4double CoulombProb = XsCoulomb/(XsCoulomb+XsElastHad); >> 159 >> 160 // G4cout<<" XselastHadron " << XsElastHad << " XsCol "<< XsCoulomb <<G4endl; >> 161 // G4cout <<" XsTotal" << XstotalHad <<G4endl; >> 162 // G4cout<<"XsInel"<< XstotalHad-XsElastHad<<G4endl; 169 163 170 if(G4UniformRand() < CoulombProb) 164 if(G4UniformRand() < CoulombProb) 171 { // Simulation of Coulomb scattering 165 { // Simulation of Coulomb scattering 172 166 173 G4double phi = twopi * G4UniformRand(); << 167 G4double phi = twopi * G4UniformRand(); 174 G4double Ksi = G4UniformRand(); << 168 G4double Ksi = G4UniformRand(); 175 169 176 G4double par1 = 2.*(1.+Am)/(1.+ctet1); << 170 G4double par1 = 2.*(1.+Am)/(1.+ctet1); 177 171 178 // ////sample ThetaCMS in Coulomb part << 172 // ////sample ThetaCMS in Coulomb part 179 173 180 G4double cosThetaCMS = (par1*ctet1- Ksi*(1 << 174 G4double cosThetaCMS = (par1*ctet1- Ksi*(1.+2.*Am))/(par1-Ksi); 181 175 182 G4double PtZ=ptot*cosThetaCMS; << 176 G4double PtZ=ptot*cosThetaCMS; 183 Fproj.setPz(PtZ); << 177 Fproj.setPz(PtZ); 184 G4double PtProjCMS = ptot*std::sqrt(1.0 - << 178 G4double PtProjCMS = ptot*std::sqrt(1.0 - cosThetaCMS*cosThetaCMS); 185 G4double PtX= PtProjCMS * std::cos(phi); << 179 G4double PtX= PtProjCMS * std::cos(phi); 186 G4double PtY= PtProjCMS * std::sin(phi); << 180 G4double PtY= PtProjCMS * std::sin(phi); 187 Fproj.setPx(PtX); << 181 Fproj.setPx(PtX); 188 Fproj.setPy(PtY); << 182 Fproj.setPy(PtY); 189 Fproj.setE(std::sqrt(PtX*PtX+PtY*PtY+PtZ*P << 183 Fproj.setE(std::sqrt(PtX*PtX+PtY*PtY+PtZ*PtZ+Mproj*Mproj)); 190 T = -(Pproj-Fproj).mag2(); << 184 T = -(Pproj-Fproj).mag2(); 191 } << 185 } else 192 else << 193 { << 194 // Simulation of strong interaction scatte << 195 186 196 G4double Qmax = 2.*ptot/197.33; // in fm << 187 { >> 188 ///////Simulation of strong interaction scattering//////////////////////////// 197 189 198 G4double Amag = 1.0; // A1 in Majora << 190 // G4double Qmax = 2.*ptot*197.33; // in fm^-1 199 G4double SlopeMag = 0.5; // A2 in Majora << 191 G4double Qmax = 2.*3.0*197.33; // in fm^-1 200 << 192 G4double Amag = 70*70; // A1 in Magora funct:A1*exp(-q*A2) 201 G4double sig_pbarp = cs->GetAntiHadronNucl << 193 G4double SlopeMag = 2.*3.0; // A2 in Magora funct:A1*exp(-q*A2) 202 << 194 203 fRa = 1.113*G4Pow::GetInstance()->Z13(A) - << 195 G4double sig_pbarp= cs->GetAntiHadronNucleonTotCrSc(particle,energy); 204 0.227/G4Pow::GetInstance()->Z13(A); << 196 205 if(A == 3) fRa=1.81; << 197 fRa = 1.113*G4Pow::GetInstance()->Z13(A) - 206 if(A == 4) fRa=1.37; << 198 0.227/G4Pow::GetInstance()->Z13(A); >> 199 if(A == 3) fRa=1.81; >> 200 if(A == 4) fRa=1.37; 207 201 208 if((A>=12.) && (A<27) ) fRa=fRa*0.85; << 202 if((A>=12.) && (A<27) ) fRa=fRa*0.85; 209 if((A>=27.) && (A<48) ) fRa=fRa*0.90; << 203 if((A>=27.) && (A<48) ) fRa=fRa*0.90; 210 if((A>=48.) && (A<65) ) fRa=fRa*0.95; << 204 if((A>=48.) && (A<65) ) fRa=fRa*0.95; 211 << 205 212 G4double Ref2 = XsTotalHadronic/10./2./pi; << 206 G4double Ref2 = 0; 213 G4double ceff2 = 0.0; << 207 G4double ceff2 =0; 214 G4double rho = 0.0; << 208 G4double rho = 0; 215 << 209 if ((theParticle == theAProton) || (theParticle == theANeutron)) 216 if ((theParticle == theAProton) || (thePa << 210 { 217 { << 211 if(theDef == theProton) 218 if(theTargetDef == theProton) << 212 { 219 { << 213 // G4double Mp2=sqr(theDef->GetPDGMass()/GeV ); 220 // Determination of the real part of P << 214 221 if(Plab < 610.) << 215 // change 30 October 222 { rho = 1.3347-10.342*Plab/1000.+22.27 << 216 223 13.634*Plab/1000.*Plab/1000.*P << 217 if(Plab < 610.) 224 if((Plab < 5500.)&&(Plab >= 610.) ) << 218 { rho = 1.3347-10.342*Plab/1000.+22.277*Plab/1000.*Plab/1000.- 225 { rho = 0.22; } << 219 13.634*Plab/1000.*Plab/1000.*Plab/1000. ;} 226 if((Plab >= 5500.)&&(Plab < 12300.) ) << 220 if((Plab < 5500.)&&(Plab >= 610.) ) 227 { rho = -0.32; } << 221 { rho = 0.22; } 228 if( Plab >= 12300.) << 222 if((Plab >= 5500.)&&(Plab < 12300.) ) 229 { rho = 0.135-2.26/(std::sqrt(S)) ;} << 223 { rho = -0.32; } 230 Ref2 = 0.35 + 0.9/std::sqrt(std::sqrt << 224 if( Plab >= 12300.) 231 ceff2 = 0.375 - 2./S + 0.44/(sqr(S-4.) << 225 { rho = 0.135-2.26/(std::sqrt(S)) ;} 232 Ref2 =Ref2*Ref2; << 226 233 ceff2 = ceff2*ceff2; << 227 Ref2 = 0.35 + 0.9/std::sqrt(std::sqrt(S-4.*0.88))+0.04*std::log(S) ; 234 } << 228 ceff2 = 0.375 - 2./S + 0.44/(sqr(S-4.)+1.5) ; 235 << 229 236 if( (Z==1)&&(A==2) ) << 230 /* 237 { << 231 Ref2=0.8/std::sqrt(std::sqrt(S-4.*Mp2)) + 0.55; 238 Ref2 = fRa*fRa - 0.28 + 0.019 * sig_pb << 232 if(S>1000.) Ref2=0.62+0.02*std::log(S) ; 239 ceff2 = 0.297 + 7.853e-04*sig_pbarp + << 233 ceff2 = 0.035/(sqr(S-4.3)+0.4) + 0.085 * std::log(S) ; 240 } << 234 if(S>1000.) ceff2 = 0.005 * std::log(S) + 0.29; 241 if( (Z==1)&&(A==3) ) << 235 */ 242 { << 236 243 Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pb << 237 Ref2=Ref2*Ref2; 244 ceff2 = 0.149 + 7.091e-04*sig_pbarp + << 238 ceff2 = ceff2*ceff2; 245 } << 239 246 if( (Z==2)&&(A==3) ) << 240 SlopeMag = 0.5; // Uzhi 247 { << 241 Amag= 1.; // Uzhi 248 Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pb << 242 } 249 ceff2 = 0.149 + 7.091e-04*sig_pbarp + << 243 250 } << 244 if(Z>2) 251 if( (Z==2)&&(A==4) ) << 245 { Ref2 = fRa*fRa +2.48*0.01*sig_pbarp*fRa - 2.23e-6*sig_pbarp*sig_pbarp*fRa*fRa; 252 { << 246 ceff2 = 0.16+3.3e-4*sig_pbarp+0.35*std::exp(-0.03*sig_pbarp); 253 Ref2 = fRa*fRa -0.46 +0.03*sig_pbarp - << 247 } 254 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0. << 248 if( (Z==2)&&(A==4) ) 255 } << 249 { Ref2 = fRa*fRa -0.46 +0.03*sig_pbarp - 2.98e-6*sig_pbarp*sig_pbarp; 256 if(Z>2) << 250 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*std::exp(-0.03*sig_pbarp); 257 { << 251 } 258 Ref2 = fRa*fRa +2.48*0.01*sig_pbarp*fR << 252 if( (Z==1)&&(A==3) ) 259 ceff2 = 0.16+3.3e-4*sig_pbarp+0.35*G4E << 253 { Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp; 260 } << 254 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp); 261 } // End of if ((theParticle == theAProto << 255 } 262 << 256 if( (Z==2)&&(A==3) ) 263 if (theParticle == theADeuteron) << 257 { Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp; 264 { << 258 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp); 265 if(theTargetDef == theProton) << 259 } 266 { << 260 if( (Z==1)&&(A==2) ) 267 ceff2 = 0.297 + 7.853e-04*sig_pbarp + << 261 { 268 } << 262 Ref2 = fRa*fRa - 0.28 + 0.019 * sig_pbarp + 2.06e-6 * sig_pbarp*sig_pbarp; 269 if(theTargetDef == theDeuteron) << 263 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*std::exp(-0.03*sig_pbarp); 270 { << 264 } 271 ceff2 = 0.65 + 3.0e-4*sig_pbarp + 0.55 << 265 } 272 } << 266 273 if( (theTargetDef == G4Triton::Triton()) << 267 if (theParticle == theADeuteron) 274 { << 268 { 275 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 << 269 sig_pbarp= cs->GetAntiHadronNucleonTotCrSc(particle,energy/2.); 276 } << 270 Ref2 = XstotalHad/10./2./pi ; 277 if(theTargetDef == theAlpha) << 271 if(Z>2) 278 { << 272 { 279 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.4 << 273 ceff2 = 0.38 + 2.0e-4 *sig_pbarp + 0.5 * std::exp(-0.03*sig_pbarp); 280 } << 274 } 281 if(Z>2) << 275 if(theDef == theProton) 282 { << 276 { 283 ceff2 = 0.38 + 2.0e-4 *sig_pbarp + 0.5 << 277 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*std::exp(-0.03*sig_pbarp); 284 } << 278 } 285 } << 279 if(theDef == theDeuteron) 286 << 280 { 287 if( (theParticle ==theAHe3) || (theParticl << 281 ceff2 = 0.65 + 3.0e-4*sig_pbarp + 0.55 * std::exp(-0.03*sig_pbarp); 288 { << 282 } 289 if(theTargetDef == theProton) << 283 if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) ) 290 { << 284 { 291 ceff2 = 0.149 + 7.091e-04*sig_pbarp + << 285 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * std::exp(-0.02*sig_pbarp); 292 } << 286 } 293 if(theTargetDef == theDeuteron) << 287 if(theDef == theAlpha) 294 { << 288 { 295 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 << 289 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * std::exp(-0.02*sig_pbarp); 296 } << 290 } 297 if( (theTargetDef == G4Triton::Triton()) << 291 } 298 { << 292 299 ceff2 = 0.39 + 2.7e-4*sig_pbarp + 0.7 << 293 if( (theParticle ==theAHe3) || (theParticle ==theATriton) ) 300 } << 294 { 301 if(theTargetDef == theAlpha) << 295 sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(particle,energy/3.); 302 { << 296 Ref2 = XstotalHad/10./2./pi ; 303 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 << 297 if(Z>2) 304 } << 298 { 305 if(Z>2) << 299 ceff2 = 0.26 + 2.2e-4*sig_pbarp + 0.33*std::exp(-0.03*sig_pbarp); 306 { << 300 } 307 ceff2 = 0.26 + 2.2e-4*sig_pbarp + 0.33 << 301 if(theDef == theProton) 308 } << 302 { 309 } << 303 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp); 310 << 304 } 311 if ( (theParticle == theAAlpha) || (ceff2 << 305 if(theDef == theDeuteron) 312 { << 306 { 313 if(theTargetDef == theProton) << 307 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * std::exp(-0.02*sig_pbarp); 314 { << 308 } 315 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0. << 309 if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) ) 316 } << 310 { 317 if(theTargetDef == theDeuteron) << 311 ceff2 = 0.39 + 2.7e-4*sig_pbarp + 0.7 * std::exp(-0.02*sig_pbarp); 318 { << 312 } 319 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.4 << 313 if(theDef == theAlpha) 320 } << 314 { 321 if( (theTargetDef == G4Triton::Triton()) << 315 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * std::exp(-0.03*sig_pbarp); 322 { << 316 } 323 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 << 317 } 324 } << 318 325 if(theTargetDef == theAlpha) << 319 326 { << 320 if (theParticle == theAAlpha) 327 ceff2 = 0.17 + 3.5e-4*sig_pbarp + 0.45 << 321 { 328 } << 322 sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(particle,energy/3.); 329 if(Z>2) << 323 Ref2 = XstotalHad/10./2./pi ; 330 { << 324 if(Z>2) 331 ceff2 = 0.22 + 2.0e-4*sig_pbarp + 0.2 << 325 { 332 } << 326 ceff2 = 0.22 + 2.0e-4*sig_pbarp + 0.2 * std::exp(-0.03*sig_pbarp); 333 } << 327 } 334 << 328 if(theDef == theProton) 335 fRef=std::sqrt(Ref2); << 329 { 336 fceff = std::sqrt(ceff2); << 330 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*std::exp(-0.03*sig_pbarp); 337 << 338 G4double Q = 0.0 ; << 339 G4double BracFunct; << 340 << 341 const G4int maxNumberOfLoops = 10000; << 342 G4int loopCounter = 0; << 343 do << 344 { << 345 Q = -G4Log(1.-(1.- G4Exp(-SlopeMag * Qma << 346 G4double x = fRef * Q; << 347 BracFunct = ( ( sqr(BesselOneByArg(x))+s << 348 * sqr(DampFactor(pi*fceff*Q))) /(Ama << 349 BracFunct = BracFunct * Q; << 350 } << 351 while ( (G4UniformRand()>BracFunct) && << 352 ++loopCounter < maxNumberOfLoops ) << 353 if ( loopCounter >= maxNumberOfLoops ) { << 354 fTetaCMS = 0.0; << 355 return 0.0; << 356 } 331 } >> 332 if(theDef == theDeuteron) >> 333 { >> 334 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * std::exp(-0.02*sig_pbarp); >> 335 } >> 336 if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) ) >> 337 { >> 338 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * std::exp(-0.03*sig_pbarp); >> 339 } >> 340 if(theDef == theAlpha) >> 341 { >> 342 ceff2 = 0.17 + 3.5e-4*sig_pbarp + 0.45 * std::exp(-0.03*sig_pbarp); >> 343 } >> 344 } >> 345 >> 346 fRef=std::sqrt(Ref2); >> 347 fceff = std::sqrt(ceff2); >> 348 // G4cout<<" Ref "<<fRef<<" c_eff "<<fceff<< " rho "<< rho<<G4endl; 357 349 358 T= sqr(Q); << 350 359 T*=3.893913e+4; // fm^(-2) -> MeV^2 << 351 G4double Q = 0.0 ; >> 352 G4double BracFunct; >> 353 do >> 354 { >> 355 Q = -std::log(1.-(1.- std::exp(-SlopeMag * Qmax))* G4UniformRand() )/SlopeMag; >> 356 G4double x = fRef * Q; >> 357 BracFunct = ( ( sqr(BesselOneByArg(x))+sqr(rho/2. * BesselJzero(x)) ) >> 358 * sqr(DampFactor(pi*fceff*Q))) /(Amag*std::exp(-SlopeMag*Q)); >> 359 >> 360 BracFunct = BracFunct * Q * sqr(sqr(fRef)); >> 361 } >> 362 while (G4UniformRand()>BracFunct); >> 363 >> 364 T= sqr(Q); >> 365 T*=3.893913e+4; // fm -> MeV^2 >> 366 } 360 367 361 } // End of simulation of strong interactio << 368 G4double cosTet=1.0-T/(2.*ptot*ptot); >> 369 >> 370 fTetaCMS=std::acos(cosTet); 362 371 363 return T; << 372 return T; 364 } 373 } 365 374 366 ////////////////////////////////////////////// 375 ///////////////////////////////////////////////////////////////////// 367 // Sample of Theta in CMS 376 // Sample of Theta in CMS 368 G4double G4AntiNuclElastic::SampleThetaCMS(co 377 G4double G4AntiNuclElastic::SampleThetaCMS(const G4ParticleDefinition* p, G4double plab, 369 378 G4int Z, G4int A) 370 { 379 { 371 G4double T; 380 G4double T; 372 T = SampleInvariantT( p, plab, Z, A); 381 T = SampleInvariantT( p, plab, Z, A); 373 382 374 // NaN finder 383 // NaN finder 375 if(!(T < 0.0 || T >= 0.0)) 384 if(!(T < 0.0 || T >= 0.0)) 376 { 385 { 377 if (verboseLevel > 0) 386 if (verboseLevel > 0) 378 { 387 { 379 G4cout << "G4DiffuseElastic:WARNING: A = 388 G4cout << "G4DiffuseElastic:WARNING: A = " << A 380 << " mom(GeV)= " << plab/GeV 389 << " mom(GeV)= " << plab/GeV 381 << " S-wave will be sampled" 390 << " S-wave will be sampled" 382 << G4endl; 391 << G4endl; 383 } 392 } 384 T = G4UniformRand()*fTmax; 393 T = G4UniformRand()*fTmax; 385 394 386 } 395 } 387 396 388 if(fptot > 0.) << 397 if(fptot > 0.) // Uzhi 24 Nov. 2011 389 { 398 { 390 G4double cosTet=1.0-T/(2.*fptot*fptot); 399 G4double cosTet=1.0-T/(2.*fptot*fptot); 391 if(cosTet > 1.0 ) cosTet= 1.; << 392 if(cosTet < -1.0 ) cosTet=-1.; << 393 fTetaCMS=std::acos(cosTet); 400 fTetaCMS=std::acos(cosTet); 394 return fTetaCMS; 401 return fTetaCMS; 395 } else << 402 } else // Uzhi 24 Nov. 2011 396 { << 403 { // Uzhi 24 Nov. 2011 397 return 2.*G4UniformRand()-1.; << 404 return 2.*G4UniformRand()-1.; // Uzhi 24 Nov. 2011 398 } << 405 } // Uzhi 24 Nov. 2011 399 } 406 } 400 407 401 408 402 ////////////////////////////////////////////// 409 ///////////////////////////////////////////////////////////////////// 403 // Sample of Theta in Lab System 410 // Sample of Theta in Lab System 404 G4double G4AntiNuclElastic::SampleThetaLab(co 411 G4double G4AntiNuclElastic::SampleThetaLab(const G4ParticleDefinition* p, G4double plab, 405 412 G4int Z, G4int A) 406 { 413 { 407 G4double T; 414 G4double T; 408 T = SampleInvariantT( p, plab, Z, A); 415 T = SampleInvariantT( p, plab, Z, A); 409 416 410 // NaN finder 417 // NaN finder 411 if(!(T < 0.0 || T >= 0.0)) 418 if(!(T < 0.0 || T >= 0.0)) 412 { 419 { 413 if (verboseLevel > 0) 420 if (verboseLevel > 0) 414 { 421 { 415 G4cout << "G4DiffuseElastic:WARNING: A = 422 G4cout << "G4DiffuseElastic:WARNING: A = " << A 416 << " mom(GeV)= " << plab/GeV 423 << " mom(GeV)= " << plab/GeV 417 << " S-wave will be sampled" 424 << " S-wave will be sampled" 418 << G4endl; 425 << G4endl; 419 } 426 } 420 T = G4UniformRand()*fTmax; 427 T = G4UniformRand()*fTmax; 421 } 428 } 422 429 423 G4double phi = G4UniformRand()*twopi; 430 G4double phi = G4UniformRand()*twopi; 424 431 425 G4double cost(1.); 432 G4double cost(1.); 426 if(fTmax > 0.) {cost = 1. - 2.0*T/fTmax;} << 433 if(fTmax > 0.) {cost = 1. - 2.0*T/fTmax;} // Uzhi 24 Nov. 2011 427 434 428 G4double sint; 435 G4double sint; 429 if( cost >= 1.0 ) 436 if( cost >= 1.0 ) 430 { 437 { 431 cost = 1.0; 438 cost = 1.0; 432 sint = 0.0; 439 sint = 0.0; 433 } 440 } 434 else if( cost <= -1.0) 441 else if( cost <= -1.0) 435 { 442 { 436 cost = -1.0; 443 cost = -1.0; 437 sint = 0.0; 444 sint = 0.0; 438 } 445 } 439 else 446 else 440 { 447 { 441 sint = std::sqrt((1.0-cost)*(1.0+cost)); 448 sint = std::sqrt((1.0-cost)*(1.0+cost)); 442 } 449 } 443 450 444 G4double m1 = p->GetPDGMass(); 451 G4double m1 = p->GetPDGMass(); 445 G4ThreeVector v(sint*std::cos(phi),sint*std: 452 G4ThreeVector v(sint*std::cos(phi),sint*std::sin(phi),cost); 446 v *= fptot; 453 v *= fptot; 447 G4LorentzVector nlv(v.x(),v.y(),v.z(),std::s 454 G4LorentzVector nlv(v.x(),v.y(),v.z(),std::sqrt(fptot*fptot + m1*m1)); 448 455 449 nlv.boost(fbst); 456 nlv.boost(fbst); 450 457 451 G4ThreeVector np = nlv.vect(); 458 G4ThreeVector np = nlv.vect(); 452 G4double theta = np.theta(); 459 G4double theta = np.theta(); 453 fThetaLab = theta; 460 fThetaLab = theta; 454 461 455 return theta; 462 return theta; 456 } 463 } 457 464 458 ////////////////////////////////////////////// 465 //////////////////////////////////////////////////////////////////// 459 // Calculation of Damp factor 466 // Calculation of Damp factor 460 G4double G4AntiNuclElastic::DampFactor(G4doub 467 G4double G4AntiNuclElastic::DampFactor(G4double x) 461 { 468 { 462 G4double df; 469 G4double df; 463 G4double f3 = 6.; // first factorials 470 G4double f3 = 6.; // first factorials 464 471 465 if( std::fabs(x) < 0.01 ) 472 if( std::fabs(x) < 0.01 ) 466 { 473 { 467 df=1./(1.+x*x/f3); 474 df=1./(1.+x*x/f3); 468 } 475 } 469 else 476 else 470 { 477 { 471 df = x/std::sinh(x); 478 df = x/std::sinh(x); 472 } 479 } 473 return df; 480 return df; 474 } 481 } 475 482 476 483 477 ////////////////////////////////////////////// 484 ///////////////////////////////////////////////////////////////////////////////// 478 // Calculation of particle velocity Beta 485 // Calculation of particle velocity Beta 479 486 480 G4double G4AntiNuclElastic::CalculateParticle 487 G4double G4AntiNuclElastic::CalculateParticleBeta( const G4ParticleDefinition* particle, 481 G4double mom 488 G4double momentum ) 482 { 489 { 483 G4double mass = particle->GetPDGMass(); 490 G4double mass = particle->GetPDGMass(); 484 G4double a = momentum/mass; 491 G4double a = momentum/mass; 485 fBeta = a/std::sqrt(1+a*a); 492 fBeta = a/std::sqrt(1+a*a); 486 493 487 return fBeta; 494 return fBeta; 488 } 495 } 489 496 490 497 491 ////////////////////////////////////////////// 498 /////////////////////////////////////////////////////////////////////////////////// 492 // Calculation of parameter Zommerfeld 499 // Calculation of parameter Zommerfeld 493 500 494 G4double G4AntiNuclElastic::CalculateZommerfe 501 G4double G4AntiNuclElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 ) 495 { 502 { 496 fZommerfeld = fine_structure_const*Z1*Z2/bet 503 fZommerfeld = fine_structure_const*Z1*Z2/beta; 497 504 498 return fZommerfeld; 505 return fZommerfeld; 499 } 506 } 500 507 501 ////////////////////////////////////////////// 508 //////////////////////////////////////////////////////////////////////////////////// 502 // 509 // 503 G4double G4AntiNuclElastic::CalculateAm( G4dou 510 G4double G4AntiNuclElastic::CalculateAm( G4double momentum, G4double n, G4double Z) 504 { 511 { 505 G4double k = momentum/hbarc; 512 G4double k = momentum/hbarc; 506 G4double ch = 1.13 + 3.76*n*n; 513 G4double ch = 1.13 + 3.76*n*n; 507 G4double zn = 1.77*k/G4Pow::GetInstance()-> 514 G4double zn = 1.77*k/G4Pow::GetInstance()->A13(Z)*Bohr_radius; 508 G4double zn2 = zn*zn; 515 G4double zn2 = zn*zn; 509 fAm = ch/zn2; 516 fAm = ch/zn2; 510 517 511 return fAm; 518 return fAm; 512 } 519 } 513 520 514 ////////////////////////////////////////////// 521 ///////////////////////////////////////////////////////////// 515 // 522 // 516 // Bessel J0 function based on rational approx 523 // Bessel J0 function based on rational approximation from 517 // J.F. Hart, Computer Approximations, New Yor 524 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 518 525 519 G4double G4AntiNuclElastic::BesselJzero(G4doub 526 G4double G4AntiNuclElastic::BesselJzero(G4double value) 520 { 527 { 521 G4double modvalue, value2, fact1, fact2, arg 528 G4double modvalue, value2, fact1, fact2, arg, shift, bessel; 522 529 523 modvalue = std::fabs(value); 530 modvalue = std::fabs(value); 524 531 525 if ( value < 8.0 && value > -8.0 ) 532 if ( value < 8.0 && value > -8.0 ) 526 { 533 { 527 value2 = value*value; 534 value2 = value*value; 528 535 529 fact1 = 57568490574.0 + value2*(-13362590 536 fact1 = 57568490574.0 + value2*(-13362590354.0 530 + value2*( 65161964 537 + value2*( 651619640.7 531 + value2*(-11214424 538 + value2*(-11214424.18 532 + value2*( 77392.33 539 + value2*( 77392.33017 533 + value2*(-184.9052 540 + value2*(-184.9052456 ) ) ) ) ); 534 541 535 fact2 = 57568490411.0 + value2*( 10295329 542 fact2 = 57568490411.0 + value2*( 1029532985.0 536 + value2*( 9494680. 543 + value2*( 9494680.718 537 + value2*(59272.648 544 + value2*(59272.64853 538 + value2*(267.85327 545 + value2*(267.8532712 539 + value2*1.0 546 + value2*1.0 ) ) ) ); 540 547 541 bessel = fact1/fact2; 548 bessel = fact1/fact2; 542 } 549 } 543 else 550 else 544 { 551 { 545 arg = 8.0/modvalue; 552 arg = 8.0/modvalue; 546 553 547 value2 = arg*arg; 554 value2 = arg*arg; 548 555 549 shift = modvalue-0.785398164; 556 shift = modvalue-0.785398164; 550 557 551 fact1 = 1.0 + value2*(-0.1098628627e-2 558 fact1 = 1.0 + value2*(-0.1098628627e-2 552 + value2*(0.2734510407e-4 559 + value2*(0.2734510407e-4 553 + value2*(-0.2073370639e-5 560 + value2*(-0.2073370639e-5 554 + value2*0.2093887211e-6 ) 561 + value2*0.2093887211e-6 ) ) ); 555 fact2 = -0.1562499995e-1 + value2*(0.143048 562 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3 556 + value2*(-0.691 563 + value2*(-0.6911147651e-5 557 + value2*(0.7621 564 + value2*(0.7621095161e-6 558 - value2*0.93494 565 - value2*0.934945152e-7 ) ) ); 559 566 560 bessel = std::sqrt(0.636619772/modvalue)*( 567 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2); 561 } 568 } 562 return bessel; 569 return bessel; 563 } 570 } 564 571 565 572 566 ////////////////////////////////////////////// 573 ////////////////////////////////////////////////////////////////////////////// 567 // Bessel J1 function based on rational approx 574 // Bessel J1 function based on rational approximation from 568 // J.F. Hart, Computer Approximations, New Yor 575 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 569 576 570 G4double G4AntiNuclElastic::BesselJone(G4doub 577 G4double G4AntiNuclElastic::BesselJone(G4double value) 571 { 578 { 572 G4double modvalue, value2, fact1, fact2, arg 579 G4double modvalue, value2, fact1, fact2, arg, shift, bessel; 573 580 574 modvalue = std::fabs(value); 581 modvalue = std::fabs(value); 575 582 576 if ( modvalue < 8.0 ) 583 if ( modvalue < 8.0 ) 577 { 584 { 578 value2 = value*value; 585 value2 = value*value; 579 fact1 = value*(72362614232.0 + value2*(-7 586 fact1 = value*(72362614232.0 + value2*(-7895059235.0 580 + value2*( 2 587 + value2*( 242396853.1 581 + value2*(-2 588 + value2*(-2972611.439 582 + value2*( 1 589 + value2*( 15704.48260 583 + value2*(-3 590 + value2*(-30.16036606 ) ) ) ) ) ); 584 591 585 fact2 = 144725228442.0 + value2*(23005351 592 fact2 = 144725228442.0 + value2*(2300535178.0 586 + value2*(18583304 593 + value2*(18583304.74 587 + value2*(99447.43 594 + value2*(99447.43394 588 + value2*(376.9991 595 + value2*(376.9991397 589 + value2*1.0 596 + value2*1.0 ) ) ) ); 590 bessel = fact1/fact2; 597 bessel = fact1/fact2; 591 } 598 } 592 else 599 else 593 { 600 { 594 arg = 8.0/modvalue; 601 arg = 8.0/modvalue; 595 value2 = arg*arg; 602 value2 = arg*arg; 596 603 597 shift = modvalue - 2.356194491; 604 shift = modvalue - 2.356194491; 598 605 599 fact1 = 1.0 + value2*( 0.183105e-2 606 fact1 = 1.0 + value2*( 0.183105e-2 600 + value2*(-0.3516396496e-4 607 + value2*(-0.3516396496e-4 601 + value2*(0.2457520174e-5 608 + value2*(0.2457520174e-5 602 + value2*(-0.240337019e-6 609 + value2*(-0.240337019e-6 ) ) ) ); 603 610 604 fact2 = 0.04687499995 + value2*(-0.2002690 611 fact2 = 0.04687499995 + value2*(-0.2002690873e-3 605 + value2*( 0.8449199 612 + value2*( 0.8449199096e-5 606 + value2*(-0.8822898 613 + value2*(-0.88228987e-6 607 + value2*0.105787412 614 + value2*0.105787412e-6 ) ) ); 608 615 609 bessel = std::sqrt( 0.636619772/modvalue)* 616 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2); 610 if (value < 0.0) bessel = -bessel; 617 if (value < 0.0) bessel = -bessel; 611 } 618 } 612 return bessel; 619 return bessel; 613 } 620 } 614 621 615 ////////////////////////////////////////////// 622 //////////////////////////////////////////////////////////////////////////////// 616 // return J1(x)/x with special case for small 623 // return J1(x)/x with special case for small x 617 G4double G4AntiNuclElastic::BesselOneByArg(G4d 624 G4double G4AntiNuclElastic::BesselOneByArg(G4double x) 618 { 625 { 619 G4double x2, result; 626 G4double x2, result; 620 627 621 if( std::fabs(x) < 0.01 ) 628 if( std::fabs(x) < 0.01 ) 622 { 629 { 623 x *= 0.5; 630 x *= 0.5; 624 x2 = x*x; 631 x2 = x*x; 625 result = (2.- x2 + x2*x2/6.)/4.; 632 result = (2.- x2 + x2*x2/6.)/4.; 626 } 633 } 627 else 634 else 628 { 635 { 629 result = BesselJone(x)/x; 636 result = BesselJone(x)/x; 630 } 637 } 631 return result; 638 return result; 632 } 639 } 633 640 634 ////////////////////////////////////////////// 641 ///////////////////////////////////////////////////////////////////////////////// 635 // return angle from which Coulomb scattering 642 // return angle from which Coulomb scattering is calculated 636 G4double G4AntiNuclElastic::GetcosTeta1(G4doub 643 G4double G4AntiNuclElastic::GetcosTeta1(G4double plab, G4int A) 637 { 644 { 638 645 639 // G4double p0 =G4LossTableManager::Instance() 646 // G4double p0 =G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi; 640 G4double p0 = 1.*hbarc/fermi; 647 G4double p0 = 1.*hbarc/fermi; 641 //G4double cteta1 = 1.0 - p0*p0/2.0 * pow(A,2. 648 //G4double cteta1 = 1.0 - p0*p0/2.0 * pow(A,2./3.)/(plab*plab); 642 G4double cteta1 = 1.0 - p0*p0/2.0 * G4Pow::G 649 G4double cteta1 = 1.0 - p0*p0/2.0 * G4Pow::GetInstance()->Z23(A)/(plab*plab); 643 ////////////////// 650 ////////////////// 644 if(cteta1 < -1.) cteta1 = -1.0; 651 if(cteta1 < -1.) cteta1 = -1.0; 645 return cteta1; 652 return cteta1; 646 } 653 } 647 654 648 655 649 656 650 657 651 658 652 659 653 660 654 661