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