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 "G4QKaonPlusElasticCrossSection.hh"
00046 #include "G4SystemOfUnits.hh"
00047
00048
00049 const G4int G4QKaonPlusElasticCrossSection::nPoints=128;
00050 const G4int G4QKaonPlusElasticCrossSection::nLast=nPoints-1;
00051 G4double G4QKaonPlusElasticCrossSection::lPMin=-8.;
00052 G4double G4QKaonPlusElasticCrossSection::lPMax= 8.;
00053 G4double G4QKaonPlusElasticCrossSection::dlnP=(lPMax-lPMin)/nLast;
00054 G4bool G4QKaonPlusElasticCrossSection::onlyCS=true;
00055 G4double G4QKaonPlusElasticCrossSection::lastSIG=0.;
00056 G4double G4QKaonPlusElasticCrossSection::lastLP=-10.;
00057 G4double G4QKaonPlusElasticCrossSection::lastTM=0.;
00058 G4double G4QKaonPlusElasticCrossSection::theSS=0.;
00059 G4double G4QKaonPlusElasticCrossSection::theS1=0.;
00060 G4double G4QKaonPlusElasticCrossSection::theB1=0.;
00061 G4double G4QKaonPlusElasticCrossSection::theS2=0.;
00062 G4double G4QKaonPlusElasticCrossSection::theB2=0.;
00063 G4double G4QKaonPlusElasticCrossSection::theS3=0.;
00064 G4double G4QKaonPlusElasticCrossSection::theB3=0.;
00065 G4double G4QKaonPlusElasticCrossSection::theS4=0.;
00066 G4double G4QKaonPlusElasticCrossSection::theB4=0.;
00067 G4int G4QKaonPlusElasticCrossSection::lastTZ=0;
00068 G4int G4QKaonPlusElasticCrossSection::lastTN=0;
00069 G4double G4QKaonPlusElasticCrossSection::lastPIN=0.;
00070 G4double* G4QKaonPlusElasticCrossSection::lastCST=0;
00071 G4double* G4QKaonPlusElasticCrossSection::lastPAR=0;
00072 G4double* G4QKaonPlusElasticCrossSection::lastSST=0;
00073 G4double* G4QKaonPlusElasticCrossSection::lastS1T=0;
00074 G4double* G4QKaonPlusElasticCrossSection::lastB1T=0;
00075 G4double* G4QKaonPlusElasticCrossSection::lastS2T=0;
00076 G4double* G4QKaonPlusElasticCrossSection::lastB2T=0;
00077 G4double* G4QKaonPlusElasticCrossSection::lastS3T=0;
00078 G4double* G4QKaonPlusElasticCrossSection::lastB3T=0;
00079 G4double* G4QKaonPlusElasticCrossSection::lastS4T=0;
00080 G4double* G4QKaonPlusElasticCrossSection::lastB4T=0;
00081 G4int G4QKaonPlusElasticCrossSection::lastN=0;
00082 G4int G4QKaonPlusElasticCrossSection::lastZ=0;
00083 G4double G4QKaonPlusElasticCrossSection::lastP=0.;
00084 G4double G4QKaonPlusElasticCrossSection::lastTH=0.;
00085 G4double G4QKaonPlusElasticCrossSection::lastCS=0.;
00086 G4int G4QKaonPlusElasticCrossSection::lastI=0;
00087
00088 std::vector<G4double*> G4QKaonPlusElasticCrossSection::PAR;
00089 std::vector<G4double*> G4QKaonPlusElasticCrossSection::CST;
00090 std::vector<G4double*> G4QKaonPlusElasticCrossSection::SST;
00091 std::vector<G4double*> G4QKaonPlusElasticCrossSection::S1T;
00092 std::vector<G4double*> G4QKaonPlusElasticCrossSection::B1T;
00093 std::vector<G4double*> G4QKaonPlusElasticCrossSection::S2T;
00094 std::vector<G4double*> G4QKaonPlusElasticCrossSection::B2T;
00095 std::vector<G4double*> G4QKaonPlusElasticCrossSection::S3T;
00096 std::vector<G4double*> G4QKaonPlusElasticCrossSection::B3T;
00097 std::vector<G4double*> G4QKaonPlusElasticCrossSection::S4T;
00098 std::vector<G4double*> G4QKaonPlusElasticCrossSection::B4T;
00099
00100 G4QKaonPlusElasticCrossSection::G4QKaonPlusElasticCrossSection()
00101 {
00102 }
00103
00104 G4QKaonPlusElasticCrossSection::~G4QKaonPlusElasticCrossSection()
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* G4QKaonPlusElasticCrossSection::GetPointer()
00144 {
00145 static G4QKaonPlusElasticCrossSection theCrossSection;
00146 return &theCrossSection;
00147 }
00148
00149
00150
00151 G4double G4QKaonPlusElasticCrossSection::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<<"G4QKaPlusElCS::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*G4QKaonPlusElaCS::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<<"G4QKPElCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",i="<<i<<G4endl;
00186
00187 #endif
00188 if(pEn<=lastTH)
00189 {
00190 #ifdef pdebug
00191 G4cout<<"G4QKPElCS::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<<"G4QKaonPlusElastCS::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<<"G4QKPElCS::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<<"G4QKaonPlusElCS::GetCrosSec:****>New(inDB) Calculated CS="<<lastCS<<G4endl;
00215
00216 #endif
00217 if(lastCS<=0. && pEn>lastTH)
00218 {
00219 #ifdef pdebug
00220 G4cout<<"G4QKaonPlusElCS::GetCS: T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00221 #endif
00222 lastTH=pEn;
00223 }
00224 break;
00225 }
00226 #ifdef pdebug
00227 G4cout<<"---G4QKaonPlusElasticCrossSection::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<<"G4QKaonPlusElCS::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<<"G4QKaonPlusElasticCrossSection::GetCS:NewThr="<<lastTH<<",T="<<pEn<<G4endl;
00244 #endif
00245 if(pEn>lastTH)
00246 {
00247 #ifdef pdebug
00248 G4cout<<"G4QKaonPlusElCS::GetCS:1st T="<<pEn<<"(CS=0)>Threshold="<<lastTH<<G4endl;
00249 #endif
00250 lastTH=pEn;
00251 }
00252 }
00253 #ifdef pdebug
00254 G4cout<<"G4QKaonPlusElCS::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<<"G4QKPElCS::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<<"G4QKaonPlusElasticCrossSection::GetCS: Update lastI="<<lastI<<G4endl;
00272 #endif
00273 colP[lastI]=pMom;
00274 colCS[lastI]=lastCS;
00275 }
00276 #ifdef pdebug
00277 G4cout<<"G4QKPElCS::GetCSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00278
00279 G4cout<<"G4QKaonPlusElasticCrossSection::GetCrosSec:**End**, onlyCS="<<onlyCS<<G4endl;
00280 #endif
00281 return lastCS*millibarn;
00282 }
00283
00284
00285
00286 G4double G4QKaonPlusElasticCrossSection::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<<"G4QKaonPlusElasticCS::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<<"G4QKaonPlusElCS::CalCS:DB updated for I="<<I<<",*,PIN4="<<PIN[4]<<G4endl;
00316 #endif
00317 }
00318 #ifdef pdebug
00319 G4cout<<"G4QKaonPlusElasticCS::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<<"G4QKPElCS::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<<"G4QKaonPlusElasticCS::CalcCS:*ini*,lstLP="<<lastLP<<",min="<<lPMin<<G4endl;
00346 #endif
00347 lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN);
00348 #ifdef pdebug
00349 G4cout<<"G4QKaPlElCS::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<<"G4QKPElCS::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<<"G4QKaonPlusElCS::CalcCS:*updated(O)*,LP="<<lastLP<<"<IN="<<lastPIN<<G4endl;
00373 #endif
00374 }
00375 #ifdef pdebug
00376 G4cout<<"G4QKaPlElCS::CalcCS:lastLP="<<lastLP<<",lPM="<<lPMin<<",lPIN="<<lastPIN<<G4endl;
00377 #endif
00378 if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom);
00379 #ifdef pdebug
00380 G4cout<<"G4QKaPlElCrosSec::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<<"G4QKPElCS::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<<"G4QKaonPlusElasticCS::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<<"G4QKaonPlusElasticCrossSection::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<<"G4QKaonPlusElasticCS::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<<"G4QKaonPlusElasticCrossSection::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<<"G4QKaonPlusElasticCrossSection::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<<"G4QKaonPlusElasticCS::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<<"G4QKaonPlusElasticCrossSection::CalculateCS: END, onlyCS="<<onlyCS<<G4endl;
00462 #endif
00463 return lastSIG;
00464 }
00465
00466
00467 G4double G4QKaonPlusElasticCrossSection::GetPTables(G4double LP,G4double ILP, G4int PDG,
00468 G4int tgZ, G4int tgN)
00469 {
00470
00471 static const G4double pwd=2727;
00472 const G4int n_kppel=35;
00473
00474 G4double kpp_el[n_kppel]={.7,.38,.0676,.0557,3.5,2.23,.7,.1,2.,1.,.372,5.,74.,3.,3.4,
00475 .2,.17,.001,8.,.055,3.64,5.e-5,4000.,1500.,.46,1.2e6,3.5e6,
00476 5.e-5,1.e10,8.5e8,1.e10,1.1,3.4e6,6.8e6,0.};
00477
00478
00479
00480 if(PDG == 321)
00481 {
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493 if(lastPAR[nLast]!=pwd)
00494 {
00495 if ( tgZ == 1 && tgN == 0 )
00496 {
00497 for (G4int ip=0; ip<n_kppel; ip++) lastPAR[ip]=kpp_el[ip];
00498 }
00499 else
00500 {
00501 G4double a=tgZ+tgN;
00502 G4double sa=std::sqrt(a);
00503 G4double ssa=std::sqrt(sa);
00504 G4double asa=a*sa;
00505 G4double a2=a*a;
00506 G4double a3=a2*a;
00507 G4double a4=a3*a;
00508 G4double a5=a4*a;
00509 G4double a6=a4*a2;
00510 G4double a7=a6*a;
00511 G4double a8=a7*a;
00512 G4double a9=a8*a;
00513 G4double a10=a5*a5;
00514 G4double a12=a6*a6;
00515 G4double a14=a7*a7;
00516 G4double a16=a8*a8;
00517 G4double a17=a16*a;
00518
00519 G4double a32=a16*a16;
00520
00521 lastPAR[0]=.06*asa/(1.+a*(.01+.1/ssa));
00522 lastPAR[1]=.75*asa/(1.+.009*a);
00523 lastPAR[2]=.9*asa*ssa/(1.+.03*a);
00524 lastPAR[3]=3.;
00525 lastPAR[4]=4.2;
00526 lastPAR[5]=0.;
00527 lastPAR[6]=0.;
00528 lastPAR[7]=0.;
00529 lastPAR[8]=0.;
00530
00531 if(a<6.5)
00532 {
00533 G4double a28=a16*a12;
00534
00535 lastPAR[ 9]=4000*a;
00536 lastPAR[10]=1.2e7*a8+380*a17;
00537 lastPAR[11]=.7/(1.+4.e-12*a16);
00538 lastPAR[12]=2.5/a8/(a4+1.e-16*a32);
00539 lastPAR[13]=.28*a;
00540 lastPAR[14]=1.2*a2+2.3;
00541 lastPAR[15]=3.8/a;
00542
00543 lastPAR[16]=.01/(1.+.0024*a5);
00544 lastPAR[17]=.2*a;
00545 lastPAR[18]=9.e-7/(1.+.035*a5);
00546 lastPAR[19]=(42.+2.7e-11*a16)/(1.+.14*a);
00547
00548 lastPAR[20]=2.25*a3;
00549 lastPAR[21]=18.;
00550 lastPAR[22]=2.4e-3*a8/(1.+2.6e-4*a7);
00551 lastPAR[23]=3.5e-36*a32*a8/(1.+5.e-15*a32/a);
00552
00553 lastPAR[24]=1.e5/(a8+2.5e12/a16);
00554 lastPAR[25]=8.e7/(a12+1.e-27*a28*a28);
00555 lastPAR[26]=.0006*a3;
00556
00557 lastPAR[27]=10.+4.e-8*a12*a;
00558 lastPAR[28]=.114;
00559 lastPAR[29]=.003;
00560 lastPAR[30]=2.e-23;
00561
00562 lastPAR[31]=1./(1.+.0001*a8);
00563 lastPAR[32]=1.5e-4/(1.+5.e-6*a12);
00564 lastPAR[33]=.03;
00565
00566 lastPAR[34]=a/2;
00567 lastPAR[35]=2.e-7*a4;
00568 lastPAR[36]=4.;
00569 lastPAR[37]=64./a3;
00570
00571 lastPAR[38]=1.e8*std::exp(.32*asa);
00572 lastPAR[39]=20.*std::exp(.45*asa);
00573 lastPAR[40]=7.e3+2.4e6/a5;
00574 lastPAR[41]=2.5e5*std::exp(.085*a3);
00575 lastPAR[42]=2.5*a;
00576
00577 lastPAR[43]=920.+.03*a8*a3;
00578 lastPAR[44]=93.+.0023*a12;
00579 #ifdef pdebug
00580 G4cout<<"G4QKPElCS::CalcCS:la "<<lastPAR[38]<<", "<<lastPAR[39]<<", "<<lastPAR[40]
00581 <<", "<<lastPAR[42]<<", "<<lastPAR[43]<<", "<<lastPAR[44]<<G4endl;
00582 #endif
00583 }
00584 else
00585 {
00586 G4double p1a10=2.2e-28*a10;
00587 G4double r4a16=6.e14/a16;
00588 G4double s4a16=r4a16*r4a16;
00589
00590
00591
00592 lastPAR[ 9]=4.5*std::pow(a,1.15);
00593 lastPAR[10]=.06*std::pow(a,.6);
00594 lastPAR[11]=.6*a/(1.+2.e15/a16);
00595 lastPAR[12]=.17/(a+9.e5/a3+1.5e33/a32);
00596 lastPAR[13]=(.001+7.e-11*a5)/(1.+4.4e-11*a5);
00597 lastPAR[14]=(p1a10*p1a10+2.e-29)/(1.+2.e-22*a12);
00598
00599 lastPAR[15]=400./a12+2.e-22*a9;
00600 lastPAR[16]=1.e-32*a12/(1.+5.e22/a14);
00601 lastPAR[17]=1000./a2+9.5*sa*ssa;
00602 lastPAR[18]=4.e-6*a*asa+1.e11/a16;
00603 lastPAR[19]=(120./a+.002*a2)/(1.+2.e14/a16);
00604 lastPAR[20]=9.+100./a;
00605
00606 lastPAR[21]=.002*a3+3.e7/a6;
00607 lastPAR[22]=7.e-15*a4*asa;
00608 lastPAR[23]=9000./a4;
00609
00610 lastPAR[24]=.0011*asa/(1.+3.e34/a32/a4);
00611 lastPAR[25]=1.e-5*a2+2.e14/a16;
00612 lastPAR[26]=1.2e-11*a2/(1.+1.5e19/a12);
00613 lastPAR[27]=.016*asa/(1.+5.e16/a16);
00614
00615 lastPAR[28]=.002*a4/(1.+7.e7/std::pow(a-6.83,14));
00616 lastPAR[29]=2.e6/a6+7.2/std::pow(a,.11);
00617 lastPAR[30]=11.*a3/(1.+7.e23/a16/a8);
00618 lastPAR[31]=100./asa;
00619
00620 lastPAR[32]=(.1+4.4e-5*a2)/(1.+5.e5/a4);
00621 lastPAR[33]=3.5e-4*a2/(1.+1.e8/a8);
00622 lastPAR[34]=1.3+3.e5/a4;
00623 lastPAR[35]=500./(a2+50.)+3;
00624 lastPAR[36]=1.e-9/a+s4a16*s4a16;
00625
00626 lastPAR[37]=.4*asa+3.e-9*a6;
00627 lastPAR[38]=.0005*a5;
00628 lastPAR[39]=.002*a5;
00629 lastPAR[40]=10.;
00630
00631 lastPAR[41]=.05+.005*a;
00632 lastPAR[42]=7.e-8/sa;
00633 lastPAR[43]=.8*sa;
00634 lastPAR[44]=.02*sa;
00635 lastPAR[45]=1.e8/a3;
00636 lastPAR[46]=3.e32/(a32+1.e32);
00637
00638 lastPAR[47]=24.;
00639 lastPAR[48]=20./sa;
00640 lastPAR[49]=7.e3*a/(sa+1.);
00641 lastPAR[50]=900.*sa/(1.+500./a3);
00642 #ifdef pdebug
00643 G4cout<<"G4QKPElCS::CalcCS:ha "<<lastPAR[41]<<", "<<lastPAR[42]<<", "<<lastPAR[43]
00644 <<", "<<lastPAR[44]<<", "<<lastPAR[45]<<", "<<lastPAR[46]<<G4endl;
00645 #endif
00646 }
00647
00648 lastPAR[51]=1.e15+2.e27/a4/(1.+2.e-18*a16);
00649 }
00650 lastPAR[nLast]=pwd;
00651
00652 G4double lp=lPMin;
00653 G4bool memCS=onlyCS;
00654 onlyCS=false;
00655 lastCST[0]=GetTabValues(lp, PDG, tgZ, tgN);
00656 onlyCS=memCS;
00657 lastSST[0]=theSS;
00658 lastS1T[0]=theS1;
00659 lastB1T[0]=theB1;
00660 lastS2T[0]=theS2;
00661 lastB2T[0]=theB2;
00662 lastS3T[0]=theS3;
00663 lastB3T[0]=theB3;
00664 lastS4T[0]=theS4;
00665 lastB4T[0]=theB4;
00666 #ifdef pdebug
00667 G4cout<<"G4QKaonPlusElasticCrossSection::GetPTables:ip=0(init), lp="<<lp<<",S1="
00668 <<theS1<<",B1="<<theB1<<",S2="<<theS2<<",B2="<<theB3<<",S3="<<theS3
00669 <<",B3="<<theB3<<",S4="<<theS4<<",B4="<<theB4<<G4endl;
00670 #endif
00671 }
00672 if(LP>ILP)
00673 {
00674 G4int ini = static_cast<int>((ILP-lPMin+.000001)/dlnP)+1;
00675 if(ini<0) ini=0;
00676 if(ini<nPoints)
00677 {
00678 G4int fin = static_cast<int>((LP-lPMin)/dlnP)+1;
00679 if(fin>=nPoints) fin=nLast;
00680 if(fin>=ini)
00681 {
00682 G4double lp=0.;
00683 for(G4int ip=ini; ip<=fin; ip++)
00684 {
00685 lp=lPMin+ip*dlnP;
00686 G4bool memCS=onlyCS;
00687 onlyCS=false;
00688 lastCST[ip]=GetTabValues(lp, PDG, tgZ, tgN);
00689 onlyCS=memCS;
00690 lastSST[ip]=theSS;
00691 lastS1T[ip]=theS1;
00692 lastB1T[ip]=theB1;
00693 lastS2T[ip]=theS2;
00694 lastB2T[ip]=theB2;
00695 lastS3T[ip]=theS3;
00696 lastB3T[ip]=theB3;
00697 lastS4T[ip]=theS4;
00698 lastB4T[ip]=theB4;
00699 #ifdef pdebug
00700 G4cout<<"G4QKaonPlusElasticCrossSection::GetPTables:ip="<<ip<<",lp="<<lp
00701 <<",S1="<<theS1<<",B1="<<theB1<<",S2="<<theS2<<",B2="<<theB2<<",S3="
00702 <<theS3<<",B3="<<theB3<<",S4="<<theS4<<",B4="<<theB4<<G4endl;
00703 #endif
00704 }
00705 return lp;
00706 }
00707 else G4cout<<"*Warning*G4QKaonPlusElasticCrossSection::GetPTables: PDG="<<PDG
00708 <<", Z="<<tgZ<<", N="<<tgN<<", i="<<ini<<" > fin="<<fin<<", LP="<<LP
00709 <<" > ILP="<<ILP<<" nothing is done!"<<G4endl;
00710 }
00711 else G4cout<<"*Warning*G4QKaonPlusElasticCrossSection::GetPTables: PDG="<<PDG
00712 <<", Z="<<tgZ<<", N="<<tgN<<", i="<<ini<<">= max="<<nPoints<<", LP="<<LP
00713 <<" > ILP="<<ILP<<", lPMax="<<lPMax<<" nothing is done!"<<G4endl;
00714 }
00715 #ifdef pdebug
00716 else G4cout<<"*Warning*G4QKaonPlusElasticCrossSection::GetPTa:PDG="<<PDG<<",Z="<<tgZ
00717 <<", N="<<tgN<<", LP="<<LP<<" <= ILP="<<ILP<<" nothing is done!"<<G4endl;
00718 #endif
00719 }
00720 else
00721 {
00722
00723
00724
00725 G4ExceptionDescription ed;
00726 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
00727 << ", while it is defined only for PDG=321 (K+) " << G4endl;
00728 G4Exception("G4QKaonPlusElasticCrossSection::GetPTables()", "HAD_CHPS_0000",
00729 FatalException, ed);
00730 }
00731 return ILP;
00732 }
00733
00734
00735 G4double G4QKaonPlusElasticCrossSection::GetExchangeT(G4int tgZ, G4int tgN, G4int PDG)
00736 {
00737 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
00738 static const G4double third=1./3.;
00739 static const G4double fifth=1./5.;
00740 static const G4double sevth=1./7.;
00741 #ifdef tdebug
00742 G4cout<<"G4QKaPlElCS::GetExchT:F="<<onlyCS<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
00743 #endif
00744 if(PDG!=321) G4cout<<"*Warning*G4QKaonPlusElasticCrossSection::GetExT:PDG="<<PDG<<G4endl;
00745 if(onlyCS) G4cout<<"*Warning*G4QKaonPlusElasticCrossSection::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<<"G4QKaonPlusElCS::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*G4QKaonPlusElasticCrossSection::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*G4QKaonPlusElasticCrossSection::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*G4QKaonPlusElasticCrossSection::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<<"G4QKaonPlusElasticCrossSection::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*G4QKaonPlusElasticCrossSection::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*G4QKaonPlusElasticCrossSection::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*G4QKaonPlusElasticCrossSection::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*G4QKaonPlusElasticCrossSection::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<<"G4QKPElCS::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<<"G4QKPElCS::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<<"G4QKPElCS::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<<"G4QKPElCS::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<<"G4QKPElCS::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*G4QKaonPlusElasticCS::GetExchT: -t="<<q2<<G4endl;
00915 if(q2>lastTM)
00916 {
00917 #ifdef tdebug
00918 G4cout<<"*Warning*G4QKaonPlusElasticCrossSection::GET:-t="<<q2<<">"<<lastTM<<G4endl;
00919 #endif
00920 q2=lastTM;
00921 }
00922 return q2*GeVSQ;
00923 }
00924
00925
00926 G4double G4QKaonPlusElasticCrossSection::GetSlope(G4int tgZ, G4int tgN, G4int PDG)
00927 {
00928 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
00929 #ifdef tdebug
00930 G4cout<<"G4QKaonPlusECS::GetS:"<<onlyCS<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
00931 #endif
00932 if(onlyCS)G4cout<<"*Warning*G4QKaonPlusElasticCrossSection::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 G4Exception("G4QKaonPlusElasticCrossSection::GetSlope()", "HAD_CHPS_0000",
00943 FatalException, ed);
00944 }
00945 if(theB1<0.) theB1=0.;
00946 if(!(theB1>=-1.||theB1<=1.))G4cout<<"*NAN*G4QKaonPlusElCS::GetSlope:B1="<<theB1<<G4endl;
00947 return theB1/GeVSQ;
00948 }
00949
00950
00951 G4double G4QKaonPlusElasticCrossSection::GetHMaxT()
00952 {
00953 static const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
00954 return lastTM*HGeVSQ;
00955 }
00956
00957
00958 G4double G4QKaonPlusElasticCrossSection::GetTabValues(G4double lp, G4int PDG, G4int tgZ,
00959 G4int tgN)
00960 {
00961 if(PDG!=321)G4cout<<"*Warning*G4QKaonPlusElasticCrossSection::GetTaV:PDG="<<PDG<<G4endl;
00962 if(tgZ<0 || tgZ>92)
00963 {
00964 G4cout<<"*Warning*G4QKaonPlusElasticCS::GetTabV:(1-92)NoIsotopes for Z="<<tgZ<<G4endl;
00965 return 0.;
00966 }
00967 G4int iZ=tgZ-1;
00968 if(iZ<0)
00969 {
00970 iZ=0;
00971 tgZ=1;
00972 tgN=0;
00973 }
00974
00975
00976 #ifdef isodebug
00977
00978 #endif
00979
00980
00981 #ifdef pdebug
00982 G4cout<<"G4QKaonPlusECS::GetTV: l="<<lp<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
00983 #endif
00984 G4double p=std::exp(lp);
00985 G4double sp=std::sqrt(p);
00986 G4double p2=p*p;
00987 G4double p3=p2*p;
00988 G4double p4=p3*p;
00989 if ( tgZ == 1 && tgN == 0 )
00990 {
00991 G4double dl2=lp-lastPAR[11];
00992 theSS=lastPAR[34];
00993 theS1=(lastPAR[12]+lastPAR[13]*dl2*dl2)/(1.+lastPAR[14]/p4/p)+
00994 (lastPAR[15]/p2+lastPAR[16]*p)/(p4+lastPAR[17]*sp);
00995 theB1=lastPAR[18]*std::pow(p,lastPAR[19])/(1.+lastPAR[20]/p3);
00996 theS2=lastPAR[21]+lastPAR[22]/(p4+lastPAR[23]*p);
00997 theB2=lastPAR[24]+lastPAR[25]/(p4+lastPAR[26]/sp);
00998 theS3=lastPAR[27]+lastPAR[28]/(p4*p4+lastPAR[29]*p2+lastPAR[30]);
00999 theB3=lastPAR[31]+lastPAR[32]/(p4+lastPAR[33]);
01000 theS4=0.;
01001 theB4=0.;
01002 #ifdef tdebug
01003 G4cout<<"G4QKaonPlusElasticCrossSection::GetTabV:TM="<<lastTM<<",S1="<<theS1<<",B1="
01004 <<theB1<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS1<<",B3="<<theB1<<G4endl;
01005 #endif
01006
01007 G4double dp=lp-lastPAR[4];
01008
01009 return lastPAR[0]/(lastPAR[2]+sqr(p-lastPAR[1]))+(lastPAR[3]*dp*dp+lastPAR[5])/
01010 (1.-lastPAR[6]/sp+lastPAR[7]/p4)
01011 +lastPAR[8]/(sqr(p-lastPAR[9])+lastPAR[10]);
01012
01013 }
01014 else
01015 {
01016 G4double p5=p4*p;
01017 G4double p6=p5*p;
01018 G4double p8=p6*p2;
01019 G4double p10=p8*p2;
01020 G4double p12=p10*p2;
01021 G4double p16=p8*p8;
01022
01023 G4double dl=lp-5.;
01024 G4double a=tgZ+tgN;
01025 G4double pah=std::pow(p,a/2);
01026 G4double pa=pah*pah;
01027 G4double pa2=pa*pa;
01028 if(a<6.5)
01029 {
01030 theS1=lastPAR[9]/(1.+lastPAR[10]*p4*pa)+lastPAR[11]/(p4+lastPAR[12]*p4/pa2)+
01031 (lastPAR[13]*dl*dl+lastPAR[14])/(1.+lastPAR[15]/p2);
01032 theB1=(lastPAR[16]+lastPAR[17]*p2)/(p4+lastPAR[18]/pah)+lastPAR[19];
01033 theSS=lastPAR[20]/(1.+lastPAR[21]/p2)+lastPAR[22]/(p6/pa+lastPAR[23]/p16);
01034 theS2=lastPAR[24]/(pa/p2+lastPAR[25]/p4)+lastPAR[26];
01035 theB2=lastPAR[27]*std::pow(p,lastPAR[28])+lastPAR[29]/(p8+lastPAR[30]/p16);
01036 theS3=lastPAR[31]/(pa*p+lastPAR[32]/pa)+lastPAR[33];
01037 theB3=lastPAR[34]/(p3+lastPAR[35]/p6)+lastPAR[36]/(1.+lastPAR[37]/p2);
01038 theS4=p2*(pah*lastPAR[38]*std::exp(-pah*lastPAR[39])+
01039 lastPAR[40]/(1.+lastPAR[41]*std::pow(p,lastPAR[42])));
01040 theB4=lastPAR[43]*pa/p2/(1.+pa*lastPAR[44]);
01041 #ifdef tdebug
01042 G4cout<<"G4QKaonPlusElasticCS::GetTabV: lA, p="<<p<<",S1="<<theS1<<",B1="<<theB1
01043 <<",SS="<<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3
01044 <<",S4="<<theS4<<",B4="<<theB4<<G4endl;
01045 #endif
01046 }
01047 else
01048 {
01049 theS1=lastPAR[9]/(1.+lastPAR[10]/p4)+lastPAR[11]/(p4+lastPAR[12]/p2)+
01050 lastPAR[13]/(p5+lastPAR[14]/p16);
01051 theB1=(lastPAR[15]/p8+lastPAR[19])/(p+lastPAR[16]/std::pow(p,lastPAR[20]))+
01052 lastPAR[17]/(1.+lastPAR[18]/p4);
01053 theSS=lastPAR[21]/(p4/std::pow(p,lastPAR[23])+lastPAR[22]/p4);
01054 theS2=lastPAR[24]/p4/(std::pow(p,lastPAR[25])+lastPAR[26]/p12)+lastPAR[27];
01055 theB2=lastPAR[28]/std::pow(p,lastPAR[29])+lastPAR[30]/std::pow(p,lastPAR[31]);
01056 theS3=lastPAR[32]/std::pow(p,lastPAR[35])/(1.+lastPAR[36]/p12)+
01057 lastPAR[33]/(1.+lastPAR[34]/p6);
01058 theB3=lastPAR[37]/p8+lastPAR[38]/p2+lastPAR[39]/(1.+lastPAR[40]/p8);
01059 theS4=(lastPAR[41]/p4+lastPAR[46]/p)/(1.+lastPAR[42]/p10)+
01060 (lastPAR[43]+lastPAR[44]*dl*dl)/(1.+lastPAR[45]/p12);
01061 theB4=lastPAR[47]/(1.+lastPAR[48]/p)+lastPAR[49]*p4/(1.+lastPAR[50]*p5);
01062 #ifdef tdebug
01063 G4cout<<"G4QKaonPlusElasticCS::GetTabVal:hA,p="<<p<<",S1="<<theS1<<",B1="<<theB1
01064 <<",SS="<<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3
01065 <<",S4="<<theS4<<",B4="<<theB4<<G4endl;
01066 #endif
01067 }
01068
01069 #ifdef tdebug
01070 G4cout<<"G4QKaonPlusElCS::GTV: PDG="<<PDG<<",P="<<p<<",N="<<tgN<<",Z="<<tgZ<<G4endl;
01071 #endif
01072 G4double dlp=lp-lastPAR[4];
01073
01074 return (lastPAR[0]*dlp*dlp+lastPAR[1]+lastPAR[2]/p2)/(1.+lastPAR[3]/p2/sp);
01075 }
01076 return 0.;
01077 }
01078
01079
01080 G4double G4QKaonPlusElasticCrossSection::GetQ2max(G4int PDG, G4int tgZ, G4int tgN,
01081 G4double pP)
01082 {
01083
01084
01085 static const G4double mK= G4QPDGCode(321).GetMass()*.001;
01086
01087
01088
01089
01090
01091 static const G4double mK2= mK*mK;
01092
01093
01094
01095
01096 G4double pP2=pP*pP;
01097 if(tgZ || tgN>-1)
01098 {
01099 G4double mt=G4QPDGCode(90000000+tgZ*1000+tgN).GetMass()*.001;
01100 G4double dmt=mt+mt;
01101 G4double s_value=dmt*std::sqrt(pP2+mK2)+mK2+mt*mt;
01102 return dmt*dmt*pP2/s_value;
01103 }
01104 else
01105 {
01106
01107
01108
01109 G4ExceptionDescription ed;
01110 ed << "PDG = " << PDG << ",Z = " << tgZ << ", N = " << tgN
01111 << ", while it is defined only for p projectiles & Z_target>0" << G4endl;
01112 G4Exception("G4QKaonPlusElasticCrossSection::GetQ2max()", "HAD_CHPS_0000",
01113 FatalException, ed);
01114 return 0;
01115 }
01116 }