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