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