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 "G4QHyperonElasticCrossSection.hh"
00046 #include "G4SystemOfUnits.hh"
00047
00048
00049 const G4int G4QHyperonElasticCrossSection::nPoints=128;
00050 const G4int G4QHyperonElasticCrossSection::nLast=nPoints-1;
00051 G4double G4QHyperonElasticCrossSection::lPMin=-8.;
00052 G4double G4QHyperonElasticCrossSection::lPMax= 8.;
00053 G4double G4QHyperonElasticCrossSection::dlnP=(lPMax-lPMin)/nLast;
00054 G4bool G4QHyperonElasticCrossSection::onlyCS=true;
00055 G4double G4QHyperonElasticCrossSection::lastSIG=0.;
00056 G4double G4QHyperonElasticCrossSection::lastLP=-10.;
00057 G4double G4QHyperonElasticCrossSection::lastTM=0.;
00058 G4double G4QHyperonElasticCrossSection::theSS=0.;
00059 G4double G4QHyperonElasticCrossSection::theS1=0.;
00060 G4double G4QHyperonElasticCrossSection::theB1=0.;
00061 G4double G4QHyperonElasticCrossSection::theS2=0.;
00062 G4double G4QHyperonElasticCrossSection::theB2=0.;
00063 G4double G4QHyperonElasticCrossSection::theS3=0.;
00064 G4double G4QHyperonElasticCrossSection::theB3=0.;
00065 G4double G4QHyperonElasticCrossSection::theS4=0.;
00066 G4double G4QHyperonElasticCrossSection::theB4=0.;
00067 G4int G4QHyperonElasticCrossSection::lastTZ=0;
00068 G4int G4QHyperonElasticCrossSection::lastTN=0;
00069 G4double G4QHyperonElasticCrossSection::lastPIN=0.;
00070 G4double* G4QHyperonElasticCrossSection::lastCST=0;
00071 G4double* G4QHyperonElasticCrossSection::lastPAR=0;
00072 G4double* G4QHyperonElasticCrossSection::lastSST=0;
00073 G4double* G4QHyperonElasticCrossSection::lastS1T=0;
00074 G4double* G4QHyperonElasticCrossSection::lastB1T=0;
00075 G4double* G4QHyperonElasticCrossSection::lastS2T=0;
00076 G4double* G4QHyperonElasticCrossSection::lastB2T=0;
00077 G4double* G4QHyperonElasticCrossSection::lastS3T=0;
00078 G4double* G4QHyperonElasticCrossSection::lastB3T=0;
00079 G4double* G4QHyperonElasticCrossSection::lastS4T=0;
00080 G4double* G4QHyperonElasticCrossSection::lastB4T=0;
00081 G4int G4QHyperonElasticCrossSection::lastN=0;
00082 G4int G4QHyperonElasticCrossSection::lastZ=0;
00083 G4double G4QHyperonElasticCrossSection::lastP=0.;
00084 G4double G4QHyperonElasticCrossSection::lastTH=0.;
00085 G4double G4QHyperonElasticCrossSection::lastCS=0.;
00086 G4int G4QHyperonElasticCrossSection::lastI=0;
00087
00088 std::vector<G4double*> G4QHyperonElasticCrossSection::PAR;
00089 std::vector<G4double*> G4QHyperonElasticCrossSection::CST;
00090 std::vector<G4double*> G4QHyperonElasticCrossSection::SST;
00091 std::vector<G4double*> G4QHyperonElasticCrossSection::S1T;
00092 std::vector<G4double*> G4QHyperonElasticCrossSection::B1T;
00093 std::vector<G4double*> G4QHyperonElasticCrossSection::S2T;
00094 std::vector<G4double*> G4QHyperonElasticCrossSection::B2T;
00095 std::vector<G4double*> G4QHyperonElasticCrossSection::S3T;
00096 std::vector<G4double*> G4QHyperonElasticCrossSection::B3T;
00097 std::vector<G4double*> G4QHyperonElasticCrossSection::S4T;
00098 std::vector<G4double*> G4QHyperonElasticCrossSection::B4T;
00099
00100 G4QHyperonElasticCrossSection::G4QHyperonElasticCrossSection()
00101 {
00102 }
00103
00104 G4QHyperonElasticCrossSection::~G4QHyperonElasticCrossSection()
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* G4QHyperonElasticCrossSection::GetPointer()
00144 {
00145 static G4QHyperonElasticCrossSection theCrossSection;
00146 return &theCrossSection;
00147 }
00148
00149
00150
00151 G4double G4QHyperonElasticCrossSection::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<<"G4QHyperElCS::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==3222 || pPDG<3000 || pPDG>3334)
00168 {
00169 G4cout<<"*Warning*G4QHyperonElaCS::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<<"G4QHyElCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",i="<<i<<G4endl;
00186
00187 #endif
00188 if(pEn<=lastTH)
00189 {
00190 #ifdef pdebug
00191 G4cout<<"G4QHyElCS::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<<"G4QHyperonElasticCS::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<<"G4QHyElCS::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<<"G4QHyperElCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
00215
00216 #endif
00217 if(lastCS<=0. && pEn>lastTH)
00218 {
00219 #ifdef pdebug
00220 G4cout<<"G4QHyperonElCS::GetCS: NewT="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00221 #endif
00222 lastTH=pEn;
00223 }
00224 break;
00225 }
00226 #ifdef pdebug
00227 G4cout<<"---G4QHyperonElCrossSection::GetCrosSec:pPDG="<<pPDG<<",i="<<i<<",N="<<colN[i]
00228 <<",Z["<<i<<"]="<<colZ[i]<<G4endl;
00229
00230 #endif
00231 }
00232 if(!in)
00233 {
00234 #ifdef pdebug
00235 G4cout<<"G4QHyElCS::GetCrosSec:CalcNew P="<<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<<"G4QHyperonElasticCrossSection::GetCS:NewThres="<<lastTH<<",T="<<pEn<<G4endl;
00244 #endif
00245 if(pEn>lastTH)
00246 {
00247 #ifdef pdebug
00248 G4cout<<"G4QHyperonElCS::GetCS:1st T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00249 #endif
00250 lastTH=pEn;
00251 }
00252 }
00253 #ifdef pdebug
00254 G4cout<<"G4QHyElCS::GetCrosSec: New CS="<<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<<"G4QHyElCS::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<<"G4QHyperonElasticCrossSection::GetCS: Update lastI="<<lastI<<G4endl;
00272 #endif
00273 colP[lastI]=pMom;
00274 colCS[lastI]=lastCS;
00275 }
00276 #ifdef pdebug
00277 G4cout<<"G4QHyElCS::GetCSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00278
00279 G4cout<<"G4QHyperonElasticCrossSection::GetCrosSec:***End***,onlyCS="<<onlyCS<<G4endl;
00280 #endif
00281 return lastCS*millibarn;
00282 }
00283
00284
00285
00286 G4double G4QHyperonElasticCrossSection::CalculateCrossSection(G4bool CS,G4int F,G4int I,
00287 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<<"G4QHyperonElasticCS::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<<"G4QHyperonElCS::CalcCS: DB's updated for I="<<I<<",*,PIN4="<<PIN[4]<<G4endl;
00316 #endif
00317 }
00318 #ifdef pdebug
00319 G4cout<<"G4QHyperonElasticCS::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<<"G4QHyElCS::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<<"G4QHyperonElasticCrosS::CalcCS:*ini*,lastLP="<<lastLP<<",min="<<lPMin<<G4endl;
00346 #endif
00347 lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN);
00348 #ifdef pdebug
00349 G4cout<<"G4QHypElCS::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<<"G4QHypElCS::CalCS:?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<<"G4QHyperionElCS::CalcCS:*updated(O)*, LP="<<lastLP<<" < IN="<<lastPIN<<G4endl;
00373 #endif
00374 }
00375 #ifdef pdebug
00376 G4cout<<"G4QHypElCS::CalcCS: lastLP="<<lastLP<<",lPM="<<lPMin<<",lPIN="<<lastPIN<<G4endl;
00377 #endif
00378 if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom);
00379 #ifdef pdebug
00380 G4cout<<"G4QHyperElCrosSec::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<<"G4QHyperElCS::CCS:b="<<blast<<","<<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<<"G4QHyperonElasticCS::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<<"G4QHyperonElasticCrossSection::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<<"G4QHyperonElasticCS::CalcCrossSection: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<<"G4QHyperonElasticCrossSection::CCS:s3l="<<S3TL<<",sh3="<<shift<<",s3h="
00440 <<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<<"G4QHyperonElasticCrossSection::CCS:s4l="<<S4TL<<",sh4="<<shift<<",s4h="
00448 <<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<<"G4QHyperonElasticCS::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<<"G4QHyperonElasticCrossSection::CalculateCS: END, onlyCS="<<onlyCS<<G4endl;
00462 #endif
00463 return lastSIG;
00464 }
00465
00466
00467 G4double G4QHyperonElasticCrossSection::GetPTables(G4double LP, G4double ILP, G4int PDG,
00468 G4int tgZ, G4int tgN)
00469 {
00470
00471 static const G4double pwd=2727;
00472 const G4int n_hypel=33;
00473
00474 G4double hyp_el[n_hypel]={1.,.002,.12,.0557,3.5,6.72,99.,2.,3.,5.,74.,3.,3.4,.2,.17,
00475 .001,8.,.055,3.64,5.e-5,4000.,1500.,.46,1.2e6,3.5e6,5.e-5,
00476 1.e10,8.5e8,1.e10,1.1,3.4e6,6.8e6,0.};
00477
00478
00479 if(PDG!=3222 && PDG>3000 && PDG<3335)
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_hypel; ip++) lastPAR[ip]=hyp_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]=4./(1.+22/asa);
00521 lastPAR[1]=2.36*asa/(1.+a*.055/ssa);
00522 lastPAR[2]=(1.+.00007*a3/ssa)/(1.+.0026*a2);
00523 lastPAR[3]=1.76*a/ssa+.00003*a3;
00524 lastPAR[4]=(.03+200./a3)/(1.+1.E5/a3/sa);
00525 lastPAR[5]=5.;
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<<"G4QHyElCS::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<<"G4QHyElCS::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<<"G4QHyperonElasticCrossSection::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<<"G4QHyperonElasticCrossSection::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*G4QHyperonElasticCrossSection::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*G4QHyperonElasticCrossSection::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*G4QHyperonElasticCrossSection::GetPTab:PDG="<<PDG<<",Z="<<tgZ
00716 <<", N="<<tgN<<", LP="<<LP<<" <= ILP="<<ILP<<" nothing is done!"<<G4endl;
00717 #endif
00718 } else {
00719
00720
00721
00722 G4ExceptionDescription ed;
00723 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
00724 << ", while it is defined only for Hyperons" << G4endl;
00725 G4Exception("G4QHyperonElasticCrossSection::GetPTables()", "HAD_CHPS_0000",
00726 FatalException, ed);
00727 }
00728 return ILP;
00729 }
00730
00731
00732 G4double G4QHyperonElasticCrossSection::GetExchangeT(G4int tgZ, G4int tgN, G4int PDG)
00733 {
00734 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
00735 static const G4double third=1./3.;
00736 static const G4double fifth=1./5.;
00737 static const G4double sevth=1./7.;
00738 #ifdef tdebug
00739 G4cout<<"G4QHyperElCS::GetExcT:F="<<onlyCS<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
00740 #endif
00741 if(PDG==3222 || PDG<3000 || PDG>3334)G4cout<<"*Warning*G4QHyElCS::GET:PDG="<<PDG<<G4endl;
00742 if(onlyCS)G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExchanT: onlyCS=1"<<G4endl;
00743 if(lastLP<-4.3) return lastTM*GeVSQ*G4UniformRand();
00744 G4double q2=0.;
00745 if(tgZ==1 && tgN==0)
00746 {
00747 #ifdef tdebug
00748 G4cout<<"G4QHyperElCS::GetExchangeT: TM="<<lastTM<<",S1="<<theS1<<",B1="<<theB1<<",S2="
00749 <<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3<<",GeV2="<<GeVSQ<<G4endl;
00750 #endif
00751 G4double E1=lastTM*theB1;
00752 G4double R1=(1.-std::exp(-E1));
00753 #ifdef tdebug
00754 G4double ts1=-std::log(1.-R1)/theB1;
00755 G4double ds1=std::fabs(ts1-lastTM)/lastTM;
00756 if(ds1>.0001)
00757 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:1p "<<ts1<<"#"<<lastTM
00758 <<",d="<<ds1<<",R1="<<R1<<",E1="<<E1<<G4endl;
00759 #endif
00760 G4double E2=lastTM*theB2;
00761 G4double R2=(1.-std::exp(-E2*E2*E2));
00762 #ifdef tdebug
00763 G4double ts2=std::pow(-std::log(1.-R2),.333333333)/theB2;
00764 G4double ds2=std::fabs(ts2-lastTM)/lastTM;
00765 if(ds2>.0001)
00766 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:2p "<<ts2<<"#"<<lastTM
00767 <<",d="<<ds2<<",R2="<<R2<<",E2="<<E2<<G4endl;
00768 #endif
00769 G4double E3=lastTM*theB3;
00770 G4double R3=(1.-std::exp(-E3));
00771 #ifdef tdebug
00772 G4double ts3=-std::log(1.-R3)/theB3;
00773 G4double ds3=std::fabs(ts3-lastTM)/lastTM;
00774 if(ds3>.0001)
00775 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:3p "<<ts3<<"#"<<lastTM
00776 <<",d="<<ds3<<",R3="<<R1<<",E3="<<E3<<G4endl;
00777 #endif
00778 G4double I1=R1*theS1/theB1;
00779 G4double I2=R2*theS2;
00780 G4double I3=R3*theS3;
00781 G4double I12=I1+I2;
00782 G4double rand=(I12+I3)*G4UniformRand();
00783 if (rand<I1 )
00784 {
00785 G4double ran=R1*G4UniformRand();
00786 if(ran>1.) ran=1.;
00787 q2=-std::log(1.-ran)/theB1;
00788 }
00789 else if(rand<I12)
00790 {
00791 G4double ran=R2*G4UniformRand();
00792 if(ran>1.) ran=1.;
00793 q2=-std::log(1.-ran);
00794 if(q2<0.) q2=0.;
00795 q2=std::pow(q2,third)/theB2;
00796 }
00797 else
00798 {
00799 G4double ran=R3*G4UniformRand();
00800 if(ran>1.) ran=1.;
00801 q2=-std::log(1.-ran)/theB3;
00802 }
00803 }
00804 else
00805 {
00806 G4double a=tgZ+tgN;
00807 #ifdef tdebug
00808 G4cout<<"G4QHyperonElasticCrossSection::GetExT:a="<<a<<",t="<<lastTM<<",S1="<<theS1
00809 <<",B1="<<theB1<<",SS="<<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3
00810 <<",B3="<<theB3<<",S4="<<theS4<<",B4="<<theB4<<G4endl;
00811 #endif
00812 G4double E1=lastTM*(theB1+lastTM*theSS);
00813 G4double R1=(1.-std::exp(-E1));
00814 G4double tss=theSS+theSS;
00815 #ifdef tdebug
00816 G4double ts1=-std::log(1.-R1)/theB1;
00817 if(std::fabs(tss)>1.e-7) ts1=(std::sqrt(theB1*(theB1+(tss+tss)*ts1))-theB1)/tss;
00818 G4double ds1=(ts1-lastTM)/lastTM;
00819 if(ds1>.0001)
00820 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:1a "<<ts1<<"#"<<lastTM
00821 <<",d="<<ds1<<",R1="<<R1<<",E1="<<E1<<G4endl;
00822 #endif
00823 G4double tm2=lastTM*lastTM;
00824 G4double E2=lastTM*tm2*theB2;
00825 if(a>6.5)E2*=tm2;
00826 G4double R2=(1.-std::exp(-E2));
00827 #ifdef tdebug
00828 G4double ts2=-std::log(1.-R2)/theB2;
00829 if(a<6.5)ts2=std::pow(ts2,third);
00830 else ts2=std::pow(ts2,fifth);
00831 G4double ds2=std::fabs(ts2-lastTM)/lastTM;
00832 if(ds2>.0001)
00833 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:2a "<<ts2<<"#"<<lastTM
00834 <<",d="<<ds2<<",R2="<<R2<<",E2="<<E2<<G4endl;
00835 #endif
00836 G4double E3=lastTM*theB3;
00837 if(a>6.5)E3*=tm2*tm2*tm2;
00838 G4double R3=(1.-std::exp(-E3));
00839 #ifdef tdebug
00840 G4double ts3=-std::log(1.-R3)/theB3;
00841 if(a>6.5)ts3=std::pow(ts3,sevth);
00842 G4double ds3=std::fabs(ts3-lastTM)/lastTM;
00843 if(ds3>.0001)
00844 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:3a "<<ts3<<"#"<<lastTM
00845 <<",d="<<ds3<<",R3="<<R3<<",E3="<<E3<<G4endl;
00846 #endif
00847 G4double E4=lastTM*theB4;
00848 G4double R4=(1.-std::exp(-E4));
00849 #ifdef tdebug
00850 G4double ts4=-std::log(1.-R4)/theB4;
00851 G4double ds4=std::fabs(ts4-lastTM)/lastTM;
00852 if(ds4>.0001)
00853 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetExT:4a "<<ts4<<"#"<<lastTM
00854 <<",d="<<ds4<<",R4="<<R4<<",E4="<<E4<<G4endl;
00855 #endif
00856 G4double I1=R1*theS1;
00857 G4double I2=R2*theS2;
00858 G4double I3=R3*theS3;
00859 G4double I4=R4*theS4;
00860 G4double I12=I1+I2;
00861 G4double I13=I12+I3;
00862 G4double rand=(I13+I4)*G4UniformRand();
00863 #ifdef tdebug
00864 G4cout<<"G4QHyElCS::GExT:1="<<I1<<",2="<<I2<<",3="<<I3<<",4="<<I4<<",r="<<rand<<G4endl;
00865 #endif
00866 if(rand<I1)
00867 {
00868 G4double ran=R1*G4UniformRand();
00869 if(ran>1.) ran=1.;
00870 q2=-std::log(1.-ran)/theB1;
00871 if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
00872 #ifdef tdebug
00873 G4cout<<"G4QHyElCS::GExT:Q2="<<q2<<",ss="<<tss/2<<",b1="<<theB1<<",t1="<<ts1<<G4endl;
00874 #endif
00875 }
00876 else if(rand<I12)
00877 {
00878 G4double ran=R2*G4UniformRand();
00879 if(ran>1.) ran=1.;
00880 q2=-std::log(1.-ran)/theB2;
00881 if(q2<0.) q2=0.;
00882 if(a<6.5) q2=std::pow(q2,third);
00883 else q2=std::pow(q2,fifth);
00884 #ifdef tdebug
00885 G4cout<<"G4QHyElCS::GetExT: Q2="<<q2<<",r2="<<R2<<",b2="<<theB2<<",t2="<<ts2<<G4endl;
00886 #endif
00887 }
00888 else if(rand<I13)
00889 {
00890 G4double ran=R3*G4UniformRand();
00891 if(ran>1.) ran=1.;
00892 q2=-std::log(1.-ran)/theB3;
00893 if(q2<0.) q2=0.;
00894 if(a>6.5) q2=std::pow(q2,sevth);
00895 #ifdef tdebug
00896 G4cout<<"G4QHyElCS::GetExT: Q2="<<q2<<",r3="<<R2<<",b3="<<theB2<<",t3="<<ts2<<G4endl;
00897 #endif
00898 }
00899 else
00900 {
00901 G4double ran=R4*G4UniformRand();
00902 if(ran>1.) ran=1.;
00903 q2=-std::log(1.-ran)/theB4;
00904 if(a<6.5) q2=lastTM-q2;
00905 #ifdef tdebug
00906 G4cout<<"G4QHyElCS::GExT:Q2="<<q2<<",m="<<lastTM<<",b4="<<theB3<<",t4="<<ts3<<G4endl;
00907 #endif
00908 }
00909 }
00910 if(q2<0.) q2=0.;
00911 if(!(q2>=-1.||q2<=1.))G4cout<<"*NAN*G4QHyElasticCrossSect::GetExchangeT:-t="<<q2<<G4endl;
00912 if(q2>lastTM)
00913 {
00914 #ifdef tdebug
00915 G4cout<<"*Warning*G4QHyperonElasticCrossSection::GExT: -t="<<q2<<" > "<<lastTM<<G4endl;
00916 #endif
00917 q2=lastTM;
00918 }
00919 return q2*GeVSQ;
00920 }
00921
00922
00923 G4double G4QHyperonElasticCrossSection::GetSlope(G4int tgZ, G4int tgN, G4int PDG)
00924 {
00925 static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
00926 #ifdef tdebug
00927 G4cout<<"G4QHyperElCS::GetSlope:"<<onlyCS<<", Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
00928 #endif
00929 if(onlyCS)G4cout<<"*Warning*G4QHyperonElasticCrossSection::GetSlope: onlCS=true"<<G4endl;
00930 if(lastLP<-4.3) return 0.;
00931 if(PDG==3222 || PDG<3000 || PDG>3334)
00932 {
00933
00934
00935
00936 G4ExceptionDescription ed;
00937 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
00938 << ", while it is defined only for Hyperons" << G4endl;
00939 G4Exception("G4QHyperonElasticCrossSection::GetSlope()", "HAD_CHPS_0000",
00940 FatalException, ed);
00941 }
00942 if(theB1<0.) theB1=0.;
00943 if(!(theB1>=-1.||theB1<=1.)) G4cout<<"*NAN*G4QHyElasticCrossS::Getslope:"<<theB1<<G4endl;
00944 return theB1/GeVSQ;
00945 }
00946
00947
00948 G4double G4QHyperonElasticCrossSection::GetHMaxT()
00949 {
00950 static const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
00951 return lastTM*HGeVSQ;
00952 }
00953
00954
00955 G4double G4QHyperonElasticCrossSection::GetTabValues(G4double lp, G4int PDG, G4int tgZ,
00956 G4int tgN)
00957 {
00958 if(PDG==3222 || PDG<3000 || PDG>3334) G4cout<<"*Warning*G4QHypElCS::GTV:P="<<PDG<<G4endl;
00959 if(tgZ<0 || tgZ>92)
00960 {
00961 G4cout<<"*Warning*G4QHyperonElastCS::GetTabValue:(1-92) NoIsotopesFor Z="<<tgZ<<G4endl;
00962 return 0.;
00963 }
00964 G4int iZ=tgZ-1;
00965 if(iZ<0)
00966 {
00967 iZ=0;
00968 tgZ=1;
00969 tgN=0;
00970 }
00971
00972
00973 #ifdef isodebug
00974
00975 #endif
00976
00977
00978 #ifdef pdebug
00979 G4cout<<"G4QHyElasticCS::GetTabVal:l="<<lp<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
00980 #endif
00981 G4double p=std::exp(lp);
00982 G4double sp=std::sqrt(p);
00983 G4double p2=p*p;
00984 G4double p3=p2*p;
00985 G4double p4=p3*p;
00986 if ( tgZ == 1 && tgN == 0 )
00987 {
00988 G4double dl2=lp-lastPAR[9];
00989 theSS=lastPAR[32];
00990 theS1=(lastPAR[10]+lastPAR[11]*dl2*dl2)/(1.+lastPAR[12]/p4/p)+
00991 (lastPAR[13]/p2+lastPAR[14]*p)/(p4+lastPAR[15]*sp);
00992 theB1=lastPAR[16]*std::pow(p,lastPAR[17])/(1.+lastPAR[18]/p3);
00993 theS2=lastPAR[19]+lastPAR[20]/(p4+lastPAR[21]*p);
00994 theB2=lastPAR[22]+lastPAR[23]/(p4+lastPAR[24]/sp);
00995 theS3=lastPAR[25]+lastPAR[26]/(p4*p4+lastPAR[27]*p2+lastPAR[28]);
00996 theB3=lastPAR[29]+lastPAR[30]/(p4+lastPAR[31]);
00997 theS4=0.;
00998 theB4=0.;
00999 #ifdef tdebug
01000 G4cout<<"G4QHyperonElasticCrossSection::GetTableVal:TM="<<lastTM<<",S1="<<theS1<<",B1="
01001 <<theB1<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS1<<",B3="<<theB1<<G4endl;
01002 #endif
01003
01004 G4double dp=lp-lastPAR[4];
01005 return lastPAR[0]/(lastPAR[1]+p2*(lastPAR[2]+p2))+(lastPAR[3]*dp*dp+lastPAR[5]+
01006 lastPAR[6]/p2)/(1.+lastPAR[7]/sp+lastPAR[8]/p4);
01007 }
01008 else
01009 {
01010 G4double p5=p4*p;
01011 G4double p6=p5*p;
01012 G4double p8=p6*p2;
01013 G4double p10=p8*p2;
01014 G4double p12=p10*p2;
01015 G4double p16=p8*p8;
01016
01017 G4double dl=lp-5.;
01018 G4double a=tgZ+tgN;
01019 G4double pah=std::pow(p,a/2);
01020 G4double pa=pah*pah;
01021 G4double pa2=pa*pa;
01022 if(a<6.5)
01023 {
01024 theS1=lastPAR[9]/(1.+lastPAR[10]*p4*pa)+lastPAR[11]/(p4+lastPAR[12]*p4/pa2)+
01025 (lastPAR[13]*dl*dl+lastPAR[14])/(1.+lastPAR[15]/p2);
01026 theB1=(lastPAR[16]+lastPAR[17]*p2)/(p4+lastPAR[18]/pah)+lastPAR[19];
01027 theSS=lastPAR[20]/(1.+lastPAR[21]/p2)+lastPAR[22]/(p6/pa+lastPAR[23]/p16);
01028 theS2=lastPAR[24]/(pa/p2+lastPAR[25]/p4)+lastPAR[26];
01029 theB2=lastPAR[27]*std::pow(p,lastPAR[28])+lastPAR[29]/(p8+lastPAR[30]/p16);
01030 theS3=lastPAR[31]/(pa*p+lastPAR[32]/pa)+lastPAR[33];
01031 theB3=lastPAR[34]/(p3+lastPAR[35]/p6)+lastPAR[36]/(1.+lastPAR[37]/p2);
01032 theS4=p2*(pah*lastPAR[38]*std::exp(-pah*lastPAR[39])+
01033 lastPAR[40]/(1.+lastPAR[41]*std::pow(p,lastPAR[42])));
01034 theB4=lastPAR[43]*pa/p2/(1.+pa*lastPAR[44]);
01035 #ifdef tdebug
01036 G4cout<<"G4QHyperonElasticCS::GetTabV: lA, p="<<p<<",S1="<<theS1<<",B1="<<theB1
01037 <<",SS="<<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3
01038 <<",S4="<<theS4<<",B4="<<theB4<<G4endl;
01039 #endif
01040 }
01041 else
01042 {
01043 theS1=lastPAR[9]/(1.+lastPAR[10]/p4)+lastPAR[11]/(p4+lastPAR[12]/p2)+
01044 lastPAR[13]/(p5+lastPAR[14]/p16);
01045 theB1=(lastPAR[15]/p8+lastPAR[19])/(p+lastPAR[16]/std::pow(p,lastPAR[20]))+
01046 lastPAR[17]/(1.+lastPAR[18]/p4);
01047 theSS=lastPAR[21]/(p4/std::pow(p,lastPAR[23])+lastPAR[22]/p4);
01048 theS2=lastPAR[24]/p4/(std::pow(p,lastPAR[25])+lastPAR[26]/p12)+lastPAR[27];
01049 theB2=lastPAR[28]/std::pow(p,lastPAR[29])+lastPAR[30]/std::pow(p,lastPAR[31]);
01050 theS3=lastPAR[32]/std::pow(p,lastPAR[35])/(1.+lastPAR[36]/p12)+
01051 lastPAR[33]/(1.+lastPAR[34]/p6);
01052 theB3=lastPAR[37]/p8+lastPAR[38]/p2+lastPAR[39]/(1.+lastPAR[40]/p8);
01053 theS4=(lastPAR[41]/p4+lastPAR[46]/p)/(1.+lastPAR[42]/p10)+
01054 (lastPAR[43]+lastPAR[44]*dl*dl)/(1.+lastPAR[45]/p12);
01055 theB4=lastPAR[47]/(1.+lastPAR[48]/p)+lastPAR[49]*p4/(1.+lastPAR[50]*p5);
01056 #ifdef tdebug
01057 G4cout<<"G4QHyperonElasticCS::GetTabV:hA,p="<<p<<",S1="<<theS1<<",B1="<<theB1<<",SS="
01058 <<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3<<",S4="
01059 <<theS4<<",B4="<<theB4<<G4endl;
01060 #endif
01061 }
01062
01063 #ifdef tdebug
01064 G4cout<<"G4QHyperonElCS::GetTabV: PDG="<<PDG<<",P="<<p<<",N="<<tgN<<",Z="<<tgZ<<G4endl;
01065 #endif
01066 G4double dlp=lp-lastPAR[5];
01067
01068 return (lastPAR[0]*dlp*dlp+lastPAR[1])/(1.+lastPAR[2]/p)+lastPAR[3]/(p3+lastPAR[4]);
01069 }
01070 return 0.;
01071 }
01072
01073
01074 G4double G4QHyperonElasticCrossSection::GetQ2max(G4int PDG, G4int tgZ, G4int tgN,
01075 G4double pP)
01076 {
01077
01078
01079
01080 static const G4double mLamb= G4QPDGCode(3122).GetMass()*.001;
01081
01082
01083
01084
01085 static const G4double mLa2= mLamb*mLamb;
01086
01087
01088
01089 G4double pP2=pP*pP;
01090 if(tgZ || tgN>-1)
01091 {
01092 G4double mt=G4QPDGCode(90000000+tgZ*1000+tgN).GetMass()*.001;
01093 G4double dmt=mt+mt;
01094 G4double s_value=dmt*std::sqrt(pP2+mLa2)+mLa2+mt*mt;
01095 return dmt*dmt*pP2/s_value;
01096 }
01097 else
01098 {
01099
01100
01101
01102 G4ExceptionDescription ed;
01103 ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
01104 << ", while it is defined only for p projectiles & Z_target>0" << G4endl;
01105 G4Exception("G4QHyperonElasticCrossSection::GetQ2max()", "HAD_CHPS_0000",
01106 FatalException, ed);
01107 return 0;
01108 }
01109 }