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 "G4ChipsKaonMinusElasticXS.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4DynamicParticle.hh"
00043 #include "G4ParticleDefinition.hh"
00044 #include "G4KaonMinus.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(G4ChipsKaonMinusElasticXS);
00053
00054 G4ChipsKaonMinusElasticXS::G4ChipsKaonMinusElasticXS():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 G4ChipsKaonMinusElasticXS::~G4ChipsKaonMinusElasticXS()
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 G4ChipsKaonMinusElasticXS::IsIsoApplicable(const G4DynamicParticle* Pt, G4int, G4int,
00133 const G4Element*,
00134 const G4Material*)
00135 {
00136 G4ParticleDefinition* particle = Pt->GetDefinition();
00137 if (particle == G4KaonMinus::KaonMinus() ) return true;
00138 return false;
00139 }
00140
00141
00142
00143 G4double G4ChipsKaonMinusElasticXS::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 G4ChipsKaonMinusElasticXS::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
00165 G4double pEn=pMom;
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,-321,lastZ,lastN,pMom);
00189 return lastCS*millibarn;
00190 }
00191 in = true;
00192
00193 lastCS=CalculateCrossSection(fCS,-1,i,-321,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,-321,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 G4ChipsKaonMinusElasticXS::CalculateCrossSection(G4bool CS, G4int F,
00231 G4int I, 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<<"G4QKMElCS::CCS:b="<<blast<<",n="<<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 G4ChipsKaonMinusElasticXS::GetPTables(G4double LP,G4double ILP, G4int PDG,
00357 G4int tgZ, G4int tgN)
00358 {
00359
00360 static const G4double pwd=2727;
00361 const G4int n_kmpel=36;
00362
00363 G4double kmp_el[n_kmpel]={5.2,.0557,3.5,2.23,.7,.075,.004,.39,.000156,.15,1.,.0156,5.,
00364 74.,3.,3.4,.2,.17,.001,8.,.055,3.64,5.e-5,4000.,1500.,.46,
00365 1.2e6,3.5e6,5.e-5,1.e10,8.5e8,1.e10,1.1,3.4e6,6.8e6,0.};
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_kmpel; ip++) lastPAR[ip]=kmp_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]=.1*a2*ssa/(1.+.0015*a2/ssa);
00412 lastPAR[3]=1./(1.+500./a2);
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*G4ChipsKaonMinusElasticXS::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*G4ChipsKaonMinusElasticXS::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("G4ChipsKaonMinusElasticXS::GetPTables()", "HAD_CHPS_0000",
00595 FatalException, ed);
00596 }
00597 return ILP;
00598 }
00599
00600
00601 G4double G4ChipsKaonMinusElasticXS::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==310 || PDG==130) PDG=-321;
00608 if(PDG!=-321)G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetET:PDG="<<PDG<<G4endl;
00609 if(onlyCS) G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetExT: onlyCS=1"<<G4endl;
00610 if(lastLP<-4.3) return lastTM*GeVSQ*G4UniformRand();
00611 G4double q2=0.;
00612 if(tgZ==1 && tgN==0)
00613 {
00614 G4double E1=lastTM*theB1;
00615 G4double R1=(1.-std::exp(-E1));
00616 G4double E2=lastTM*theB2;
00617 G4double R2=(1.-std::exp(-E2*E2*E2));
00618 G4double E3=lastTM*theB3;
00619 G4double R3=(1.-std::exp(-E3));
00620 G4double I1=R1*theS1/theB1;
00621 G4double I2=R2*theS2;
00622 G4double I3=R3*theS3;
00623 G4double I12=I1+I2;
00624 G4double rand=(I12+I3)*G4UniformRand();
00625 if (rand<I1 )
00626 {
00627 G4double ran=R1*G4UniformRand();
00628 if(ran>1.) ran=1.;
00629 q2=-std::log(1.-ran)/theB1;
00630 }
00631 else if(rand<I12)
00632 {
00633 G4double ran=R2*G4UniformRand();
00634 if(ran>1.) ran=1.;
00635 q2=-std::log(1.-ran);
00636 if(q2<0.) q2=0.;
00637 q2=std::pow(q2,third)/theB2;
00638 }
00639 else
00640 {
00641 G4double ran=R3*G4UniformRand();
00642 if(ran>1.) ran=1.;
00643 q2=-std::log(1.-ran)/theB3;
00644 }
00645 }
00646 else
00647 {
00648 G4double a=tgZ+tgN;
00649 G4double E1=lastTM*(theB1+lastTM*theSS);
00650 G4double R1=(1.-std::exp(-E1));
00651 G4double tss=theSS+theSS;
00652 G4double tm2=lastTM*lastTM;
00653 G4double E2=lastTM*tm2*theB2;
00654 if(a>6.5)E2*=tm2;
00655 G4double R2=(1.-std::exp(-E2));
00656 G4double E3=lastTM*theB3;
00657 if(a>6.5)E3*=tm2*tm2*tm2;
00658 G4double R3=(1.-std::exp(-E3));
00659 G4double E4=lastTM*theB4;
00660 G4double R4=(1.-std::exp(-E4));
00661 G4double I1=R1*theS1;
00662 G4double I2=R2*theS2;
00663 G4double I3=R3*theS3;
00664 G4double I4=R4*theS4;
00665 G4double I12=I1+I2;
00666 G4double I13=I12+I3;
00667 G4double rand=(I13+I4)*G4UniformRand();
00668 if(rand<I1)
00669 {
00670 G4double ran=R1*G4UniformRand();
00671 if(ran>1.) ran=1.;
00672 q2=-std::log(1.-ran)/theB1;
00673 if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
00674 }
00675 else if(rand<I12)
00676 {
00677 G4double ran=R2*G4UniformRand();
00678 if(ran>1.) ran=1.;
00679 q2=-std::log(1.-ran)/theB2;
00680 if(q2<0.) q2=0.;
00681 if(a<6.5) q2=std::pow(q2,third);
00682 else q2=std::pow(q2,fifth);
00683 }
00684 else if(rand<I13)
00685 {
00686 G4double ran=R3*G4UniformRand();
00687 if(ran>1.) ran=1.;
00688 q2=-std::log(1.-ran)/theB3;
00689 if(q2<0.) q2=0.;
00690 if(a>6.5) q2=std::pow(q2,sevth);
00691 }
00692 else
00693 {
00694 G4double ran=R4*G4UniformRand();
00695 if(ran>1.) ran=1.;
00696 q2=-std::log(1.-ran)/theB4;
00697 if(a<6.5) q2=lastTM-q2;
00698 }
00699 }
00700 if(q2<0.) q2=0.;
00701 if(!(q2>=-1.||q2<=1.)) G4cout<<"*NAN*G4QKaonMinusElasticCS::GetExchT: -t="<<q2<<G4endl;
00702 if(q2>lastTM)
00703 {
00704 q2=lastTM;
00705 }
00706 return q2*GeVSQ;
00707 }
00708
00709
00710 G4double G4ChipsKaonMinusElasticXS::GetSlope(G4int tgZ, G4int tgN, G4int PDG)
00711 {
00712 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
00713 if(onlyCS)G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetSl:onlCS=true"<<G4endl;
00714 if(lastLP<-4.3) return 0.;
00715 if(PDG != -321)
00716 {
00717
00718
00719
00720 G4ExceptionDescription ed;
00721 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
00722 << ", while it is defined only for PDG=-321 (K-)" << G4endl;
00723 }
00724 if(theB1<0.) theB1=0.;
00725 if(!(theB1>=-1.||theB1<=1.))G4cout<<"*NAN*G4QKaonMinusElCS::GetSlope:B1="<<theB1<<G4endl;
00726 return theB1/GeVSQ;
00727 }
00728
00729
00730 G4double G4ChipsKaonMinusElasticXS::GetHMaxT()
00731 {
00732 static const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
00733 return lastTM*HGeVSQ;
00734 }
00735
00736
00737 G4double G4ChipsKaonMinusElasticXS::GetTabValues(G4double lp, G4int PDG, G4int tgZ,
00738 G4int tgN)
00739 {
00740 if(PDG!=-321)G4cout<<"*Warning*G4ChipsKaonMinusElasticXS::GetTV:PDG="<<PDG<<G4endl;
00741 if(tgZ<0 || tgZ>92)
00742 {
00743 G4cout<<"*Warning*G4QKaonMinusElasticCS::GetTabV:(1-92)NoIsotopes for Z="<<tgZ<<G4endl;
00744 return 0.;
00745 }
00746 G4int iZ=tgZ-1;
00747 if(iZ<0)
00748 {
00749 iZ=0;
00750 tgZ=1;
00751 tgN=0;
00752 }
00753 G4double p=std::exp(lp);
00754 G4double sp=std::sqrt(p);
00755 G4double psp=p*sp;
00756 G4double p2=p*p;
00757 G4double p3=p2*p;
00758 G4double p4=p3*p;
00759 if ( tgZ == 1 && tgN == 0 )
00760 {
00761 G4double dl2=lp-lastPAR[12];
00762 theSS=lastPAR[35];
00763 theS1=(lastPAR[13]+lastPAR[14]*dl2*dl2)/(1.+lastPAR[15]/p4/p)+
00764 (lastPAR[16]/p2+lastPAR[17]*p)/(p4+lastPAR[18]*sp);
00765 theB1=lastPAR[19]*std::pow(p,lastPAR[20])/(1.+lastPAR[21]/p3);
00766 theS2=lastPAR[22]+lastPAR[23]/(p4+lastPAR[24]*p);
00767 theB2=lastPAR[25]+lastPAR[26]/(p4+lastPAR[27]/sp);
00768 theS3=lastPAR[28]+lastPAR[29]/(p4*p4+lastPAR[30]*p2+lastPAR[31]);
00769 theB3=lastPAR[32]+lastPAR[33]/(p4+lastPAR[34]);
00770 theS4=0.;
00771 theB4=0.;
00772
00773 G4double dp=lp-lastPAR[2];
00774 return lastPAR[0]/psp+(lastPAR[1]*dp*dp+lastPAR[3])/(1.-lastPAR[4]/sp+lastPAR[5]/p4)+
00775 lastPAR[6]/(sqr(p-lastPAR[7])+lastPAR[8])+lastPAR[9]/(sqr(p-lastPAR[10])+lastPAR[11]);
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]/p3)/(1.+lastPAR[3]/p2/sp);
00825 }
00826 return 0.;
00827 }
00828
00829
00830 G4double G4ChipsKaonMinusElasticXS::GetQ2max(G4int PDG, G4int tgZ, G4int tgN,
00831 G4double pP)
00832 {
00833 static const G4double mK= G4KaonMinus::KaonMinus()->GetPDGMass()*.001;
00834
00835 static const G4double mK2= mK*mK;
00836
00837 G4double pP2=pP*pP;
00838 if(tgZ || tgN>-1)
00839 {
00840 G4double mt=G4ParticleTable::GetParticleTable()->FindIon(tgZ,tgZ+tgN,0,tgZ)->GetPDGMass()*.001;
00841
00842 G4double dmt=mt+mt;
00843 G4double mds=dmt*std::sqrt(pP2+mK2)+mK2+mt*mt;
00844 return dmt*dmt*pP2/mds;
00845 }
00846 else
00847 {
00848 G4ExceptionDescription ed;
00849 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
00850 << ", while it is defined only for p projectiles & Z_target>0" << G4endl;
00851 G4Exception("G4ChipsKaonMinusElasticXS::GetQ2max()", "HAD_CHPS_0000",
00852 FatalException, ed);
00853 return 0;
00854 }
00855 }