00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041 #include "G4ChipsAntiBaryonElasticXS.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4DynamicParticle.hh"
00044 #include "G4ParticleDefinition.hh"
00045 #include "G4AntiProton.hh"
00046 #include "G4Nucleus.hh"
00047 #include "G4ParticleTable.hh"
00048 #include "G4NucleiProperties.hh"
00049
00050
00051 #include "G4CrossSectionFactory.hh"
00052
00053 G4_DECLARE_XS_FACTORY(G4ChipsAntiBaryonElasticXS);
00054
00055 G4ChipsAntiBaryonElasticXS::G4ChipsAntiBaryonElasticXS():G4VCrossSectionDataSet(Default_Name()), nPoints(128), nLast(nPoints-1)
00056 {
00057 lPMin=-8.;
00058 lPMax= 8.;
00059 dlnP=(lPMax-lPMin)/nLast;
00060 onlyCS=true;
00061 lastSIG=0.;
00062 lastLP=-10.;
00063 lastTM=0.;
00064 theSS=0.;
00065 theS1=0.;
00066 theB1=0.;
00067 theS2=0.;
00068 theB2=0.;
00069 theS3=0.;
00070 theB3=0.;
00071 theS4=0.;
00072 theB4=0.;
00073 lastTZ=0;
00074 lastTN=0;
00075 lastPIN=0.;
00076 lastCST=0;
00077 lastPAR=0;
00078 lastSST=0;
00079 lastS1T=0;
00080 lastB1T=0;
00081 lastS2T=0;
00082 lastB2T=0;
00083 lastS3T=0;
00084 lastB3T=0;
00085 lastS4T=0;
00086 lastB4T=0;
00087 lastN=0;
00088 lastZ=0;
00089 lastP=0.;
00090 lastTH=0.;
00091 lastCS=0.;
00092 lastI=0;
00093 }
00094
00095 G4ChipsAntiBaryonElasticXS::~G4ChipsAntiBaryonElasticXS()
00096 {
00097 std::vector<G4double*>::iterator pos;
00098 for (pos=CST.begin(); pos<CST.end(); pos++)
00099 { delete [] *pos; }
00100 CST.clear();
00101 for (pos=PAR.begin(); pos<PAR.end(); pos++)
00102 { delete [] *pos; }
00103 PAR.clear();
00104 for (pos=SST.begin(); pos<SST.end(); pos++)
00105 { delete [] *pos; }
00106 SST.clear();
00107 for (pos=S1T.begin(); pos<S1T.end(); pos++)
00108 { delete [] *pos; }
00109 S1T.clear();
00110 for (pos=B1T.begin(); pos<B1T.end(); pos++)
00111 { delete [] *pos; }
00112 B1T.clear();
00113 for (pos=S2T.begin(); pos<S2T.end(); pos++)
00114 { delete [] *pos; }
00115 S2T.clear();
00116 for (pos=B2T.begin(); pos<B2T.end(); pos++)
00117 { delete [] *pos; }
00118 B2T.clear();
00119 for (pos=S3T.begin(); pos<S3T.end(); pos++)
00120 { delete [] *pos; }
00121 S3T.clear();
00122 for (pos=B3T.begin(); pos<B3T.end(); pos++)
00123 { delete [] *pos; }
00124 B3T.clear();
00125 for (pos=S4T.begin(); pos<S4T.end(); pos++)
00126 { delete [] *pos; }
00127 S4T.clear();
00128 for (pos=B4T.begin(); pos<B4T.end(); pos++)
00129 { delete [] *pos; }
00130 B4T.clear();
00131 }
00132
00133
00134 G4bool G4ChipsAntiBaryonElasticXS::IsIsoApplicable(const G4DynamicParticle* Pt, G4int, G4int,
00135 const G4Element*,
00136 const G4Material*)
00137 {
00138 G4ParticleDefinition* particle = Pt->GetDefinition();
00139
00140 if(particle == G4AntiNeutron::AntiNeutron())
00141 {
00142 return true;
00143 }
00144 else if(particle == G4AntiProton::AntiProton())
00145 {
00146 return true;
00147 }
00148 else if(particle == G4AntiLambda::AntiLambda())
00149 {
00150 return true;
00151 }
00152 else if(particle == G4AntiSigmaPlus::AntiSigmaPlus())
00153 {
00154 return true;
00155 }
00156 else if(particle == G4AntiSigmaMinus::AntiSigmaMinus())
00157 {
00158 return true;
00159 }
00160 else if(particle == G4AntiSigmaZero::AntiSigmaZero())
00161 {
00162 return true;
00163 }
00164 else if(particle == G4AntiXiMinus::AntiXiMinus())
00165 {
00166 return true;
00167 }
00168 else if(particle == G4AntiXiZero::AntiXiZero())
00169 {
00170 return true;
00171 }
00172 else if(particle == G4AntiOmegaMinus::AntiOmegaMinus())
00173 {
00174 return true;
00175 }
00176 return false;
00177 }
00178
00179
00180
00181 G4double G4ChipsAntiBaryonElasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,
00182 const G4Isotope*,
00183 const G4Element*,
00184 const G4Material*)
00185 {
00186 G4double pMom=Pt->GetTotalMomentum();
00187 G4int tgN = A - tgZ;
00188 G4int pdg = Pt->GetDefinition()->GetPDGEncoding();
00189
00190 return GetChipsCrossSection(pMom, tgZ, tgN, pdg);
00191 }
00192
00193 G4double G4ChipsAntiBaryonElasticXS::GetChipsCrossSection(G4double pMom, G4int tgZ, G4int tgN, G4int pPDG)
00194 {
00195 static std::vector <G4int> colN;
00196 static std::vector <G4int> colZ;
00197 static std::vector <G4double> colP;
00198 static std::vector <G4double> colTH;
00199 static std::vector <G4double> colCS;
00200
00201
00202 G4bool fCS = false;
00203
00204 G4double pEn=pMom;
00205 onlyCS=fCS;
00206
00207 G4bool in=false;
00208 lastP = 0.;
00209 lastN = tgN;
00210 lastZ = tgZ;
00211 lastI = colN.size();
00212 if(lastI) for(G4int i=0; i<lastI; i++)
00213 {
00214 if(colN[i]==tgN && colZ[i]==tgZ)
00215 {
00216 lastI=i;
00217 lastTH =colTH[i];
00218 if(pEn<=lastTH)
00219 {
00220 return 0.;
00221 }
00222 lastP =colP [i];
00223 lastCS =colCS[i];
00224
00225 if(lastP == pMom)
00226 {
00227 CalculateCrossSection(fCS,-1,i,pPDG,lastZ,lastN,pMom);
00228 return lastCS*millibarn;
00229 }
00230 in = true;
00231
00232 lastCS=CalculateCrossSection(fCS,-1,i,pPDG,lastZ,lastN,pMom);
00233 if(lastCS<=0. && pEn>lastTH)
00234 {
00235 lastTH=pEn;
00236 }
00237 break;
00238 }
00239 }
00240 if(!in)
00241 {
00243 lastCS=CalculateCrossSection(fCS,0,lastI,pPDG,lastZ,lastN,pMom);
00244 if(lastCS<=0.)
00245 {
00246 lastTH = 0;
00247 if(pEn>lastTH)
00248 {
00249 lastTH=pEn;
00250 }
00251 }
00252 colN.push_back(tgN);
00253 colZ.push_back(tgZ);
00254 colP.push_back(pMom);
00255 colTH.push_back(lastTH);
00256 colCS.push_back(lastCS);
00257 return lastCS*millibarn;
00258 }
00259 else
00260 {
00261 colP[lastI]=pMom;
00262 colCS[lastI]=lastCS;
00263 }
00264 return lastCS*millibarn;
00265 }
00266
00267
00268
00269 G4double G4ChipsAntiBaryonElasticXS::CalculateCrossSection(G4bool CS,G4int F,G4int I,
00270 G4int PDG, G4int tgZ, G4int tgN, G4double pIU)
00271 {
00272
00273 static std::vector <G4double> PIN;
00274
00275 G4double pMom=pIU/GeV;
00276 onlyCS=CS;
00277 lastLP=std::log(pMom);
00278 if(F)
00279 {
00280 if(F<0)
00281 {
00282 lastPIN = PIN[I];
00283 lastPAR = PAR[I];
00284 lastCST = CST[I];
00285 lastSST = SST[I];
00286 lastS1T = S1T[I];
00287 lastB1T = B1T[I];
00288 lastS2T = S2T[I];
00289 lastB2T = B2T[I];
00290 lastS3T = S3T[I];
00291 lastB3T = B3T[I];
00292 lastS4T = S4T[I];
00293 lastB4T = B4T[I];
00294 }
00295 if(lastLP>lastPIN && lastLP<lPMax)
00296 {
00297 lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
00298 PIN[I]=lastPIN;
00299 }
00300 }
00301 else
00302 {
00303 lastPAR = new G4double[nPoints];
00304 lastPAR[nLast]=0;
00305 lastCST = new G4double[nPoints];
00306 lastSST = new G4double[nPoints];
00307 lastS1T = new G4double[nPoints];
00308 lastB1T = new G4double[nPoints];
00309 lastS2T = new G4double[nPoints];
00310 lastB2T = new G4double[nPoints];
00311 lastS3T = new G4double[nPoints];
00312 lastB3T = new G4double[nPoints];
00313 lastS4T = new G4double[nPoints];
00314 lastB4T = new G4double[nPoints];
00315 lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN);
00316 PIN.push_back(lastPIN);
00317 PAR.push_back(lastPAR);
00318 CST.push_back(lastCST);
00319 SST.push_back(lastSST);
00320 S1T.push_back(lastS1T);
00321 B1T.push_back(lastB1T);
00322 S2T.push_back(lastS2T);
00323 B2T.push_back(lastB2T);
00324 S3T.push_back(lastS3T);
00325 B3T.push_back(lastB3T);
00326 S4T.push_back(lastS4T);
00327 B4T.push_back(lastB4T);
00328 }
00329
00330 if(lastLP>lastPIN && lastLP<lPMax)
00331 {
00332 lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
00333 }
00334 if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom);
00335 if(lastLP>lPMin && lastLP<=lastPIN)
00336 {
00337 if(lastLP==lastPIN)
00338 {
00339 G4double shift=(lastLP-lPMin)/dlnP+.000001;
00340 G4int blast=static_cast<int>(shift);
00341 if(blast<0 || blast>=nLast) G4cout<<"G4QaBarElCS::CCS:b="<<blast<<","<<nLast<<G4endl;
00342 lastSIG = lastCST[blast];
00343 if(!onlyCS)
00344 {
00345 theSS = lastSST[blast];
00346 theS1 = lastS1T[blast];
00347 theB1 = lastB1T[blast];
00348 theS2 = lastS2T[blast];
00349 theB2 = lastB2T[blast];
00350 theS3 = lastS3T[blast];
00351 theB3 = lastB3T[blast];
00352 theS4 = lastS4T[blast];
00353 theB4 = lastB4T[blast];
00354 }
00355 }
00356 else
00357 {
00358 G4double shift=(lastLP-lPMin)/dlnP;
00359 G4int blast=static_cast<int>(shift);
00360 if(blast<0) blast=0;
00361 if(blast>=nLast) blast=nLast-1;
00362 shift-=blast;
00363 G4int lastL=blast+1;
00364 G4double SIGL=lastCST[blast];
00365 lastSIG= SIGL+shift*(lastCST[lastL]-SIGL);
00366 if(!onlyCS)
00367 {
00368 G4double SSTL=lastSST[blast];
00369 theSS=SSTL+shift*(lastSST[lastL]-SSTL);
00370 G4double S1TL=lastS1T[blast];
00371 theS1=S1TL+shift*(lastS1T[lastL]-S1TL);
00372 G4double B1TL=lastB1T[blast];
00373 theB1=B1TL+shift*(lastB1T[lastL]-B1TL);
00374 G4double S2TL=lastS2T[blast];
00375 theS2=S2TL+shift*(lastS2T[lastL]-S2TL);
00376 G4double B2TL=lastB2T[blast];
00377 theB2=B2TL+shift*(lastB2T[lastL]-B2TL);
00378 G4double S3TL=lastS3T[blast];
00379 theS3=S3TL+shift*(lastS3T[lastL]-S3TL);
00380 G4double B3TL=lastB3T[blast];
00381 theB3=B3TL+shift*(lastB3T[lastL]-B3TL);
00382 G4double S4TL=lastS4T[blast];
00383 theS4=S4TL+shift*(lastS4T[lastL]-S4TL);
00384 G4double B4TL=lastB4T[blast];
00385 theB4=B4TL+shift*(lastB4T[lastL]-B4TL);
00386 }
00387 }
00388 }
00389 else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN);
00390 if(lastSIG<0.) lastSIG = 0.;
00391 return lastSIG;
00392 }
00393
00394
00395 G4double G4ChipsAntiBaryonElasticXS::GetPTables(G4double LP, G4double ILP, G4int PDG,
00396 G4int tgZ, G4int tgN)
00397 {
00398
00399 static const G4double pwd=2727;
00400 const G4int n_appel=30;
00401
00402 G4double app_el[n_appel]={1.25,3.5,80.,1.,.0557,6.72,5.,74.,3.,3.4,.2,.17,.001,8.,.055,
00403 3.64,5.e-5,4000.,1500.,.46,1.2e6,3.5e6,5.e-5,1.e10,8.5e8,
00404 1.e10,1.1,3.4e6,6.8e6,0.};
00405
00406
00407 if(PDG>-3334 && PDG<-1111)
00408 {
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420 if(lastPAR[nLast]!=pwd)
00421 {
00422 if ( tgZ == 1 && tgN == 0 )
00423 {
00424 for (G4int ip=0; ip<n_appel; ip++) lastPAR[ip]=app_el[ip];
00425 }
00426 else
00427 {
00428 G4double a=tgZ+tgN;
00429 G4double sa=std::sqrt(a);
00430 G4double ssa=std::sqrt(sa);
00431 G4double asa=a*sa;
00432 G4double a2=a*a;
00433 G4double a3=a2*a;
00434 G4double a4=a3*a;
00435 G4double a5=a4*a;
00436 G4double a6=a4*a2;
00437 G4double a7=a6*a;
00438 G4double a8=a7*a;
00439 G4double a9=a8*a;
00440 G4double a10=a5*a5;
00441 G4double a12=a6*a6;
00442 G4double a14=a7*a7;
00443 G4double a16=a8*a8;
00444 G4double a17=a16*a;
00445
00446 G4double a32=a16*a16;
00447
00448 lastPAR[0]=.23*asa/(1.+a*.15);
00449 lastPAR[1]=2.8*asa/(1.+a*(.015+.05/ssa));
00450 lastPAR[2]=15.*a/(1.+.005*a2);
00451 lastPAR[3]=.013*a2/(1.+a3*(.006+a*.00001));
00452 lastPAR[4]=5.;
00453 lastPAR[5]=0.;
00454 lastPAR[6]=0.;
00455 lastPAR[7]=0.;
00456 lastPAR[8]=0.;
00457
00458 if(a<6.5)
00459 {
00460 G4double a28=a16*a12;
00461
00462 lastPAR[ 9]=4000*a;
00463 lastPAR[10]=1.2e7*a8+380*a17;
00464 lastPAR[11]=.7/(1.+4.e-12*a16);
00465 lastPAR[12]=2.5/a8/(a4+1.e-16*a32);
00466 lastPAR[13]=.28*a;
00467 lastPAR[14]=1.2*a2+2.3;
00468 lastPAR[15]=3.8/a;
00469
00470 lastPAR[16]=.01/(1.+.0024*a5);
00471 lastPAR[17]=.2*a;
00472 lastPAR[18]=9.e-7/(1.+.035*a5);
00473 lastPAR[19]=(42.+2.7e-11*a16)/(1.+.14*a);
00474
00475 lastPAR[20]=2.25*a3;
00476 lastPAR[21]=18.;
00477 lastPAR[22]=2.4e-3*a8/(1.+2.6e-4*a7);
00478 lastPAR[23]=3.5e-36*a32*a8/(1.+5.e-15*a32/a);
00479
00480 lastPAR[24]=1.e5/(a8+2.5e12/a16);
00481 lastPAR[25]=8.e7/(a12+1.e-27*a28*a28);
00482 lastPAR[26]=.0006*a3;
00483
00484 lastPAR[27]=10.+4.e-8*a12*a;
00485 lastPAR[28]=.114;
00486 lastPAR[29]=.003;
00487 lastPAR[30]=2.e-23;
00488
00489 lastPAR[31]=1./(1.+.0001*a8);
00490 lastPAR[32]=1.5e-4/(1.+5.e-6*a12);
00491 lastPAR[33]=.03;
00492
00493 lastPAR[34]=a/2;
00494 lastPAR[35]=2.e-7*a4;
00495 lastPAR[36]=4.;
00496 lastPAR[37]=64./a3;
00497
00498 lastPAR[38]=1.e8*std::exp(.32*asa);
00499 lastPAR[39]=20.*std::exp(.45*asa);
00500 lastPAR[40]=7.e3+2.4e6/a5;
00501 lastPAR[41]=2.5e5*std::exp(.085*a3);
00502 lastPAR[42]=2.5*a;
00503
00504 lastPAR[43]=920.+.03*a8*a3;
00505 lastPAR[44]=93.+.0023*a12;
00506 }
00507 else
00508 {
00509 G4double p1a10=2.2e-28*a10;
00510 G4double r4a16=6.e14/a16;
00511 G4double s4a16=r4a16*r4a16;
00512
00513
00514
00515 lastPAR[ 9]=4.5*std::pow(a,1.15);
00516 lastPAR[10]=.06*std::pow(a,.6);
00517 lastPAR[11]=.6*a/(1.+2.e15/a16);
00518 lastPAR[12]=.17/(a+9.e5/a3+1.5e33/a32);
00519 lastPAR[13]=(.001+7.e-11*a5)/(1.+4.4e-11*a5);
00520 lastPAR[14]=(p1a10*p1a10+2.e-29)/(1.+2.e-22*a12);
00521
00522 lastPAR[15]=400./a12+2.e-22*a9;
00523 lastPAR[16]=1.e-32*a12/(1.+5.e22/a14);
00524 lastPAR[17]=1000./a2+9.5*sa*ssa;
00525 lastPAR[18]=4.e-6*a*asa+1.e11/a16;
00526 lastPAR[19]=(120./a+.002*a2)/(1.+2.e14/a16);
00527 lastPAR[20]=9.+100./a;
00528
00529 lastPAR[21]=.002*a3+3.e7/a6;
00530 lastPAR[22]=7.e-15*a4*asa;
00531 lastPAR[23]=9000./a4;
00532
00533 lastPAR[24]=.0011*asa/(1.+3.e34/a32/a4);
00534 lastPAR[25]=1.e-5*a2+2.e14/a16;
00535 lastPAR[26]=1.2e-11*a2/(1.+1.5e19/a12);
00536 lastPAR[27]=.016*asa/(1.+5.e16/a16);
00537
00538 lastPAR[28]=.002*a4/(1.+7.e7/std::pow(a-6.83,14));
00539 lastPAR[29]=2.e6/a6+7.2/std::pow(a,.11);
00540 lastPAR[30]=11.*a3/(1.+7.e23/a16/a8);
00541 lastPAR[31]=100./asa;
00542
00543 lastPAR[32]=(.1+4.4e-5*a2)/(1.+5.e5/a4);
00544 lastPAR[33]=3.5e-4*a2/(1.+1.e8/a8);
00545 lastPAR[34]=1.3+3.e5/a4;
00546 lastPAR[35]=500./(a2+50.)+3;
00547 lastPAR[36]=1.e-9/a+s4a16*s4a16;
00548
00549 lastPAR[37]=.4*asa+3.e-9*a6;
00550 lastPAR[38]=.0005*a5;
00551 lastPAR[39]=.002*a5;
00552 lastPAR[40]=10.;
00553
00554 lastPAR[41]=.05+.005*a;
00555 lastPAR[42]=7.e-8/sa;
00556 lastPAR[43]=.8*sa;
00557 lastPAR[44]=.02*sa;
00558 lastPAR[45]=1.e8/a3;
00559 lastPAR[46]=3.e32/(a32+1.e32);
00560
00561 lastPAR[47]=24.;
00562 lastPAR[48]=20./sa;
00563 lastPAR[49]=7.e3*a/(sa+1.);
00564 lastPAR[50]=900.*sa/(1.+500./a3);
00565 }
00566
00567 lastPAR[51]=1.e15+2.e27/a4/(1.+2.e-18*a16);
00568 }
00569 lastPAR[nLast]=pwd;
00570
00571 G4double lp=lPMin;
00572 G4bool memCS=onlyCS;
00573 onlyCS=false;
00574 lastCST[0]=GetTabValues(lp, PDG, tgZ, tgN);
00575 onlyCS=memCS;
00576 lastSST[0]=theSS;
00577 lastS1T[0]=theS1;
00578 lastB1T[0]=theB1;
00579 lastS2T[0]=theS2;
00580 lastB2T[0]=theB2;
00581 lastS3T[0]=theS3;
00582 lastB3T[0]=theB3;
00583 lastS4T[0]=theS4;
00584 lastB4T[0]=theB4;
00585 }
00586 if(LP>ILP)
00587 {
00588 G4int ini = static_cast<int>((ILP-lPMin+.000001)/dlnP)+1;
00589 if(ini<0) ini=0;
00590 if(ini<nPoints)
00591 {
00592 G4int fin = static_cast<int>((LP-lPMin)/dlnP)+1;
00593 if(fin>=nPoints) fin=nLast;
00594 if(fin>=ini)
00595 {
00596 G4double lp=0.;
00597 for(G4int ip=ini; ip<=fin; ip++)
00598 {
00599 lp=lPMin+ip*dlnP;
00600 G4bool memCS=onlyCS;
00601 onlyCS=false;
00602 lastCST[ip]=GetTabValues(lp, PDG, tgZ, tgN);
00603 onlyCS=memCS;
00604 lastSST[ip]=theSS;
00605 lastS1T[ip]=theS1;
00606 lastB1T[ip]=theB1;
00607 lastS2T[ip]=theS2;
00608 lastB2T[ip]=theB2;
00609 lastS3T[ip]=theS3;
00610 lastB3T[ip]=theB3;
00611 lastS4T[ip]=theS4;
00612 lastB4T[ip]=theB4;
00613 }
00614 return lp;
00615 }
00616 else G4cout<<"*Warning*G4ChipsAntiBaryonElasticXS::GetPTables: PDG="<<PDG
00617 <<", Z="<<tgZ<<", N="<<tgN<<", i="<<ini<<" > fin="<<fin<<", LP="<<LP
00618 <<" > ILP="<<ILP<<" nothing is done!"<<G4endl;
00619 }
00620 else G4cout<<"*Warning*G4ChipsAntiBaryonElasticXS::GetPTables: PDG="<<PDG
00621 <<", Z="<<tgZ<<", N="<<tgN<<", i="<<ini<<">= max="<<nPoints<<", LP="<<LP
00622 <<" > ILP="<<ILP<<", lPMax="<<lPMax<<" nothing is done!"<<G4endl;
00623 }
00624 }
00625 else
00626 {
00627
00628
00629
00630 G4ExceptionDescription ed;
00631 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
00632 << ", while it is defined only for Anti Baryons" << G4endl;
00633 G4Exception("G4ChipsAntiBaryonElasticXS::GetPTables()", "HAD_CHPS_0000",
00634 FatalException, ed);
00635 }
00636 return ILP;
00637 }
00638
00639
00640 G4double G4ChipsAntiBaryonElasticXS::GetExchangeT(G4int tgZ, G4int tgN, G4int PDG)
00641 {
00642 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
00643 static const G4double third=1./3.;
00644 static const G4double fifth=1./5.;
00645 static const G4double sevth=1./7.;
00646
00647 if(PDG<-3334 || PDG>-1111)G4cout<<"*Warning*G4QAntiBaryonElCS::GetExT:PDG="<<PDG<<G4endl;
00648 if(onlyCS)G4cout<<"WarningG4ChipsAntiBaryonElasticXS::GetExchanT:onlyCS=1"<<G4endl;
00649 if(lastLP<-4.3) return lastTM*GeVSQ*G4UniformRand();
00650 G4double q2=0.;
00651 if(tgZ==1 && tgN==0)
00652 {
00653 G4double E1=lastTM*theB1;
00654 G4double R1=(1.-std::exp(-E1));
00655 G4double E2=lastTM*theB2;
00656 G4double R2=(1.-std::exp(-E2*E2*E2));
00657 G4double E3=lastTM*theB3;
00658 G4double R3=(1.-std::exp(-E3));
00659 G4double I1=R1*theS1/theB1;
00660 G4double I2=R2*theS2;
00661 G4double I3=R3*theS3;
00662 G4double I12=I1+I2;
00663 G4double rand=(I12+I3)*G4UniformRand();
00664 if (rand<I1 )
00665 {
00666 G4double ran=R1*G4UniformRand();
00667 if(ran>1.) ran=1.;
00668 q2=-std::log(1.-ran)/theB1;
00669 }
00670 else if(rand<I12)
00671 {
00672 G4double ran=R2*G4UniformRand();
00673 if(ran>1.) ran=1.;
00674 q2=-std::log(1.-ran);
00675 if(q2<0.) q2=0.;
00676 q2=std::pow(q2,third)/theB2;
00677 }
00678 else
00679 {
00680 G4double ran=R3*G4UniformRand();
00681 if(ran>1.) ran=1.;
00682 q2=-std::log(1.-ran)/theB3;
00683 }
00684 }
00685 else
00686 {
00687 G4double a=tgZ+tgN;
00688 G4double E1=lastTM*(theB1+lastTM*theSS);
00689 G4double R1=(1.-std::exp(-E1));
00690 G4double tss=theSS+theSS;
00691 G4double tm2=lastTM*lastTM;
00692 G4double E2=lastTM*tm2*theB2;
00693 if(a>6.5)E2*=tm2;
00694 G4double R2=(1.-std::exp(-E2));
00695 G4double E3=lastTM*theB3;
00696 if(a>6.5)E3*=tm2*tm2*tm2;
00697 G4double R3=(1.-std::exp(-E3));
00698 G4double E4=lastTM*theB4;
00699 G4double R4=(1.-std::exp(-E4));
00700 G4double I1=R1*theS1;
00701 G4double I2=R2*theS2;
00702 G4double I3=R3*theS3;
00703 G4double I4=R4*theS4;
00704 G4double I12=I1+I2;
00705 G4double I13=I12+I3;
00706 G4double rand=(I13+I4)*G4UniformRand();
00707 if(rand<I1)
00708 {
00709 G4double ran=R1*G4UniformRand();
00710 if(ran>1.) ran=1.;
00711 q2=-std::log(1.-ran)/theB1;
00712 if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
00713 }
00714 else if(rand<I12)
00715 {
00716 G4double ran=R2*G4UniformRand();
00717 if(ran>1.) ran=1.;
00718 q2=-std::log(1.-ran)/theB2;
00719 if(q2<0.) q2=0.;
00720 if(a<6.5) q2=std::pow(q2,third);
00721 else q2=std::pow(q2,fifth);
00722 }
00723 else if(rand<I13)
00724 {
00725 G4double ran=R3*G4UniformRand();
00726 if(ran>1.) ran=1.;
00727 q2=-std::log(1.-ran)/theB3;
00728 if(q2<0.) q2=0.;
00729 if(a>6.5) q2=std::pow(q2,sevth);
00730 }
00731 else
00732 {
00733 G4double ran=R4*G4UniformRand();
00734 if(ran>1.) ran=1.;
00735 q2=-std::log(1.-ran)/theB4;
00736 if(a<6.5) q2=lastTM-q2;
00737 }
00738 }
00739 if(q2<0.) q2=0.;
00740 if(!(q2>=-1.||q2<=1.))G4cout<<"*NAN*G4QaBElasticCrossSect::GetExchangeT:-t="<<q2<<G4endl;
00741 if(q2>lastTM)
00742 {
00743 q2=lastTM;
00744 }
00745 return q2*GeVSQ;
00746 }
00747
00748
00749 G4double G4ChipsAntiBaryonElasticXS::GetSlope(G4int tgZ, G4int tgN, G4int PDG)
00750 {
00751 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
00752 if(onlyCS)G4cout<<"WarningG4ChipsAntiBaryonElasticXS::GetSlope:onlCS=true"<<G4endl;
00753 if(lastLP<-4.3) return 0.;
00754 if(PDG<-3334 || PDG>-1111)
00755 {
00756
00757
00758
00759 G4ExceptionDescription ed;
00760 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
00761 << ", while it is defined only for Anti Baryons" << G4endl;
00762 G4Exception("G4ChipsAntiBaryonElasticXS::GetSlope()", "HAD_CHPS_0000",
00763 FatalException, ed);
00764 }
00765 if(theB1<0.) theB1=0.;
00766 if(!(theB1>=-1.||theB1<=1.))G4cout<<"*NAN*G4QaBaElasticCrossS::Getslope:"<<theB1<<G4endl;
00767 return theB1/GeVSQ;
00768 }
00769
00770
00771 G4double G4ChipsAntiBaryonElasticXS::GetHMaxT()
00772 {
00773 static const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
00774 return lastTM*HGeVSQ;
00775 }
00776
00777
00778 G4double G4ChipsAntiBaryonElasticXS::GetTabValues(G4double lp, G4int PDG, G4int tgZ,
00779 G4int tgN)
00780 {
00781 if(PDG<-3334 || PDG>-1111) G4cout<<"*Warning*G4QAntiBaryElCS::GetTabV:PDG="<<PDG<<G4endl;
00782 if(tgZ<0 || tgZ>92)
00783 {
00784 G4cout<<"*Warning*G4QAntiBaryonElCS::GetTabValue:(1-92) NoIsotopesFor Z="<<tgZ<<G4endl;
00785 return 0.;
00786 }
00787 G4int iZ=tgZ-1;
00788 if(iZ<0)
00789 {
00790 iZ=0;
00791 tgZ=1;
00792 tgN=0;
00793 }
00794 G4double p=std::exp(lp);
00795 G4double sp=std::sqrt(p);
00796 G4double p2=p*p;
00797 G4double p3=p2*p;
00798 G4double p4=p3*p;
00799 if ( tgZ == 1 && tgN == 0 )
00800 {
00801 G4double dl2=lp-lastPAR[6];
00802 theSS=lastPAR[29];
00803 theS1=(lastPAR[7]+lastPAR[8]*dl2*dl2)/(1.+lastPAR[9]/p4/p)+
00804 (lastPAR[10]/p2+lastPAR[11]*p)/(p4+lastPAR[12]*sp);
00805 theB1=lastPAR[13]*std::pow(p,lastPAR[14])/(1.+lastPAR[15]/p3);
00806 theS2=lastPAR[16]+lastPAR[17]/(p4+lastPAR[18]*p);
00807 theB2=lastPAR[19]+lastPAR[20]/(p4+lastPAR[21]/sp);
00808 theS3=lastPAR[22]+lastPAR[23]/(p4*p4+lastPAR[24]*p2+lastPAR[25]);
00809 theB3=lastPAR[26]+lastPAR[27]/(p4+lastPAR[28]);
00810 theS4=0.;
00811 theB4=0.;
00812
00813 G4double ye=std::exp(lp*lastPAR[0]);
00814 G4double dp=lp-lastPAR[1];
00815 return lastPAR[2]/(ye+lastPAR[3])+lastPAR[4]*dp*dp+lastPAR[5];
00816 }
00817 else
00818 {
00819 G4double p5=p4*p;
00820 G4double p6=p5*p;
00821 G4double p8=p6*p2;
00822 G4double p10=p8*p2;
00823 G4double p12=p10*p2;
00824 G4double p16=p8*p8;
00825
00826 G4double dl=lp-5.;
00827 G4double a=tgZ+tgN;
00828 G4double pah=std::pow(p,a/2);
00829 G4double pa=pah*pah;
00830 G4double pa2=pa*pa;
00831 if(a<6.5)
00832 {
00833 theS1=lastPAR[9]/(1.+lastPAR[10]*p4*pa)+lastPAR[11]/(p4+lastPAR[12]*p4/pa2)+
00834 (lastPAR[13]*dl*dl+lastPAR[14])/(1.+lastPAR[15]/p2);
00835 theB1=(lastPAR[16]+lastPAR[17]*p2)/(p4+lastPAR[18]/pah)+lastPAR[19];
00836 theSS=lastPAR[20]/(1.+lastPAR[21]/p2)+lastPAR[22]/(p6/pa+lastPAR[23]/p16);
00837 theS2=lastPAR[24]/(pa/p2+lastPAR[25]/p4)+lastPAR[26];
00838 theB2=lastPAR[27]*std::pow(p,lastPAR[28])+lastPAR[29]/(p8+lastPAR[30]/p16);
00839 theS3=lastPAR[31]/(pa*p+lastPAR[32]/pa)+lastPAR[33];
00840 theB3=lastPAR[34]/(p3+lastPAR[35]/p6)+lastPAR[36]/(1.+lastPAR[37]/p2);
00841 theS4=p2*(pah*lastPAR[38]*std::exp(-pah*lastPAR[39])+
00842 lastPAR[40]/(1.+lastPAR[41]*std::pow(p,lastPAR[42])));
00843 theB4=lastPAR[43]*pa/p2/(1.+pa*lastPAR[44]);
00844 }
00845 else
00846 {
00847 theS1=lastPAR[9]/(1.+lastPAR[10]/p4)+lastPAR[11]/(p4+lastPAR[12]/p2)+
00848 lastPAR[13]/(p5+lastPAR[14]/p16);
00849 theB1=(lastPAR[15]/p8+lastPAR[19])/(p+lastPAR[16]/std::pow(p,lastPAR[20]))+
00850 lastPAR[17]/(1.+lastPAR[18]/p4);
00851 theSS=lastPAR[21]/(p4/std::pow(p,lastPAR[23])+lastPAR[22]/p4);
00852 theS2=lastPAR[24]/p4/(std::pow(p,lastPAR[25])+lastPAR[26]/p12)+lastPAR[27];
00853 theB2=lastPAR[28]/std::pow(p,lastPAR[29])+lastPAR[30]/std::pow(p,lastPAR[31]);
00854 theS3=lastPAR[32]/std::pow(p,lastPAR[35])/(1.+lastPAR[36]/p12)+
00855 lastPAR[33]/(1.+lastPAR[34]/p6);
00856 theB3=lastPAR[37]/p8+lastPAR[38]/p2+lastPAR[39]/(1.+lastPAR[40]/p8);
00857 theS4=(lastPAR[41]/p4+lastPAR[46]/p)/(1.+lastPAR[42]/p10)+
00858 (lastPAR[43]+lastPAR[44]*dl*dl)/(1.+lastPAR[45]/p12);
00859 theB4=lastPAR[47]/(1.+lastPAR[48]/p)+lastPAR[49]*p4/(1.+lastPAR[50]*p5);
00860 }
00861
00862 G4double dlp=lp-lastPAR[4];
00863
00864 return (lastPAR[0]*dlp*dlp+lastPAR[1]+lastPAR[2]/p)/(1.+lastPAR[3]/p);
00865 }
00866 return 0.;
00867 }
00868
00869
00870 G4double G4ChipsAntiBaryonElasticXS::GetQ2max(G4int PDG, G4int tgZ, G4int tgN,
00871 G4double pP)
00872 {
00873 static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass()*.001;
00874 static const G4double mProt= G4Proton::Proton()->GetPDGMass()*.001;
00875 static const G4double mNuc2= sqr((mProt+mNeut)/2);
00876 G4double pP2=pP*pP;
00877 if(tgZ || tgN>-1)
00878 {
00879 G4double mt=G4ParticleTable::GetParticleTable()->FindIon(tgZ,tgZ+tgN,0,tgZ)->GetPDGMass()*.001;
00880 G4double dmt=mt+mt;
00881 G4double mds=dmt*std::sqrt(pP2+mNuc2)+mNuc2+mt*mt;
00882 return dmt*dmt*pP2/mds;
00883 }
00884 else
00885 {
00886
00887
00888
00889 G4ExceptionDescription ed;
00890 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
00891 << ", while it is defined only for p projectiles & Z_target>0" << G4endl;
00892 G4Exception("G4ChipsAntiBaryonElasticXS::GetQ2max()", "HAD_CHPS_0000",
00893 FatalException, ed);
00894 return 0;
00895 }
00896 }