00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #include "G4QKaonMinusElasticCrossSection.hh"
00046 #include "G4SystemOfUnits.hh"
00047
00048
00049 const G4int G4QKaonMinusElasticCrossSection::nPoints=128;
00050 const G4int G4QKaonMinusElasticCrossSection::nLast=nPoints-1;
00051 G4double G4QKaonMinusElasticCrossSection::lPMin=-8.;
00052 G4double G4QKaonMinusElasticCrossSection::lPMax= 8.;
00053 G4double G4QKaonMinusElasticCrossSection::dlnP=(lPMax-lPMin)/nLast;
00054 G4bool G4QKaonMinusElasticCrossSection::onlyCS=true;
00055 G4double G4QKaonMinusElasticCrossSection::lastSIG=0.;
00056 G4double G4QKaonMinusElasticCrossSection::lastLP=-10.;
00057 G4double G4QKaonMinusElasticCrossSection::lastTM=0.;
00058 G4double G4QKaonMinusElasticCrossSection::theSS=0.;
00059 G4double G4QKaonMinusElasticCrossSection::theS1=0.;
00060 G4double G4QKaonMinusElasticCrossSection::theB1=0.;
00061 G4double G4QKaonMinusElasticCrossSection::theS2=0.;
00062 G4double G4QKaonMinusElasticCrossSection::theB2=0.;
00063 G4double G4QKaonMinusElasticCrossSection::theS3=0.;
00064 G4double G4QKaonMinusElasticCrossSection::theB3=0.;
00065 G4double G4QKaonMinusElasticCrossSection::theS4=0.;
00066 G4double G4QKaonMinusElasticCrossSection::theB4=0.;
00067 G4int G4QKaonMinusElasticCrossSection::lastTZ=0;
00068 G4int G4QKaonMinusElasticCrossSection::lastTN=0;
00069 G4double G4QKaonMinusElasticCrossSection::lastPIN=0.;
00070 G4double* G4QKaonMinusElasticCrossSection::lastCST=0;
00071 G4double* G4QKaonMinusElasticCrossSection::lastPAR=0;
00072 G4double* G4QKaonMinusElasticCrossSection::lastSST=0;
00073 G4double* G4QKaonMinusElasticCrossSection::lastS1T=0;
00074 G4double* G4QKaonMinusElasticCrossSection::lastB1T=0;
00075 G4double* G4QKaonMinusElasticCrossSection::lastS2T=0;
00076 G4double* G4QKaonMinusElasticCrossSection::lastB2T=0;
00077 G4double* G4QKaonMinusElasticCrossSection::lastS3T=0;
00078 G4double* G4QKaonMinusElasticCrossSection::lastB3T=0;
00079 G4double* G4QKaonMinusElasticCrossSection::lastS4T=0;
00080 G4double* G4QKaonMinusElasticCrossSection::lastB4T=0;
00081 G4int G4QKaonMinusElasticCrossSection::lastN=0;
00082 G4int G4QKaonMinusElasticCrossSection::lastZ=0;
00083 G4double G4QKaonMinusElasticCrossSection::lastP=0.;
00084 G4double G4QKaonMinusElasticCrossSection::lastTH=0.;
00085 G4double G4QKaonMinusElasticCrossSection::lastCS=0.;
00086 G4int G4QKaonMinusElasticCrossSection::lastI=0;
00087
00088 std::vector<G4double*> G4QKaonMinusElasticCrossSection::PAR;
00089 std::vector<G4double*> G4QKaonMinusElasticCrossSection::CST;
00090 std::vector<G4double*> G4QKaonMinusElasticCrossSection::SST;
00091 std::vector<G4double*> G4QKaonMinusElasticCrossSection::S1T;
00092 std::vector<G4double*> G4QKaonMinusElasticCrossSection::B1T;
00093 std::vector<G4double*> G4QKaonMinusElasticCrossSection::S2T;
00094 std::vector<G4double*> G4QKaonMinusElasticCrossSection::B2T;
00095 std::vector<G4double*> G4QKaonMinusElasticCrossSection::S3T;
00096 std::vector<G4double*> G4QKaonMinusElasticCrossSection::B3T;
00097 std::vector<G4double*> G4QKaonMinusElasticCrossSection::S4T;
00098 std::vector<G4double*> G4QKaonMinusElasticCrossSection::B4T;
00099
00100 G4QKaonMinusElasticCrossSection::G4QKaonMinusElasticCrossSection()
00101 {
00102 }
00103
00104 G4QKaonMinusElasticCrossSection::~G4QKaonMinusElasticCrossSection()
00105 {
00106 std::vector<G4double*>::iterator pos;
00107 for (pos=CST.begin(); pos<CST.end(); pos++)
00108 { delete [] *pos; }
00109 CST.clear();
00110 for (pos=PAR.begin(); pos<PAR.end(); pos++)
00111 { delete [] *pos; }
00112 PAR.clear();
00113 for (pos=SST.begin(); pos<SST.end(); pos++)
00114 { delete [] *pos; }
00115 SST.clear();
00116 for (pos=S1T.begin(); pos<S1T.end(); pos++)
00117 { delete [] *pos; }
00118 S1T.clear();
00119 for (pos=B1T.begin(); pos<B1T.end(); pos++)
00120 { delete [] *pos; }
00121 B1T.clear();
00122 for (pos=S2T.begin(); pos<S2T.end(); pos++)
00123 { delete [] *pos; }
00124 S2T.clear();
00125 for (pos=B2T.begin(); pos<B2T.end(); pos++)
00126 { delete [] *pos; }
00127 B2T.clear();
00128 for (pos=S3T.begin(); pos<S3T.end(); pos++)
00129 { delete [] *pos; }
00130 S3T.clear();
00131 for (pos=B3T.begin(); pos<B3T.end(); pos++)
00132 { delete [] *pos; }
00133 B3T.clear();
00134 for (pos=S4T.begin(); pos<S4T.end(); pos++)
00135 { delete [] *pos; }
00136 S4T.clear();
00137 for (pos=B4T.begin(); pos<B4T.end(); pos++)
00138 { delete [] *pos; }
00139 B4T.clear();
00140 }
00141
00142
00143 G4VQCrossSection* G4QKaonMinusElasticCrossSection::GetPointer()
00144 {
00145 static G4QKaonMinusElasticCrossSection theCrossSection;
00146 return &theCrossSection;
00147 }
00148
00149
00150
00151 G4double G4QKaonMinusElasticCrossSection::GetCrossSection(G4bool fCS, G4double pMom,
00152 G4int tgZ, G4int tgN, G4int pPDG)
00153 {
00154 static std::vector <G4int> colN;
00155 static std::vector <G4int> colZ;
00156 static std::vector <G4double> colP;
00157 static std::vector <G4double> colTH;
00158 static std::vector <G4double> colCS;
00159
00160 G4double pEn=pMom;
00161 onlyCS=fCS;
00162 #ifdef pdebug
00163 G4cout<<"G4QKaMiusElCS::GetCS:>>>f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="
00164 <<tgN<<"("<<lastN<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="<<colN.size()<<G4endl;
00165
00166 #endif
00167 if(pPDG != -321)
00168 {
00169 G4cout<<"*Warning*G4QKaonMinusElaCS::GetCS:*>Found pPDG="<<pPDG<<" ==> CS=0"<<G4endl;
00170
00171 return 0.;
00172 }
00173 G4bool in=false;
00174 lastP = 0.;
00175 lastN = tgN;
00176 lastZ = tgZ;
00177 lastI = colN.size();
00178 if(lastI) for(G4int i=0; i<lastI; i++)
00179 {
00180 if(colN[i]==tgN && colZ[i]==tgZ)
00181 {
00182 lastI=i;
00183 lastTH =colTH[i];
00184 #ifdef pdebug
00185 G4cout<<"G4QKMElCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",i="<<i<<G4endl;
00186
00187 #endif
00188 if(pEn<=lastTH)
00189 {
00190 #ifdef pdebug
00191 G4cout<<"G4QKMElCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",CS=0"<<G4endl;
00192
00193 #endif
00194 return 0.;
00195 }
00196 lastP =colP [i];
00197 lastCS =colCS[i];
00198
00199 if(lastP == pMom)
00200 {
00201 #ifdef pdebug
00202 G4cout<<"G4QKaonMinusElastCS::GetCS: P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00203 #endif
00204 CalculateCrossSection(fCS,-1,i,pPDG,lastZ,lastN,pMom);
00205 return lastCS*millibarn;
00206 }
00207 in = true;
00208
00209 #ifdef pdebug
00210 G4cout<<"G4QKMElCS::G:UpdateDB P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",i="<<i<<G4endl;
00211 #endif
00212 lastCS=CalculateCrossSection(fCS,-1,i,pPDG,lastZ,lastN,pMom);
00213 #ifdef pdebug
00214 G4cout<<"G4QKaonMinusElCS::GetCrosSec:****>New(inDB) Calculated CS="<<lastCS<<G4endl;
00215
00216 #endif
00217 if(lastCS<=0. && pEn>lastTH)
00218 {
00219 #ifdef pdebug
00220 G4cout<<"G4QKaonMinusElCS::GetCS: T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00221 #endif
00222 lastTH=pEn;
00223 }
00224 break;
00225 }
00226 #ifdef pdebug
00227 G4cout<<"---G4QKaonMinusElasticCrossSection::GetCrosSec:pPDG="<<pPDG<<",i="<<i<<",N="
00228 <<colN[i]<<",Z["<<i<<"]="<<colZ[i]<<G4endl;
00229
00230 #endif
00231 }
00232 if(!in)
00233 {
00234 #ifdef pdebug
00235 G4cout<<"G4QKaonMinusElCS::GetCS:CalNewP="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
00236 #endif
00238 lastCS=CalculateCrossSection(fCS,0,lastI,pPDG,lastZ,lastN,pMom);//calculate&create
00239 if(lastCS<=0.)
00240 {
00241 lastTH = ThresholdEnergy(tgZ, tgN);
00242 #ifdef pdebug
00243 G4cout<<"G4QKaonMinusElasticCrossSection::GetCS:NewThr="<<lastTH<<",T="<<pEn<<G4endl;
00244 #endif
00245 if(pEn>lastTH)
00246 {
00247 #ifdef pdebug
00248 G4cout<<"G4QKaonMinusElCS::GetCS:1st T="<<pEn<<"(CS=0)>Threshold="<<lastTH<<G4endl;
00249 #endif
00250 lastTH=pEn;
00251 }
00252 }
00253 #ifdef pdebug
00254 G4cout<<"G4QKaonMinusElCS::GetCS:NCS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
00255
00256 #endif
00257 colN.push_back(tgN);
00258 colZ.push_back(tgZ);
00259 colP.push_back(pMom);
00260 colTH.push_back(lastTH);
00261 colCS.push_back(lastCS);
00262 #ifdef pdebug
00263 G4cout<<"G4QKMElCS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00264
00265 #endif
00266 return lastCS*millibarn;
00267 }
00268 else
00269 {
00270 #ifdef pdebug
00271 G4cout<<"G4QKaonMinusElasticCrossSection::GetCS: Update lastI="<<lastI<<G4endl;
00272 #endif
00273 colP[lastI]=pMom;
00274 colCS[lastI]=lastCS;
00275 }
00276 #ifdef pdebug
00277 G4cout<<"G4QKMElCS::GetCSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00278
00279 G4cout<<"G4QKaonMinusElasticCrossSection::GetCrosSec:**End**, onlyCS="<<onlyCS<<G4endl;
00280 #endif
00281 return lastCS*millibarn;
00282 }
00283
00284
00285
00286 G4double G4QKaonMinusElasticCrossSection::CalculateCrossSection(G4bool CS, G4int F,
00287 G4int I, G4int PDG, G4int tgZ, G4int tgN, G4double pIU)
00288 {
00289
00290 static std::vector <G4double> PIN;
00291
00292 G4double pMom=pIU/GeV;
00293 onlyCS=CS;
00294 #ifdef pdebug
00295 G4cout<<"G4QKaonMinusElasticCS::CalcCS: onlyCS="<<onlyCS<<",F="<<F<<",p="<<pIU<<G4endl;
00296 #endif
00297 lastLP=std::log(pMom);
00298 if(F)
00299 {
00300 if(F<0)
00301 {
00302 lastPIN = PIN[I];
00303 lastPAR = PAR[I];
00304 lastCST = CST[I];
00305 lastSST = SST[I];
00306 lastS1T = S1T[I];
00307 lastB1T = B1T[I];
00308 lastS2T = S2T[I];
00309 lastB2T = B2T[I];
00310 lastS3T = S3T[I];
00311 lastB3T = B3T[I];
00312 lastS4T = S4T[I];
00313 lastB4T = B4T[I];
00314 #ifdef pdebug
00315 G4cout<<"G4QKaonMinusElCS::CalCS:DB updated for I="<<I<<",*,PIN4="<<PIN[4]<<G4endl;
00316 #endif
00317 }
00318 #ifdef pdebug
00319 G4cout<<"G4QKaonMinusElasticCS::CalcCS:*read*,LP="<<lastLP<<",PIN="<<lastPIN<<G4endl;
00320 #endif
00321 if(lastLP>lastPIN && lastLP<lPMax)
00322 {
00323 lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
00324 #ifdef pdebug
00325 G4cout<<"G4QKMElCS::CalcCS:updated(I),LP="<<lastLP<<"<IN["<<I<<"]="<<lastPIN<<G4endl;
00326 #endif
00327 PIN[I]=lastPIN;
00328 }
00329 }
00330 else
00331 {
00332 lastPAR = new G4double[nPoints];
00333 lastPAR[nLast]=0;
00334 lastCST = new G4double[nPoints];
00335 lastSST = new G4double[nPoints];
00336 lastS1T = new G4double[nPoints];
00337 lastB1T = new G4double[nPoints];
00338 lastS2T = new G4double[nPoints];
00339 lastB2T = new G4double[nPoints];
00340 lastS3T = new G4double[nPoints];
00341 lastB3T = new G4double[nPoints];
00342 lastS4T = new G4double[nPoints];
00343 lastB4T = new G4double[nPoints];
00344 #ifdef pdebug
00345 G4cout<<"G4QKaonMinusElasticCS::CalcCS:*ini*,lstLP="<<lastLP<<",min="<<lPMin<<G4endl;
00346 #endif
00347 lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN);
00348 #ifdef pdebug
00349 G4cout<<"G4QKaMiElCS::CCS:i,Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<",LP"<<lastPIN<<G4endl;
00350 #endif
00351 PIN.push_back(lastPIN);
00352 PAR.push_back(lastPAR);
00353 CST.push_back(lastCST);
00354 SST.push_back(lastSST);
00355 S1T.push_back(lastS1T);
00356 B1T.push_back(lastB1T);
00357 S2T.push_back(lastS2T);
00358 B2T.push_back(lastB2T);
00359 S3T.push_back(lastS3T);
00360 B3T.push_back(lastB3T);
00361 S4T.push_back(lastS4T);
00362 B4T.push_back(lastB4T);
00363 }
00364
00365 #ifdef pdebug
00366 G4cout<<"G4QKMElCS::CalcCS:?update?,LP="<<lastLP<<",IN="<<lastPIN<<",ML="<<lPMax<<G4endl;
00367 #endif
00368 if(lastLP>lastPIN && lastLP<lPMax)
00369 {
00370 lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
00371 #ifdef pdebug
00372 G4cout<<"G4QKaonMinusElCS::CalcCS:*updated(O)*,LP="<<lastLP<<"<IN="<<lastPIN<<G4endl;
00373 #endif
00374 }
00375 #ifdef pdebug
00376 G4cout<<"G4QKaMiElCS::CalcCS:lastLP="<<lastLP<<",lPM="<<lPMin<<",lPIN="<<lastPIN<<G4endl;
00377 #endif
00378 if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom);
00379 #ifdef pdebug
00380 G4cout<<"G4QKaMiElCrosSec::CalcCS: oCS="<<onlyCS<<",-t="<<lastTM<<", p="<<lastLP<<G4endl;
00381 #endif
00382 if(lastLP>lPMin && lastLP<=lastPIN)
00383 {
00384 if(lastLP==lastPIN)
00385 {
00386 G4double shift=(lastLP-lPMin)/dlnP+.000001;
00387 G4int blast=static_cast<int>(shift);
00388 if(blast<0 || blast>=nLast) G4cout<<"G4QKMElCS::CCS:b="<<blast<<",n="<<nLast<<G4endl;
00389 lastSIG = lastCST[blast];
00390 if(!onlyCS)
00391 {
00392 theSS = lastSST[blast];
00393 theS1 = lastS1T[blast];
00394 theB1 = lastB1T[blast];
00395 theS2 = lastS2T[blast];
00396 theB2 = lastB2T[blast];
00397 theS3 = lastS3T[blast];
00398 theB3 = lastB3T[blast];
00399 theS4 = lastS4T[blast];
00400 theB4 = lastB4T[blast];
00401 }
00402 #ifdef pdebug
00403 G4cout<<"G4QKaonMinusElasticCS::CalculateCS:(E) S1="<<theS1<<",B1="<<theB1<<G4endl;
00404 #endif
00405 }
00406 else
00407 {
00408 G4double shift=(lastLP-lPMin)/dlnP;
00409 G4int blast=static_cast<int>(shift);
00410 if(blast<0) blast=0;
00411 if(blast>=nLast) blast=nLast-1;
00412 shift-=blast;
00413 G4int lastL=blast+1;
00414 G4double SIGL=lastCST[blast];
00415 lastSIG= SIGL+shift*(lastCST[lastL]-SIGL);
00416 #ifdef pdebug
00417 G4cout<<"G4QKaonMinusElasticCrossSection::CalcCrossSection:Sig="<<lastSIG<<",P="
00418 <<pMom<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<",onlyCS="<<onlyCS<<G4endl;
00419 #endif
00420 if(!onlyCS)
00421 {
00422 G4double SSTL=lastSST[blast];
00423 theSS=SSTL+shift*(lastSST[lastL]-SSTL);
00424 G4double S1TL=lastS1T[blast];
00425 theS1=S1TL+shift*(lastS1T[lastL]-S1TL);
00426 G4double B1TL=lastB1T[blast];
00427 #ifdef pdebug
00428 G4cout<<"G4QKaonMinusElasticCS::CalcCrossSect: b="<<blast<<", ls="<<lastL<<",SL="
00429 <<S1TL<<",SU="<<lastS1T[lastL]<<",BL="<<B1TL<<",BU="<<lastB1T[lastL]<<G4endl;
00430 #endif
00431 theB1=B1TL+shift*(lastB1T[lastL]-B1TL);
00432 G4double S2TL=lastS2T[blast];
00433 theS2=S2TL+shift*(lastS2T[lastL]-S2TL);
00434 G4double B2TL=lastB2T[blast];
00435 theB2=B2TL+shift*(lastB2T[lastL]-B2TL);
00436 G4double S3TL=lastS3T[blast];
00437 theS3=S3TL+shift*(lastS3T[lastL]-S3TL);
00438 #ifdef pdebug
00439 G4cout<<"G4QKaonMinusElasticCrossSection::CalcCS: s3l="<<S3TL<<", sh3="<<shift
00440 <<", s3h="<<lastS3T[lastL]<<", b="<<blast<<", l="<<lastL<<G4endl;
00441 #endif
00442 G4double B3TL=lastB3T[blast];
00443 theB3=B3TL+shift*(lastB3T[lastL]-B3TL);
00444 G4double S4TL=lastS4T[blast];
00445 theS4=S4TL+shift*(lastS4T[lastL]-S4TL);
00446 #ifdef pdebug
00447 G4cout<<"G4QKaonMinusElasticCrossSection::CalcCS: s4l="<<S4TL<<", sh4="<<shift
00448 <<",s4h="<<lastS4T[lastL]<<",b="<<blast<<",l="<<lastL<<G4endl;
00449 #endif
00450 G4double B4TL=lastB4T[blast];
00451 theB4=B4TL+shift*(lastB4T[lastL]-B4TL);
00452 }
00453 #ifdef pdebug
00454 G4cout<<"G4QKaonMinusElasticCS::CalculateCS:(I) S1="<<theS1<<",B1="<<theB1<<G4endl;
00455 #endif
00456 }
00457 }
00458 else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN);
00459 if(lastSIG<0.) lastSIG = 0.;
00460 #ifdef pdebug
00461 G4cout<<"G4QKaonMinusElasticCrossSection::CalculateCS: END, onlyCS="<<onlyCS<<G4endl;
00462 #endif
00463 return lastSIG;
00464 }
00465
00466
00467 G4double G4QKaonMinusElasticCrossSection::GetPTables(G4double LP,G4double ILP, G4int PDG,
00468 G4int tgZ, G4int tgN)
00469 {
00470
00471 static const G4double pwd=2727;
00472 const G4int n_kmpel=36;
00473
00474 G4double kmp_el[n_kmpel]={5.2,.0557,3.5,2.23,.7,.075,.004,.39,.000156,.15,1.,.0156,5.,
00475 74.,3.,3.4,.2,.17,.001,8.,.055,3.64,5.e-5,4000.,1500.,.46,
00476 1.2e6,3.5e6,5.e-5,1.e10,8.5e8,1.e10,1.1,3.4e6,6.8e6,0.};
00477
00478
00479 if(PDG == -321)
00480 {
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492 if(lastPAR[nLast]!=pwd)
00493 {
00494 if ( tgZ == 1 && tgN == 0 )
00495 {
00496 for (G4int ip=0; ip<n_kmpel; ip++) lastPAR[ip]=kmp_el[ip];
00497 }
00498 else
00499 {
00500 G4double a=tgZ+tgN;
00501 G4double sa=std::sqrt(a);
00502 G4double ssa=std::sqrt(sa);
00503 G4double asa=a*sa;
00504 G4double a2=a*a;
00505 G4double a3=a2*a;
00506 G4double a4=a3*a;
00507 G4double a5=a4*a;
00508 G4double a6=a4*a2;
00509 G4double a7=a6*a;
00510 G4double a8=a7*a;
00511 G4double a9=a8*a;
00512 G4double a10=a5*a5;
00513 G4double a12=a6*a6;
00514 G4double a14=a7*a7;
00515 G4double a16=a8*a8;
00516 G4double a17=a16*a;
00517
00518 G4double a32=a16*a16;
00519
00520 lastPAR[0]=.06*asa/(1.+a*(.01+.1/ssa));
00521 lastPAR[1]=.75*asa/(1.+.009*a);
00522 lastPAR[2]=.1*a2*ssa/(1.+.0015*a2/ssa);
00523 lastPAR[3]=1./(1.+500./a2);
00524 lastPAR[4]=4.2;
00525 lastPAR[5]=0.;
00526 lastPAR[6]=0.;
00527 lastPAR[7]=0.;
00528 lastPAR[8]=0.;
00529
00530 if(a<6.5)
00531 {
00532 G4double a28=a16*a12;
00533
00534 lastPAR[ 9]=4000*a;
00535 lastPAR[10]=1.2e7*a8+380*a17;
00536 lastPAR[11]=.7/(1.+4.e-12*a16);
00537 lastPAR[12]=2.5/a8/(a4+1.e-16*a32);
00538 lastPAR[13]=.28*a;
00539 lastPAR[14]=1.2*a2+2.3;
00540 lastPAR[15]=3.8/a;
00541
00542 lastPAR[16]=.01/(1.+.0024*a5);
00543 lastPAR[17]=.2*a;
00544 lastPAR[18]=9.e-7/(1.+.035*a5);
00545 lastPAR[19]=(42.+2.7e-11*a16)/(1.+.14*a);
00546
00547 lastPAR[20]=2.25*a3;
00548 lastPAR[21]=18.;
00549 lastPAR[22]=2.4e-3*a8/(1.+2.6e-4*a7);
00550 lastPAR[23]=3.5e-36*a32*a8/(1.+5.e-15*a32/a);
00551
00552 lastPAR[24]=1.e5/(a8+2.5e12/a16);
00553 lastPAR[25]=8.e7/(a12+1.e-27*a28*a28);
00554 lastPAR[26]=.0006*a3;
00555
00556 lastPAR[27]=10.+4.e-8*a12*a;
00557 lastPAR[28]=.114;
00558 lastPAR[29]=.003;
00559 lastPAR[30]=2.e-23;
00560
00561 lastPAR[31]=1./(1.+.0001*a8);
00562 lastPAR[32]=1.5e-4/(1.+5.e-6*a12);
00563 lastPAR[33]=.03;
00564
00565 lastPAR[34]=a/2;
00566 lastPAR[35]=2.e-7*a4;
00567 lastPAR[36]=4.;
00568 lastPAR[37]=64./a3;
00569
00570 lastPAR[38]=1.e8*std::exp(.32*asa);
00571 lastPAR[39]=20.*std::exp(.45*asa);
00572 lastPAR[40]=7.e3+2.4e6/a5;
00573 lastPAR[41]=2.5e5*std::exp(.085*a3);
00574 lastPAR[42]=2.5*a;
00575
00576 lastPAR[43]=920.+.03*a8*a3;
00577 lastPAR[44]=93.+.0023*a12;
00578 #ifdef pdebug
00579 G4cout<<"G4QKMElCS::CalcCS:la "<<lastPAR[38]<<", "<<lastPAR[39]<<", "<<lastPAR[40]
00580 <<", "<<lastPAR[42]<<", "<<lastPAR[43]<<", "<<lastPAR[44]<<G4endl;
00581 #endif
00582 }
00583 else
00584 {
00585 G4double p1a10=2.2e-28*a10;
00586 G4double r4a16=6.e14/a16;
00587 G4double s4a16=r4a16*r4a16;
00588
00589
00590
00591 lastPAR[ 9]=4.5*std::pow(a,1.15);
00592 lastPAR[10]=.06*std::pow(a,.6);
00593 lastPAR[11]=.6*a/(1.+2.e15/a16);
00594 lastPAR[12]=.17/(a+9.e5/a3+1.5e33/a32);
00595 lastPAR[13]=(.001+7.e-11*a5)/(1.+4.4e-11*a5);
00596 lastPAR[14]=(p1a10*p1a10+2.e-29)/(1.+2.e-22*a12);
00597
00598 lastPAR[15]=400./a12+2.e-22*a9;
00599 lastPAR[16]=1.e-32*a12/(1.+5.e22/a14);
00600 lastPAR[17]=1000./a2+9.5*sa*ssa;
00601 lastPAR[18]=4.e-6*a*asa+1.e11/a16;
00602 lastPAR[19]=(120./a+.002*a2)/(1.+2.e14/a16);
00603 lastPAR[20]=9.+100./a;
00604
00605 lastPAR[21]=.002*a3+3.e7/a6;
00606 lastPAR[22]=7.e-15*a4*asa;
00607 lastPAR[23]=9000./a4;
00608
00609 lastPAR[24]=.0011*asa/(1.+3.e34/a32/a4);
00610 lastPAR[25]=1.e-5*a2+2.e14/a16;
00611 lastPAR[26]=1.2e-11*a2/(1.+1.5e19/a12);
00612 lastPAR[27]=.016*asa/(1.+5.e16/a16);
00613
00614 lastPAR[28]=.002*a4/(1.+7.e7/std::pow(a-6.83,14));
00615 lastPAR[29]=2.e6/a6+7.2/std::pow(a,.11);
00616 lastPAR[30]=11.*a3/(1.+7.e23/a16/a8);
00617 lastPAR[31]=100./asa;
00618
00619 lastPAR[32]=(.1+4.4e-5*a2)/(1.+5.e5/a4);
00620 lastPAR[33]=3.5e-4*a2/(1.+1.e8/a8);
00621 lastPAR[34]=1.3+3.e5/a4;
00622 lastPAR[35]=500./(a2+50.)+3;
00623 lastPAR[36]=1.e-9/a+s4a16*s4a16;
00624
00625 lastPAR[37]=.4*asa+3.e-9*a6;
00626 lastPAR[38]=.0005*a5;
00627 lastPAR[39]=.002*a5;
00628 lastPAR[40]=10.;
00629
00630 lastPAR[41]=.05+.005*a;
00631 lastPAR[42]=7.e-8/sa;
00632 lastPAR[43]=.8*sa;
00633 lastPAR[44]=.02*sa;
00634 lastPAR[45]=1.e8/a3;
00635 lastPAR[46]=3.e32/(a32+1.e32);
00636
00637 lastPAR[47]=24.;
00638 lastPAR[48]=20./sa;
00639 lastPAR[49]=7.e3*a/(sa+1.);
00640 lastPAR[50]=900.*sa/(1.+500./a3);
00641 #ifdef pdebug
00642 G4cout<<"G4QKMElCS::CalcCS:ha "<<lastPAR[41]<<", "<<lastPAR[42]<<", "<<lastPAR[43]
00643 <<", "<<lastPAR[44]<<", "<<lastPAR[45]<<", "<<lastPAR[46]<<G4endl;
00644 #endif
00645 }
00646
00647 lastPAR[51]=1.e15+2.e27/a4/(1.+2.e-18*a16);
00648 }
00649 lastPAR[nLast]=pwd;
00650
00651 G4double lp=lPMin;
00652 G4bool memCS=onlyCS;
00653 onlyCS=false;
00654 lastCST[0]=GetTabValues(lp, PDG, tgZ, tgN);
00655 onlyCS=memCS;
00656 lastSST[0]=theSS;
00657 lastS1T[0]=theS1;
00658 lastB1T[0]=theB1;
00659 lastS2T[0]=theS2;
00660 lastB2T[0]=theB2;
00661 lastS3T[0]=theS3;
00662 lastB3T[0]=theB3;
00663 lastS4T[0]=theS4;
00664 lastB4T[0]=theB4;
00665 #ifdef pdebug
00666 G4cout<<"G4QKaonMinusElasticCrossSection::GetPTables:ip=0(init), lp="<<lp<<",S1="
00667 <<theS1<<",B1="<<theB1<<",S2="<<theS2<<",B2="<<theB3<<",S3="<<theS3
00668 <<",B3="<<theB3<<",S4="<<theS4<<",B4="<<theB4<<G4endl;
00669 #endif
00670 }
00671 if(LP>ILP)
00672 {
00673 G4int ini = static_cast<int>((ILP-lPMin+.000001)/dlnP)+1;
00674 if(ini<0) ini=0;
00675 if(ini<nPoints)
00676 {
00677 G4int fin = static_cast<int>((LP-lPMin)/dlnP)+1;
00678 if(fin>=nPoints) fin=nLast;
00679 if(fin>=ini)
00680 {
00681 G4double lp=0.;
00682 for(G4int ip=ini; ip<=fin; ip++)
00683 {
00684 lp=lPMin+ip*dlnP;
00685 G4bool memCS=onlyCS;
00686 onlyCS=false;
00687 lastCST[ip]=GetTabValues(lp, PDG, tgZ, tgN);
00688 onlyCS=memCS;
00689 lastSST[ip]=theSS;
00690 lastS1T[ip]=theS1;
00691 lastB1T[ip]=theB1;
00692 lastS2T[ip]=theS2;
00693 lastB2T[ip]=theB2;
00694 lastS3T[ip]=theS3;
00695 lastB3T[ip]=theB3;
00696 lastS4T[ip]=theS4;
00697 lastB4T[ip]=theB4;
00698 #ifdef pdebug
00699 G4cout<<"G4QKaonMinusElasticCrossSection::GetPTables:ip="<<ip<<",lp="<<lp
00700 <<",S1="<<theS1<<",B1="<<theB1<<",S2="<<theS2<<",B2="<<theB2<<",S3="
00701 <<theS3<<",B3="<<theB3<<",S4="<<theS4<<",B4="<<theB4<<G4endl;
00702 #endif
00703 }
00704 return lp;
00705 }
00706 else G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetPTables: PDG="<<PDG
00707 <<", Z="<<tgZ<<", N="<<tgN<<", i="<<ini<<" > fin="<<fin<<", LP="<<LP
00708 <<" > ILP="<<ILP<<" nothing is done!"<<G4endl;
00709 }
00710 else G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetPTables: PDG="<<PDG
00711 <<", Z="<<tgZ<<", N="<<tgN<<", i="<<ini<<">= max="<<nPoints<<", LP="<<LP
00712 <<" > ILP="<<ILP<<", lPMax="<<lPMax<<" nothing is done!"<<G4endl;
00713 }
00714 #ifdef pdebug
00715 else G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetPTa:PDG="<<PDG<<",Z="<<tgZ
00716 <<", N="<<tgN<<", LP="<<LP<<" <= ILP="<<ILP<<" nothing is done!"<<G4endl;
00717 #endif
00718 }
00719 else
00720 {
00721
00722
00723
00724 G4ExceptionDescription ed;
00725 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
00726 << ", while it is defined only for PDG=-321 (K-) " << G4endl;
00727 G4Exception("G4QKaonMinusElasticCrossSection::GetPTables()", "HAD_CHPS_0000",
00728 FatalException, ed);
00729 }
00730 return ILP;
00731 }
00732
00733
00734 G4double G4QKaonMinusElasticCrossSection::GetExchangeT(G4int tgZ, G4int tgN, G4int PDG)
00735 {
00736 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
00737 static const G4double third=1./3.;
00738 static const G4double fifth=1./5.;
00739 static const G4double sevth=1./7.;
00740 #ifdef tdebug
00741 G4cout<<"G4QKaMiElCS::GetExchT:F="<<onlyCS<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
00742 #endif
00743 if(PDG==310 || PDG==130) PDG=-321;
00744 if(PDG!=-321)G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetET:PDG="<<PDG<<G4endl;
00745 if(onlyCS) G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT: onlyCS=1"<<G4endl;
00746 if(lastLP<-4.3) return lastTM*GeVSQ*G4UniformRand();
00747 G4double q2=0.;
00748 if(tgZ==1 && tgN==0)
00749 {
00750 #ifdef tdebug
00751 G4cout<<"G4QKaonMinusElCS::GetExT: TM="<<lastTM<<",S1="<<theS1<<",B1="<<theB1<<",S2="
00752 <<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3<<",GeV2="<<GeVSQ<<G4endl;
00753 #endif
00754 G4double E1=lastTM*theB1;
00755 G4double R1=(1.-std::exp(-E1));
00756 #ifdef tdebug
00757 G4double ts1=-std::log(1.-R1)/theB1;
00758 G4double ds1=std::fabs(ts1-lastTM)/lastTM;
00759 if(ds1>.0001)
00760 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:1p "<<ts1<<"#"<<lastTM
00761 <<",d="<<ds1<<",R1="<<R1<<",E1="<<E1<<G4endl;
00762 #endif
00763 G4double E2=lastTM*theB2;
00764 G4double R2=(1.-std::exp(-E2*E2*E2));
00765 #ifdef tdebug
00766 G4double ts2=std::pow(-std::log(1.-R2),.333333333)/theB2;
00767 G4double ds2=std::fabs(ts2-lastTM)/lastTM;
00768 if(ds2>.0001)
00769 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:2p "<<ts2<<"#"<<lastTM
00770 <<",d="<<ds2<<",R2="<<R2<<",E2="<<E2<<G4endl;
00771 #endif
00772 G4double E3=lastTM*theB3;
00773 G4double R3=(1.-std::exp(-E3));
00774 #ifdef tdebug
00775 G4double ts3=-std::log(1.-R3)/theB3;
00776 G4double ds3=std::fabs(ts3-lastTM)/lastTM;
00777 if(ds3>.0001)
00778 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:3p "<<ts3<<"#"<<lastTM
00779 <<",d="<<ds3<<",R3="<<R1<<",E3="<<E3<<G4endl;
00780 #endif
00781 G4double I1=R1*theS1/theB1;
00782 G4double I2=R2*theS2;
00783 G4double I3=R3*theS3;
00784 G4double I12=I1+I2;
00785 G4double rand=(I12+I3)*G4UniformRand();
00786 if (rand<I1 )
00787 {
00788 G4double ran=R1*G4UniformRand();
00789 if(ran>1.) ran=1.;
00790 q2=-std::log(1.-ran)/theB1;
00791 }
00792 else if(rand<I12)
00793 {
00794 G4double ran=R2*G4UniformRand();
00795 if(ran>1.) ran=1.;
00796 q2=-std::log(1.-ran);
00797 if(q2<0.) q2=0.;
00798 q2=std::pow(q2,third)/theB2;
00799 }
00800 else
00801 {
00802 G4double ran=R3*G4UniformRand();
00803 if(ran>1.) ran=1.;
00804 q2=-std::log(1.-ran)/theB3;
00805 }
00806 }
00807 else
00808 {
00809 G4double a=tgZ+tgN;
00810 #ifdef tdebug
00811 G4cout<<"G4QKaonMinusElasticCrossSection::GetExT:a="<<a<<",t="<<lastTM<<",S1="<<theS1
00812 <<",B1="<<theB1<<",SS="<<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3
00813 <<",B3="<<theB3<<",S4="<<theS4<<",B4="<<theB4<<G4endl;
00814 #endif
00815 G4double E1=lastTM*(theB1+lastTM*theSS);
00816 G4double R1=(1.-std::exp(-E1));
00817 G4double tss=theSS+theSS;
00818 #ifdef tdebug
00819 G4double ts1=-std::log(1.-R1)/theB1;
00820 if(std::fabs(tss)>1.e-7) ts1=(std::sqrt(theB1*(theB1+(tss+tss)*ts1))-theB1)/tss;
00821 G4double ds1=(ts1-lastTM)/lastTM;
00822 if(ds1>.0001)
00823 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:1a "<<ts1<<"#"<<lastTM
00824 <<",d="<<ds1<<",R1="<<R1<<",E1="<<E1<<G4endl;
00825 #endif
00826 G4double tm2=lastTM*lastTM;
00827 G4double E2=lastTM*tm2*theB2;
00828 if(a>6.5)E2*=tm2;
00829 G4double R2=(1.-std::exp(-E2));
00830 #ifdef tdebug
00831 G4double ts2=-std::log(1.-R2)/theB2;
00832 if(a<6.5)ts2=std::pow(ts2,third);
00833 else ts2=std::pow(ts2,fifth);
00834 G4double ds2=std::fabs(ts2-lastTM)/lastTM;
00835 if(ds2>.0001)
00836 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:2a "<<ts2<<"#"<<lastTM
00837 <<",d="<<ds2<<",R2="<<R2<<",E2="<<E2<<G4endl;
00838 #endif
00839 G4double E3=lastTM*theB3;
00840 if(a>6.5)E3*=tm2*tm2*tm2;
00841 G4double R3=(1.-std::exp(-E3));
00842 #ifdef tdebug
00843 G4double ts3=-std::log(1.-R3)/theB3;
00844 if(a>6.5)ts3=std::pow(ts3,sevth);
00845 G4double ds3=std::fabs(ts3-lastTM)/lastTM;
00846 if(ds3>.0001)
00847 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:3a "<<ts3<<"#"<<lastTM
00848 <<",d="<<ds3<<",R3="<<R3<<",E3="<<E3<<G4endl;
00849 #endif
00850 G4double E4=lastTM*theB4;
00851 G4double R4=(1.-std::exp(-E4));
00852 #ifdef tdebug
00853 G4double ts4=-std::log(1.-R4)/theB4;
00854 G4double ds4=std::fabs(ts4-lastTM)/lastTM;
00855 if(ds4>.0001)
00856 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetExT:4a "<<ts4<<"#"<<lastTM
00857 <<",d="<<ds4<<",R4="<<R4<<",E4="<<E4<<G4endl;
00858 #endif
00859 G4double I1=R1*theS1;
00860 G4double I2=R2*theS2;
00861 G4double I3=R3*theS3;
00862 G4double I4=R4*theS4;
00863 G4double I12=I1+I2;
00864 G4double I13=I12+I3;
00865 G4double rand=(I13+I4)*G4UniformRand();
00866 #ifdef tdebug
00867 G4cout<<"G4QKMElCS::GExT:1="<<I1<<",2="<<I2<<",3="<<I3<<",4="<<I4<<",r="<<rand<<G4endl;
00868 #endif
00869 if(rand<I1)
00870 {
00871 G4double ran=R1*G4UniformRand();
00872 if(ran>1.) ran=1.;
00873 q2=-std::log(1.-ran)/theB1;
00874 if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
00875 #ifdef tdebug
00876 G4cout<<"G4QKMElCS::GExT:Q2="<<q2<<",ss="<<tss/2<<",b1="<<theB1<<",t1="<<ts1<<G4endl;
00877 #endif
00878 }
00879 else if(rand<I12)
00880 {
00881 G4double ran=R2*G4UniformRand();
00882 if(ran>1.) ran=1.;
00883 q2=-std::log(1.-ran)/theB2;
00884 if(q2<0.) q2=0.;
00885 if(a<6.5) q2=std::pow(q2,third);
00886 else q2=std::pow(q2,fifth);
00887 #ifdef tdebug
00888 G4cout<<"G4QKMElCS::GetExT: Q2="<<q2<<",r2="<<R2<<",b2="<<theB2<<",t2="<<ts2<<G4endl;
00889 #endif
00890 }
00891 else if(rand<I13)
00892 {
00893 G4double ran=R3*G4UniformRand();
00894 if(ran>1.) ran=1.;
00895 q2=-std::log(1.-ran)/theB3;
00896 if(q2<0.) q2=0.;
00897 if(a>6.5) q2=std::pow(q2,sevth);
00898 #ifdef tdebug
00899 G4cout<<"G4QKMElCS::GetExT: Q2="<<q2<<",r3="<<R2<<",b3="<<theB2<<",t3="<<ts2<<G4endl;
00900 #endif
00901 }
00902 else
00903 {
00904 G4double ran=R4*G4UniformRand();
00905 if(ran>1.) ran=1.;
00906 q2=-std::log(1.-ran)/theB4;
00907 if(a<6.5) q2=lastTM-q2;
00908 #ifdef tdebug
00909 G4cout<<"G4QKMElCS::GET:Q2="<<q2<<",m="<<lastTM<<",b4="<<theB3<<",t4="<<ts3<<G4endl;
00910 #endif
00911 }
00912 }
00913 if(q2<0.) q2=0.;
00914 if(!(q2>=-1.||q2<=1.)) G4cout<<"*NAN*G4QKaonMinusElasticCS::GetExchT: -t="<<q2<<G4endl;
00915 if(q2>lastTM)
00916 {
00917 #ifdef tdebug
00918 G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GET:-t="<<q2<<">"<<lastTM<<G4endl;
00919 #endif
00920 q2=lastTM;
00921 }
00922 return q2*GeVSQ;
00923 }
00924
00925
00926 G4double G4QKaonMinusElasticCrossSection::GetSlope(G4int tgZ, G4int tgN, G4int PDG)
00927 {
00928 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
00929 #ifdef tdebug
00930 G4cout<<"G4QKaonMinusECS::GetS:"<<onlyCS<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
00931 #endif
00932 if(onlyCS)G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetSl:onlCS=true"<<G4endl;
00933 if(lastLP<-4.3) return 0.;
00934 if(PDG != -321)
00935 {
00936
00937
00938
00939 G4ExceptionDescription ed;
00940 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
00941 << ", while it is defined only for PDG=-321 (K-)" << G4endl;
00942 }
00943 if(theB1<0.) theB1=0.;
00944 if(!(theB1>=-1.||theB1<=1.))G4cout<<"*NAN*G4QKaonMinusElCS::GetSlope:B1="<<theB1<<G4endl;
00945 return theB1/GeVSQ;
00946 }
00947
00948
00949 G4double G4QKaonMinusElasticCrossSection::GetHMaxT()
00950 {
00951 static const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
00952 return lastTM*HGeVSQ;
00953 }
00954
00955
00956 G4double G4QKaonMinusElasticCrossSection::GetTabValues(G4double lp, G4int PDG, G4int tgZ,
00957 G4int tgN)
00958 {
00959 if(PDG!=-321)G4cout<<"*Warning*G4QKaonMinusElasticCrossSection::GetTV:PDG="<<PDG<<G4endl;
00960 if(tgZ<0 || tgZ>92)
00961 {
00962 G4cout<<"*Warning*G4QKaonMinusElasticCS::GetTabV:(1-92)NoIsotopes for Z="<<tgZ<<G4endl;
00963 return 0.;
00964 }
00965 G4int iZ=tgZ-1;
00966 if(iZ<0)
00967 {
00968 iZ=0;
00969 tgZ=1;
00970 tgN=0;
00971 }
00972
00973
00974 #ifdef isodebug
00975
00976 #endif
00977
00978
00979 #ifdef pdebug
00980 G4cout<<"G4QKaonMinusECS::GetTV: l="<<lp<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
00981 #endif
00982 G4double p=std::exp(lp);
00983 G4double sp=std::sqrt(p);
00984 G4double psp=p*sp;
00985 G4double p2=p*p;
00986 G4double p3=p2*p;
00987 G4double p4=p3*p;
00988 if ( tgZ == 1 && tgN == 0 )
00989 {
00990 G4double dl2=lp-lastPAR[12];
00991 theSS=lastPAR[35];
00992 theS1=(lastPAR[13]+lastPAR[14]*dl2*dl2)/(1.+lastPAR[15]/p4/p)+
00993 (lastPAR[16]/p2+lastPAR[17]*p)/(p4+lastPAR[18]*sp);
00994 theB1=lastPAR[19]*std::pow(p,lastPAR[20])/(1.+lastPAR[21]/p3);
00995 theS2=lastPAR[22]+lastPAR[23]/(p4+lastPAR[24]*p);
00996 theB2=lastPAR[25]+lastPAR[26]/(p4+lastPAR[27]/sp);
00997 theS3=lastPAR[28]+lastPAR[29]/(p4*p4+lastPAR[30]*p2+lastPAR[31]);
00998 theB3=lastPAR[32]+lastPAR[33]/(p4+lastPAR[34]);
00999 theS4=0.;
01000 theB4=0.;
01001 #ifdef tdebug
01002 G4cout<<"G4QKaonMinusElasticCrossSection::GetTabV:TM="<<lastTM<<",S1="<<theS1<<",B1="
01003 <<theB1<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS1<<",B3="<<theB1<<G4endl;
01004 #endif
01005
01006 G4double dp=lp-lastPAR[2];
01007 return lastPAR[0]/psp+(lastPAR[1]*dp*dp+lastPAR[3])/(1.-lastPAR[4]/sp+lastPAR[5]/p4)+
01008 lastPAR[6]/(sqr(p-lastPAR[7])+lastPAR[8])+lastPAR[9]/(sqr(p-lastPAR[10])+lastPAR[11]);
01009 }
01010 else
01011 {
01012 G4double p5=p4*p;
01013 G4double p6=p5*p;
01014 G4double p8=p6*p2;
01015 G4double p10=p8*p2;
01016 G4double p12=p10*p2;
01017 G4double p16=p8*p8;
01018
01019 G4double dl=lp-5.;
01020 G4double a=tgZ+tgN;
01021 G4double pah=std::pow(p,a/2);
01022 G4double pa=pah*pah;
01023 G4double pa2=pa*pa;
01024 if(a<6.5)
01025 {
01026 theS1=lastPAR[9]/(1.+lastPAR[10]*p4*pa)+lastPAR[11]/(p4+lastPAR[12]*p4/pa2)+
01027 (lastPAR[13]*dl*dl+lastPAR[14])/(1.+lastPAR[15]/p2);
01028 theB1=(lastPAR[16]+lastPAR[17]*p2)/(p4+lastPAR[18]/pah)+lastPAR[19];
01029 theSS=lastPAR[20]/(1.+lastPAR[21]/p2)+lastPAR[22]/(p6/pa+lastPAR[23]/p16);
01030 theS2=lastPAR[24]/(pa/p2+lastPAR[25]/p4)+lastPAR[26];
01031 theB2=lastPAR[27]*std::pow(p,lastPAR[28])+lastPAR[29]/(p8+lastPAR[30]/p16);
01032 theS3=lastPAR[31]/(pa*p+lastPAR[32]/pa)+lastPAR[33];
01033 theB3=lastPAR[34]/(p3+lastPAR[35]/p6)+lastPAR[36]/(1.+lastPAR[37]/p2);
01034 theS4=p2*(pah*lastPAR[38]*std::exp(-pah*lastPAR[39])+
01035 lastPAR[40]/(1.+lastPAR[41]*std::pow(p,lastPAR[42])));
01036 theB4=lastPAR[43]*pa/p2/(1.+pa*lastPAR[44]);
01037 #ifdef tdebug
01038 G4cout<<"G4QKaonMinusElasticCS::GetTabV: lA, p="<<p<<",S1="<<theS1<<",B1="<<theB1
01039 <<",SS="<<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3
01040 <<",S4="<<theS4<<",B4="<<theB4<<G4endl;
01041 #endif
01042 }
01043 else
01044 {
01045 theS1=lastPAR[9]/(1.+lastPAR[10]/p4)+lastPAR[11]/(p4+lastPAR[12]/p2)+
01046 lastPAR[13]/(p5+lastPAR[14]/p16);
01047 theB1=(lastPAR[15]/p8+lastPAR[19])/(p+lastPAR[16]/std::pow(p,lastPAR[20]))+
01048 lastPAR[17]/(1.+lastPAR[18]/p4);
01049 theSS=lastPAR[21]/(p4/std::pow(p,lastPAR[23])+lastPAR[22]/p4);
01050 theS2=lastPAR[24]/p4/(std::pow(p,lastPAR[25])+lastPAR[26]/p12)+lastPAR[27];
01051 theB2=lastPAR[28]/std::pow(p,lastPAR[29])+lastPAR[30]/std::pow(p,lastPAR[31]);
01052 theS3=lastPAR[32]/std::pow(p,lastPAR[35])/(1.+lastPAR[36]/p12)+
01053 lastPAR[33]/(1.+lastPAR[34]/p6);
01054 theB3=lastPAR[37]/p8+lastPAR[38]/p2+lastPAR[39]/(1.+lastPAR[40]/p8);
01055 theS4=(lastPAR[41]/p4+lastPAR[46]/p)/(1.+lastPAR[42]/p10)+
01056 (lastPAR[43]+lastPAR[44]*dl*dl)/(1.+lastPAR[45]/p12);
01057 theB4=lastPAR[47]/(1.+lastPAR[48]/p)+lastPAR[49]*p4/(1.+lastPAR[50]*p5);
01058 #ifdef tdebug
01059 G4cout<<"G4QKaonMinusElasticCS::GetTabVal:hA,p="<<p<<",S1="<<theS1<<",B1="<<theB1
01060 <<",SS="<<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3
01061 <<",S4="<<theS4<<",B4="<<theB4<<G4endl;
01062 #endif
01063 }
01064
01065 #ifdef tdebug
01066 G4cout<<"G4QKaonMinusElCS::GTV: PDG="<<PDG<<",P="<<p<<",N="<<tgN<<",Z="<<tgZ<<G4endl;
01067 #endif
01068 G4double dlp=lp-lastPAR[4];
01069
01070 return (lastPAR[0]*dlp*dlp+lastPAR[1]+lastPAR[2]/p3)/(1.+lastPAR[3]/p2/sp);
01071 }
01072 return 0.;
01073 }
01074
01075
01076 G4double G4QKaonMinusElasticCrossSection::GetQ2max(G4int PDG, G4int tgZ, G4int tgN,
01077 G4double pP)
01078 {
01079
01080
01081 static const G4double mK= G4QPDGCode(321).GetMass()*.001;
01082
01083
01084
01085
01086
01087 static const G4double mK2= mK*mK;
01088
01089
01090
01091
01092 G4double pP2=pP*pP;
01093 if(tgZ || tgN>-1)
01094 {
01095 G4double mt=G4QPDGCode(90000000+tgZ*1000+tgN).GetMass()*.001;
01096 G4double dmt=mt+mt;
01097 G4double s_value=dmt*std::sqrt(pP2+mK2)+mK2+mt*mt;
01098 return dmt*dmt*pP2/s_value;
01099 }
01100 else
01101 {
01102
01103
01104
01105 G4ExceptionDescription ed;
01106 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
01107 << ", while it is defined only for p projectiles & Z_target>0" << G4endl;
01108 G4Exception("G4QKaonMinusElasticCrossSection::GetQ2max()", "HAD_CHPS_0000",
01109 FatalException, ed);
01110 return 0;
01111 }
01112 }