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