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 "G4ChipsKaonPlusElasticXS.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4DynamicParticle.hh"
00043 #include "G4ParticleDefinition.hh"
00044 #include "G4KaonPlus.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(G4ChipsKaonPlusElasticXS);
00053
00054 G4ChipsKaonPlusElasticXS::G4ChipsKaonPlusElasticXS():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 G4ChipsKaonPlusElasticXS::~G4ChipsKaonPlusElasticXS()
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 G4ChipsKaonPlusElasticXS::IsIsoApplicable(const G4DynamicParticle* Pt, G4int, G4int,
00133 const G4Element*,
00134 const G4Material*)
00135 {
00136 G4ParticleDefinition* particle = Pt->GetDefinition();
00137 if (particle == G4KaonPlus::KaonPlus() ) return true;
00138 return false;
00139 }
00140
00141
00142
00143 G4double G4ChipsKaonPlusElasticXS::GetIsoCrossSection(const G4DynamicParticle* Pt, G4int tgZ, G4int A,
00144 const G4Isotope*,
00145 const G4Element*,
00146 const G4Material*)
00147 {
00148 G4double pMom=Pt->GetTotalMomentum();
00149 G4int tgN = A - tgZ;
00150
00151 return GetChipsCrossSection(pMom, tgZ, tgN, 321);
00152 }
00153
00154 G4double G4ChipsKaonPlusElasticXS::GetChipsCrossSection(G4double pMom, G4int tgZ, G4int tgN, G4int)
00155 {
00156 static std::vector <G4int> colN;
00157 static std::vector <G4int> colZ;
00158 static std::vector <G4double> colP;
00159 static std::vector <G4double> colTH;
00160 static std::vector <G4double> colCS;
00161
00162
00163 G4bool fCS = false;
00164 G4double pEn=pMom;
00165 onlyCS=fCS;
00166
00167 G4bool in=false;
00168 lastP = 0.;
00169 lastN = tgN;
00170 lastZ = tgZ;
00171 lastI = colN.size();
00172 if(lastI) for(G4int i=0; i<lastI; i++)
00173 {
00174 if(colN[i]==tgN && colZ[i]==tgZ)
00175 {
00176 lastI=i;
00177 lastTH =colTH[i];
00178 if(pEn<=lastTH)
00179 {
00180 return 0.;
00181 }
00182 lastP =colP [i];
00183 lastCS =colCS[i];
00184
00185 if(lastP == pMom)
00186 {
00187 CalculateCrossSection(fCS,-1,i,321,lastZ,lastN,pMom);
00188 return lastCS*millibarn;
00189 }
00190 in = true;
00191
00192 lastCS=CalculateCrossSection(fCS,-1,i,321,lastZ,lastN,pMom);
00193 if(lastCS<=0. && pEn>lastTH)
00194 {
00195 lastTH=pEn;
00196 }
00197 break;
00198 }
00199 }
00200 if(!in)
00201 {
00203 lastCS=CalculateCrossSection(fCS,0,lastI,321,lastZ,lastN,pMom);
00204 if(lastCS<=0.)
00205 {
00206 lastTH = 0;
00207 if(pEn>lastTH)
00208 {
00209 lastTH=pEn;
00210 }
00211 }
00212 colN.push_back(tgN);
00213 colZ.push_back(tgZ);
00214 colP.push_back(pMom);
00215 colTH.push_back(lastTH);
00216 colCS.push_back(lastCS);
00217 return lastCS*millibarn;
00218 }
00219 else
00220 {
00221 colP[lastI]=pMom;
00222 colCS[lastI]=lastCS;
00223 }
00224 return lastCS*millibarn;
00225 }
00226
00227
00228
00229 G4double G4ChipsKaonPlusElasticXS::CalculateCrossSection(G4bool CS, G4int F,
00230 G4int I, G4int PDG, G4int tgZ, G4int tgN, G4double pIU)
00231 {
00232
00233 static std::vector <G4double> PIN;
00234
00235 G4double pMom=pIU/GeV;
00236 onlyCS=CS;
00237 lastLP=std::log(pMom);
00238 if(F)
00239 {
00240 if(F<0)
00241 {
00242 lastPIN = PIN[I];
00243 lastPAR = PAR[I];
00244 lastCST = CST[I];
00245 lastSST = SST[I];
00246 lastS1T = S1T[I];
00247 lastB1T = B1T[I];
00248 lastS2T = S2T[I];
00249 lastB2T = B2T[I];
00250 lastS3T = S3T[I];
00251 lastB3T = B3T[I];
00252 lastS4T = S4T[I];
00253 lastB4T = B4T[I];
00254 }
00255 if(lastLP>lastPIN && lastLP<lPMax)
00256 {
00257 lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
00258 PIN[I]=lastPIN;
00259 }
00260 }
00261 else
00262 {
00263 lastPAR = new G4double[nPoints];
00264 lastPAR[nLast]=0;
00265 lastCST = new G4double[nPoints];
00266 lastSST = new G4double[nPoints];
00267 lastS1T = new G4double[nPoints];
00268 lastB1T = new G4double[nPoints];
00269 lastS2T = new G4double[nPoints];
00270 lastB2T = new G4double[nPoints];
00271 lastS3T = new G4double[nPoints];
00272 lastB3T = new G4double[nPoints];
00273 lastS4T = new G4double[nPoints];
00274 lastB4T = new G4double[nPoints];
00275 lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN);
00276 PIN.push_back(lastPIN);
00277 PAR.push_back(lastPAR);
00278 CST.push_back(lastCST);
00279 SST.push_back(lastSST);
00280 S1T.push_back(lastS1T);
00281 B1T.push_back(lastB1T);
00282 S2T.push_back(lastS2T);
00283 B2T.push_back(lastB2T);
00284 S3T.push_back(lastS3T);
00285 B3T.push_back(lastB3T);
00286 S4T.push_back(lastS4T);
00287 B4T.push_back(lastB4T);
00288 }
00289
00290 if(lastLP>lastPIN && lastLP<lPMax)
00291 {
00292 lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
00293 }
00294 if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom);
00295 if(lastLP>lPMin && lastLP<=lastPIN)
00296 {
00297 if(lastLP==lastPIN)
00298 {
00299 G4double shift=(lastLP-lPMin)/dlnP+.000001;
00300 G4int blast=static_cast<int>(shift);
00301 if(blast<0 || blast>=nLast) G4cout<<"G4QKPElCS::CCS:b="<<blast<<",n="<<nLast<<G4endl;
00302 lastSIG = lastCST[blast];
00303 if(!onlyCS)
00304 {
00305 theSS = lastSST[blast];
00306 theS1 = lastS1T[blast];
00307 theB1 = lastB1T[blast];
00308 theS2 = lastS2T[blast];
00309 theB2 = lastB2T[blast];
00310 theS3 = lastS3T[blast];
00311 theB3 = lastB3T[blast];
00312 theS4 = lastS4T[blast];
00313 theB4 = lastB4T[blast];
00314 }
00315 }
00316 else
00317 {
00318 G4double shift=(lastLP-lPMin)/dlnP;
00319 G4int blast=static_cast<int>(shift);
00320 if(blast<0) blast=0;
00321 if(blast>=nLast) blast=nLast-1;
00322 shift-=blast;
00323 G4int lastL=blast+1;
00324 G4double SIGL=lastCST[blast];
00325 lastSIG= SIGL+shift*(lastCST[lastL]-SIGL);
00326 if(!onlyCS)
00327 {
00328 G4double SSTL=lastSST[blast];
00329 theSS=SSTL+shift*(lastSST[lastL]-SSTL);
00330 G4double S1TL=lastS1T[blast];
00331 theS1=S1TL+shift*(lastS1T[lastL]-S1TL);
00332 G4double B1TL=lastB1T[blast];
00333 theB1=B1TL+shift*(lastB1T[lastL]-B1TL);
00334 G4double S2TL=lastS2T[blast];
00335 theS2=S2TL+shift*(lastS2T[lastL]-S2TL);
00336 G4double B2TL=lastB2T[blast];
00337 theB2=B2TL+shift*(lastB2T[lastL]-B2TL);
00338 G4double S3TL=lastS3T[blast];
00339 theS3=S3TL+shift*(lastS3T[lastL]-S3TL);
00340 G4double B3TL=lastB3T[blast];
00341 theB3=B3TL+shift*(lastB3T[lastL]-B3TL);
00342 G4double S4TL=lastS4T[blast];
00343 theS4=S4TL+shift*(lastS4T[lastL]-S4TL);
00344 G4double B4TL=lastB4T[blast];
00345 theB4=B4TL+shift*(lastB4T[lastL]-B4TL);
00346 }
00347 }
00348 }
00349 else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN);
00350 if(lastSIG<0.) lastSIG = 0.;
00351 return lastSIG;
00352 }
00353
00354
00355 G4double G4ChipsKaonPlusElasticXS::GetPTables(G4double LP,G4double ILP, G4int PDG,
00356 G4int tgZ, G4int tgN)
00357 {
00358
00359 static const G4double pwd=2727;
00360 const G4int n_kppel=35;
00361
00362 G4double kpp_el[n_kppel]={.7,.38,.0676,.0557,3.5,2.23,.7,.1,2.,1.,.372,5.,74.,3.,3.4,
00363 .2,.17,.001,8.,.055,3.64,5.e-5,4000.,1500.,.46,1.2e6,3.5e6,
00364 5.e-5,1.e10,8.5e8,1.e10,1.1,3.4e6,6.8e6,0.};
00365
00366
00367
00368 if(PDG == 321)
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_kppel; ip++) lastPAR[ip]=kpp_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]=.06*asa/(1.+a*(.01+.1/ssa));
00410 lastPAR[1]=.75*asa/(1.+.009*a);
00411 lastPAR[2]=.9*asa*ssa/(1.+.03*a);
00412 lastPAR[3]=3.;
00413 lastPAR[4]=4.2;
00414 lastPAR[5]=0.;
00415 lastPAR[6]=0.;
00416 lastPAR[7]=0.;
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*G4ChipsKaonPlusElasticXS::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*G4ChipsKaonPlusElasticXS::GetPTables: PDG="<<PDG
00582 <<", Z="<<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=321 (K+) " << G4endl;
00594 G4Exception("G4ChipsKaonPlusElasticXS::GetPTables()", "HAD_CHPS_0000",
00595 FatalException, ed);
00596 }
00597 return ILP;
00598 }
00599
00600
00601 G4double G4ChipsKaonPlusElasticXS::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!=321) G4cout<<"*Warning*G4ChipsKaonPlusElasticXS::GetExT:PDG="<<PDG<<G4endl;
00608 if(onlyCS) G4cout<<"*Warning*G4ChipsKaonPlusElasticXS::GetExT: 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*G4QKaonPlusElasticCS::GetExchT: -t="<<q2<<G4endl;
00701 if(q2>lastTM)
00702 {
00703 q2=lastTM;
00704 }
00705 return q2*GeVSQ;
00706 }
00707
00708
00709 G4double G4ChipsKaonPlusElasticXS::GetSlope(G4int tgZ, G4int tgN, G4int PDG)
00710 {
00711 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
00712 if(onlyCS)G4cout<<"*Warning*G4ChipsKaonPlusElasticXS::GetSl:onlCS=true"<<G4endl;
00713 if(lastLP<-4.3) return 0.;
00714 if(PDG != 321)
00715 {
00716 G4ExceptionDescription ed;
00717 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
00718 << ", while it is defined only for PDG=321 (K+)" << G4endl;
00719 G4Exception("G4ChipsKaonPlusElasticXS::GetSlope()", "HAD_CHPS_0000",
00720 FatalException, ed);
00721 }
00722 if(theB1<0.) theB1=0.;
00723 if(!(theB1>=-1.||theB1<=1.))G4cout<<"*NAN*G4QKaonPlusElCS::GetSlope:B1="<<theB1<<G4endl;
00724 return theB1/GeVSQ;
00725 }
00726
00727
00728 G4double G4ChipsKaonPlusElasticXS::GetHMaxT()
00729 {
00730 static const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
00731 return lastTM*HGeVSQ;
00732 }
00733
00734
00735 G4double G4ChipsKaonPlusElasticXS::GetTabValues(G4double lp, G4int PDG, G4int tgZ,
00736 G4int tgN)
00737 {
00738 if(PDG!=321)G4cout<<"*Warning*G4ChipsKaonPlusElasticXS::GetTaV:PDG="<<PDG<<G4endl;
00739 if(tgZ<0 || tgZ>92)
00740 {
00741 G4cout<<"*Warning*G4QKaonPlusElasticCS::GetTabV:(1-92)NoIsotopes for Z="<<tgZ<<G4endl;
00742 return 0.;
00743 }
00744 G4int iZ=tgZ-1;
00745 if(iZ<0)
00746 {
00747 iZ=0;
00748 tgZ=1;
00749 tgN=0;
00750 }
00751 G4double p=std::exp(lp);
00752 G4double sp=std::sqrt(p);
00753 G4double p2=p*p;
00754 G4double p3=p2*p;
00755 G4double p4=p3*p;
00756 if ( tgZ == 1 && tgN == 0 )
00757 {
00758 G4double dl2=lp-lastPAR[11];
00759 theSS=lastPAR[34];
00760 theS1=(lastPAR[12]+lastPAR[13]*dl2*dl2)/(1.+lastPAR[14]/p4/p)+
00761 (lastPAR[15]/p2+lastPAR[16]*p)/(p4+lastPAR[17]*sp);
00762 theB1=lastPAR[18]*std::pow(p,lastPAR[19])/(1.+lastPAR[20]/p3);
00763 theS2=lastPAR[21]+lastPAR[22]/(p4+lastPAR[23]*p);
00764 theB2=lastPAR[24]+lastPAR[25]/(p4+lastPAR[26]/sp);
00765 theS3=lastPAR[27]+lastPAR[28]/(p4*p4+lastPAR[29]*p2+lastPAR[30]);
00766 theB3=lastPAR[31]+lastPAR[32]/(p4+lastPAR[33]);
00767 theS4=0.;
00768 theB4=0.;
00769
00770 G4double dp=lp-lastPAR[4];
00771
00772 return lastPAR[0]/(lastPAR[2]+sqr(p-lastPAR[1]))+(lastPAR[3]*dp*dp+lastPAR[5])/
00773 (1.-lastPAR[6]/sp+lastPAR[7]/p4)
00774 +lastPAR[8]/(sqr(p-lastPAR[9])+lastPAR[10]);
00775
00776 }
00777 else
00778 {
00779 G4double p5=p4*p;
00780 G4double p6=p5*p;
00781 G4double p8=p6*p2;
00782 G4double p10=p8*p2;
00783 G4double p12=p10*p2;
00784 G4double p16=p8*p8;
00785
00786 G4double dl=lp-5.;
00787 G4double a=tgZ+tgN;
00788 G4double pah=std::pow(p,a/2);
00789 G4double pa=pah*pah;
00790 G4double pa2=pa*pa;
00791 if(a<6.5)
00792 {
00793 theS1=lastPAR[9]/(1.+lastPAR[10]*p4*pa)+lastPAR[11]/(p4+lastPAR[12]*p4/pa2)+
00794 (lastPAR[13]*dl*dl+lastPAR[14])/(1.+lastPAR[15]/p2);
00795 theB1=(lastPAR[16]+lastPAR[17]*p2)/(p4+lastPAR[18]/pah)+lastPAR[19];
00796 theSS=lastPAR[20]/(1.+lastPAR[21]/p2)+lastPAR[22]/(p6/pa+lastPAR[23]/p16);
00797 theS2=lastPAR[24]/(pa/p2+lastPAR[25]/p4)+lastPAR[26];
00798 theB2=lastPAR[27]*std::pow(p,lastPAR[28])+lastPAR[29]/(p8+lastPAR[30]/p16);
00799 theS3=lastPAR[31]/(pa*p+lastPAR[32]/pa)+lastPAR[33];
00800 theB3=lastPAR[34]/(p3+lastPAR[35]/p6)+lastPAR[36]/(1.+lastPAR[37]/p2);
00801 theS4=p2*(pah*lastPAR[38]*std::exp(-pah*lastPAR[39])+
00802 lastPAR[40]/(1.+lastPAR[41]*std::pow(p,lastPAR[42])));
00803 theB4=lastPAR[43]*pa/p2/(1.+pa*lastPAR[44]);
00804 }
00805 else
00806 {
00807 theS1=lastPAR[9]/(1.+lastPAR[10]/p4)+lastPAR[11]/(p4+lastPAR[12]/p2)+
00808 lastPAR[13]/(p5+lastPAR[14]/p16);
00809 theB1=(lastPAR[15]/p8+lastPAR[19])/(p+lastPAR[16]/std::pow(p,lastPAR[20]))+
00810 lastPAR[17]/(1.+lastPAR[18]/p4);
00811 theSS=lastPAR[21]/(p4/std::pow(p,lastPAR[23])+lastPAR[22]/p4);
00812 theS2=lastPAR[24]/p4/(std::pow(p,lastPAR[25])+lastPAR[26]/p12)+lastPAR[27];
00813 theB2=lastPAR[28]/std::pow(p,lastPAR[29])+lastPAR[30]/std::pow(p,lastPAR[31]);
00814 theS3=lastPAR[32]/std::pow(p,lastPAR[35])/(1.+lastPAR[36]/p12)+
00815 lastPAR[33]/(1.+lastPAR[34]/p6);
00816 theB3=lastPAR[37]/p8+lastPAR[38]/p2+lastPAR[39]/(1.+lastPAR[40]/p8);
00817 theS4=(lastPAR[41]/p4+lastPAR[46]/p)/(1.+lastPAR[42]/p10)+
00818 (lastPAR[43]+lastPAR[44]*dl*dl)/(1.+lastPAR[45]/p12);
00819 theB4=lastPAR[47]/(1.+lastPAR[48]/p)+lastPAR[49]*p4/(1.+lastPAR[50]*p5);
00820 }
00821
00822 G4double dlp=lp-lastPAR[4];
00823
00824 return (lastPAR[0]*dlp*dlp+lastPAR[1]+lastPAR[2]/p2)/(1.+lastPAR[3]/p2/sp);
00825 }
00826 return 0.;
00827 }
00828
00829
00830 G4double G4ChipsKaonPlusElasticXS::GetQ2max(G4int PDG, G4int tgZ, G4int tgN,
00831 G4double pP)
00832 {
00833 static const G4double mK= G4KaonPlus::KaonPlus()->GetPDGMass()*.001;
00834 static const G4double mK2= mK*mK;
00835
00836 G4double pP2=pP*pP;
00837 if(tgZ || tgN>-1)
00838 {
00839 G4double mt=G4ParticleTable::GetParticleTable()->FindIon(tgZ,tgZ+tgN,0,tgZ)->GetPDGMass()*.001;
00840
00841 G4double dmt=mt+mt;
00842 G4double mds=dmt*std::sqrt(pP2+mK2)+mK2+mt*mt;
00843 return dmt*dmt*pP2/mds;
00844 }
00845 else
00846 {
00847 G4ExceptionDescription ed;
00848 ed << "PDG = " << PDG << ",Z = " << tgZ << ", N = " << tgN
00849 << ", while it is defined only for p projectiles & Z_target>0" << G4endl;
00850 G4Exception("G4ChipsKaonPlusElasticXS::GetQ2max()", "HAD_CHPS_0000",
00851 FatalException, ed);
00852 return 0;
00853 }
00854 }