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 #include "G4QFreeScattering.hh"
00044 #include "G4SystemOfUnits.hh"
00045
00046
00047 std::vector<std::pair<G4double,G4double>*> G4QFreeScattering::vX;
00048
00049 G4QFreeScattering::G4QFreeScattering()
00050 {
00051 #ifdef pdebug
00052 G4cout<<"***^^^*** G4QFreeScattering singletone is created ***^^^***"<<G4endl;
00053 #endif
00054 }
00055
00056 G4QFreeScattering::~G4QFreeScattering()
00057 {
00058 std::vector<std::pair<G4double,G4double>*>::iterator pos2;
00059 for(pos2=vX.begin(); pos2<vX.end(); pos2++) delete [] * pos2;
00060 vX.clear();
00061 }
00062
00063
00064 G4QFreeScattering* G4QFreeScattering::GetPointer()
00065 {
00066 static G4QFreeScattering theQFS;
00067 return &theQFS;
00068 }
00069
00070
00071 std::pair<G4double,G4double> G4QFreeScattering::CalcElTot(G4double p, G4int I)
00072 {
00073
00074 static const G4double lmi=3.5;
00075 static const G4double pbe=.0557;
00076 static const G4double pbt=.3;
00077 static const G4double pmi=.1;
00078 static const G4double pma=1000.;
00079 G4double El=0.;
00080 G4double To=0.;
00081 if(p<=0.)
00082 {
00083
00084 return std::make_pair(El,To);
00085 }
00086 if (!I)
00087 {
00088 #ifdef debug
00089 G4cout<<"G4QFreeScatter::CalcElTot:I=0, p="<<p<<", pmi="<<pmi<<", pma="<<pma<<G4endl;
00090 #endif
00091 if(p<pmi)
00092 {
00093 G4double p2=p*p;
00094 El=1./(.00012+p2*.2);
00095 To=El;
00096 #ifdef debug
00097 G4cout<<"G4QFreeScatter::CalcElTot:I=0i, El="<<El<<", To="<<To<<", p2="<<p2<<G4endl;
00098 #endif
00099 }
00100 else if(p>pma)
00101 {
00102 G4double lp=std::log(p)-lmi;
00103 G4double lp2=lp*lp;
00104 El=pbe*lp2+6.72;
00105 To=pbt*lp2+38.2;
00106 #ifdef debug
00107 G4cout<<"G4QFreeScat::CalcElTot:I=0a, El="<<El<<", To="<<To<<", lp2="<<lp2<<G4endl;
00108 #endif
00109 }
00110 else
00111 {
00112 G4double p2=p*p;
00113 G4double LE=1./(.00012+p2*.2);
00114 G4double lp=std::log(p)-lmi;
00115 G4double lp2=lp*lp;
00116 G4double rp2=1./p2;
00117 El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p);
00118 To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2);
00119 #ifdef debug
00120 G4cout<<"G4QFreeScat::CalcElTot:0,E="<<El<<",T="<<To<<",s="<<p2<<",l="<<lp2<<G4endl;
00121 #endif
00122 }
00123 }
00124 else if(I==1)
00125 {
00126 if(p<pmi)
00127 {
00128 G4double p2=p*p;
00129 El=1./(.00012+p2*(.051+.1*p2));
00130 To=El;
00131 }
00132 else if(p>pma)
00133 {
00134 G4double lp=std::log(p)-lmi;
00135 G4double lp2=lp*lp;
00136 El=pbe*lp2+6.72;
00137 To=pbt*lp2+38.2;
00138 }
00139 else
00140 {
00141 G4double p2=p*p;
00142 G4double LE=1./(.00012+p2*(.051+.1*p2));
00143 G4double lp=std::log(p)-lmi;
00144 G4double lp2=lp*lp;
00145 G4double rp2=1./p2;
00146 El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p);
00147 To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2);
00148 }
00149 }
00150 else if(I==2)
00151 {
00152 G4double lp=std::log(p);
00153 if(p<pmi)
00154 {
00155 G4double lr=lp+1.27;
00156 El=1.53/(lr*lr+.0676);
00157 To=El*3;
00158 }
00159 else if(p>pma)
00160 {
00161 G4double ld=lp-lmi;
00162 G4double ld2=ld*ld;
00163 G4double sp=std::sqrt(p);
00164 El=pbe*ld2+2.4+7./sp;
00165 To=pbt*ld2+22.3+12./sp;
00166 }
00167 else
00168 {
00169 G4double lr=lp+1.27;
00170 G4double LE=1.53/(lr*lr+.0676);
00171 G4double ld=lp-lmi;
00172 G4double ld2=ld*ld;
00173 G4double p2=p*p;
00174 G4double p4=p2*p2;
00175 G4double sp=std::sqrt(p);
00176 G4double lm=lp+.36;
00177 G4double md=lm*lm+.04;
00178 G4double lh=lp-.017;
00179 G4double hd=lh*lh+.0025;
00180 El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd;
00181 To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+.4/p4)+1./md+.06/hd;
00182 }
00183 }
00184 else if(I==3)
00185 {
00186 G4double lp=std::log(p);
00187 if(p<pmi)
00188 {
00189 G4double lr=lp+1.27;
00190 G4double lr2=lr*lr;
00191 El=13./(lr2+lr2*lr2+.0676);
00192 To=El;
00193 }
00194 else if(p>pma)
00195 {
00196 G4double ld=lp-lmi;
00197 G4double ld2=ld*ld;
00198 G4double sp=std::sqrt(p);
00199 El=pbe*ld2+2.4+6./sp;
00200 To=pbt*ld2+22.3+5./sp;
00201 }
00202 else
00203 {
00204 G4double lr=lp+1.27;
00205 G4double lr2=lr*lr;
00206 G4double LE=13./(lr2+lr2*lr2+.0676);
00207 G4double ld=lp-lmi;
00208 G4double ld2=ld*ld;
00209 G4double p2=p*p;
00210 G4double p4=p2*p2;
00211 G4double sp=std::sqrt(p);
00212 G4double lm=lp-.32;
00213 G4double md=lm*lm+.0576;
00214 El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md;
00215 To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md;
00216 }
00217 }
00218 else if(I==4)
00219 {
00220
00221 if(p<pmi)
00222 {
00223 G4double psp=p*std::sqrt(p);
00224 El=5.2/psp;
00225 To=14./psp;
00226 }
00227 else if(p>pma)
00228 {
00229 G4double ld=std::log(p)-lmi;
00230 G4double ld2=ld*ld;
00231 El=pbe*ld2+2.23;
00232 To=pbt*ld2+19.5;
00233 }
00234 else
00235 {
00236 G4double ld=std::log(p)-lmi;
00237 G4double ld2=ld*ld;
00238 G4double sp=std::sqrt(p);
00239 G4double psp=p*sp;
00240 G4double p2=p*p;
00241 G4double p4=p2*p2;
00242 G4double lm=p-.39;
00243 G4double md=lm*lm+.000156;
00244 G4double lh=p-1.;
00245 G4double hd=lh*lh+.0156;
00246 El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd;
00247 To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd;
00248 }
00249 }
00250 else if(I==5)
00251 {
00252 if(p<pmi)
00253 {
00254 G4double lr=p-.38;
00255 G4double lm=p-1.;
00256 G4double md=lm*lm+.372;
00257 El=.7/(lr*lr+.0676)+2./md;
00258 To=El+.6/md;
00259 }
00260 else if(p>pma)
00261 {
00262 G4double ld=std::log(p)-lmi;
00263 G4double ld2=ld*ld;
00264 El=pbe*ld2+2.23;
00265 To=pbt*ld2+19.5;
00266 }
00267 else
00268 {
00269 G4double ld=std::log(p)-lmi;
00270 G4double ld2=ld*ld;
00271 G4double lr=p-.38;
00272 G4double LE=.7/(lr*lr+.0676);
00273 G4double sp=std::sqrt(p);
00274 G4double p2=p*p;
00275 G4double p4=p2*p2;
00276 G4double lm=p-1.;
00277 G4double md=lm*lm+.372;
00278 El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md;
00279 To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md;
00280 }
00281 }
00282 else if(I==6)
00283 {
00284 if(p<pmi)
00285 {
00286 G4double p2=p*p;
00287 El=1./(.002+p2*(.12+p2));
00288 To=El;
00289 }
00290 else if(p>pma)
00291 {
00292 G4double lp=std::log(p)-lmi;
00293 G4double lp2=lp*lp;
00294 G4double sp=std::sqrt(p);
00295 El=(pbe*lp2+6.72)/(1.+2./sp);
00296 To=(pbt*lp2+38.2+900./sp)/(1.+27./sp);
00297 }
00298 else
00299 {
00300 G4double p2=p*p;
00301 G4double LE=1./(.002+p2*(.12+p2));
00302 G4double lp=std::log(p)-lmi;
00303 G4double lp2=lp*lp;
00304 G4double p4=p2*p2;
00305 G4double sp=std::sqrt(p);
00306 El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4);
00307 To=LE+(pbt*lp2+38.2+900./sp)/(1.+27./sp+3./p4);
00308 }
00309 }
00310 else if(I==7)
00311 {
00312 if(p>pma)
00313 {
00314 G4double lp=std::log(p)-lmi;
00315 G4double lp2=lp*lp;
00316 El=pbe*lp2+6.72;
00317 To=pbt*lp2+38.2;
00318 }
00319 else
00320 {
00321 G4double ye=std::pow(p,1.25);
00322 G4double yt=std::pow(p,.35);
00323 G4double lp=std::log(p)-lmi;
00324 G4double lp2=lp*lp;
00325 El=80./(ye+1.)+pbe*lp2+6.72;
00326 To=(80./yt+.3)/yt+pbt*lp2+38.2;
00327 }
00328 }
00329 else
00330 {
00331 G4cout<<"*Error*G4QFreeScattering::CalcElTot:ind="<<I<<" is not defined (0-7)"<<G4endl;
00332 G4Exception("G4QFreeScattering::CalcElTot:","23",FatalException,"CHIPScrash");
00333 }
00334 if(El>To) El=To;
00335 return std::make_pair(El,To);
00336 }
00337
00338
00339
00340 std::pair<G4double,G4double> G4QFreeScattering::GetElTotXS(G4double p, G4int PDG, G4bool F)
00341 {
00342 G4int ind=0;
00343 G4bool kfl=true;
00344 G4bool kf=false;
00345 if(PDG==90001000) PDG=2212;
00346 if(PDG==90000001) PDG=2112;
00347 if(PDG==91000000) PDG=3122;
00348 if(PDG==130 || PDG==310)
00349 {
00350 kf=true;
00351 if(G4UniformRand()>.5) kfl=false;
00352 }
00353 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0;
00354 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1;
00355 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2;
00356 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3;
00357 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4;
00358 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5;
00359 else if ( PDG > 3000 && PDG < 3335) ind=6;
00360 else if ( PDG > -3335 && PDG < -2000) ind=7;
00361 else {
00362
00363
00364 ind = 0;
00365 G4cout<<"*Error*G4QFreeScattering::CalcElTotXS: PDG="<<PDG
00366 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
00367 G4Exception("G4QFreeScattering::CalcElTotXS:","22",JustWarning,"CHIPS_crash");
00368 }
00369 return CalcElTot(p,ind);
00370 }
00371
00372
00373 std::pair<G4double,G4double> G4QFreeScattering::FetchElTot(G4double p, G4int PDG, G4bool F)
00374 {
00375 static const G4int nlp=300;
00376 static const G4int mlp=nlp+1;
00377 static const G4double lpi=-5.;
00378 static const G4double lpa=10.;
00379 static const G4double mi=std::exp(lpi);
00380 static const G4double ma=std::exp(lpa);
00381 static const G4double dl=(lpa-lpi)/nlp;
00382 static const G4double edl=std::exp(dl);
00383
00384 static G4double lastP=0.;
00385 static G4int lastH=0;
00386 static G4bool lastF=true;
00387 static std::pair<G4double,G4double> lastR(0.,0.);
00388
00389 static std::vector<G4int> vI;
00390 static std::vector<G4double> vM;
00391 static std::vector<G4int> vK;
00392
00393 static G4int lastI=0;
00394 static G4double lastM=0.;
00395 static G4int lastK=0;
00396 static std::pair<G4double,G4double>* lastX=0;
00397
00398 G4int nDB=vI.size();
00399 if(PDG==90001000) PDG=2212;
00400 if(PDG==90000001) PDG=2112;
00401 if(PDG==91000000) PDG=3122;
00402 #ifdef pdebug
00403 G4cout<<"G4QFreeScatter::FetchElTot:p="<<p<<",PDG="<<PDG<<",F="<<F<<",nDB="<<nDB<<G4endl;
00404 #endif
00405 if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP) return lastR;
00406
00407 lastH=PDG;
00408 lastF=F;
00409 G4int ind=-1;
00410
00411
00412 G4bool kfl=true;
00413 G4bool kf=false;
00414 if(PDG==130 || PDG==310)
00415 {
00416 kf=true;
00417 if(G4UniformRand()>.5) kfl=false;
00418 }
00419 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0;
00420 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1;
00421 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2;
00422 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3;
00423 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4;
00424 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5;
00425 else if ( PDG > 3000 && PDG < 3335) ind=6;
00426 else if ( PDG > -3335 && PDG < -2000) ind=7;
00427 else
00428 {
00429
00430
00431 ind = 0;
00432 G4cout<<"*Error*G4QFreeScattering::FetchElTot: PDG="<<PDG
00433 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
00434 G4Exception("G4QFreeScattering::FetchElTot:","22",JustWarning,"CHIPS problem");
00435 }
00436 if(nDB && lastI==ind && p>0. && p==lastP) return lastR;
00437
00438 if(p<=mi || p>=ma) return CalcElTot(p,ind);
00439 G4bool found=false;
00440 G4int i=-1;
00441 if(nDB) for (i=0; i<nDB; ++i) if(ind==vI[i])
00442 {
00443 found=true;
00444 break;
00445 }
00446 G4double lp=std::log(p);
00447 #ifdef pdebug
00448 G4cout<<"G4QFreeScat::FetchElTot: I="<<ind<<", i="<<i<<",fd="<<found<<",lp="<<lp<<G4endl;
00449 #endif
00450 if(!nDB || !found)
00451 {
00452 #ifdef pdebug
00453 G4cout<<"G4QFreeScattering::FetchElTot: NewX, ind="<<ind<<", nDB="<<nDB<<G4endl;
00454 #endif
00455 lastX = new std::pair<G4double,G4double>[mlp];
00456 lastI = ind;
00457 lastK = static_cast<int>((lp-lpi)/dl)+1;
00458 if(lastK>nlp)
00459 {
00460 lastK=nlp;
00461 lastM=lpa-lpi;
00462 }
00463 else lastM = lastK*dl;
00464 G4double pv=mi;
00465 for(G4int j=0; j<=lastK; ++j)
00466 {
00467 lastX[j]=CalcElTot(pv,ind);
00468 #ifdef pdebug
00469 G4cout<<"G4QFreeScat::FetchElTot:I,j="<<j<<",pv="<<pv<<",E="<<lastX[j].first<<",T="
00470 <<lastX[j].second<<G4endl;
00471 #endif
00472 if(j!=lastK) pv*=edl;
00473 }
00474 i++;
00475 vI.push_back(lastI);
00476 vM.push_back(lastM);
00477 vK.push_back(lastK);
00478 vX.push_back(lastX);
00479 }
00480 else
00481 {
00482 lastI=vI[i];
00483 lastM=vM[i];
00484 lastK=vK[i];
00485 lastX=vX[i];
00486 G4int nextK=lastK+1;
00487 G4double lpM=lastM+lpi;
00488 #ifdef pdebug
00489 G4cout<<"G4QFreeScatt::FetchElTo:M="<<lpM<<",l="<<lp<<",K="<<lastK<<",n="<<nlp<<G4endl;
00490 #endif
00491 if(lp>lpM && lastK<nlp)
00492 {
00493 lastK = static_cast<int>((lp-lpi)/dl)+1;
00494 #ifdef pdebug
00495 G4cout<<"G4QFreeScat::FetET: K="<<lastK<<",lp="<<lp<<",li="<<lpi<<",dl="<<dl<<G4endl;
00496 #endif
00497 if(lastK>nlp)
00498 {
00499 lastK=nlp;
00500 lastM=lpa-lpi;
00501 }
00502 else lastM = lastK*dl;
00503 G4double pv=std::exp(lpM);
00504 for(G4int j=nextK; j<=lastK; ++j)
00505 {
00506 pv*=edl;
00507 lastX[j]=CalcElTot(pv,ind);
00508 #ifdef pdebug
00509 G4cout<<"G4QFreeScat::FetchElTot: U:j="<<j<<",p="<<pv<<",E="<<lastX[j].first<<",T="
00510 <<lastX[j].second<<G4endl;
00511 #endif
00512 }
00513 }
00514 if(lastK >= nextK)
00515 {
00516 vM[i]=lastM;
00517 vK[i]=lastK;
00518 }
00519 }
00520
00521 G4double dlp=lp-lpi;
00522 G4int n=static_cast<int>(dlp/dl);
00523 G4double d=dlp-n*dl;
00524 G4double e=lastX[n].first;
00525 lastR.first=e+d*(lastX[n+1].first-e)/dl;
00526 if(lastR.first<0.) lastR.first = 0.;
00527 G4double t=lastX[n].second;
00528 lastR.second=t+d*(lastX[n+1].second-t)/dl;
00529 if(lastR.second<0.) lastR.second= 0.;
00530 #ifdef pdebug
00531 G4cout<<"=O=>G4QFreeScat::FetchElTot:1st="<<lastR.first<<", 2nd="<<lastR.second<<G4endl;
00532 #endif
00533 if(lastR.first>lastR.second) lastR.first = lastR.second;
00534 return lastR;
00535 }
00536
00537
00538 std::pair<G4double,G4double> G4QFreeScattering::GetElTotMean(G4double pIU, G4int hPDG,
00539 G4int Z, G4int N)
00540 {
00541 G4double pGeV=pIU/gigaelectronvolt;
00542 if(hPDG==90001000) hPDG=2212;
00543 if(hPDG==90000001) hPDG=2112;
00544 if(hPDG==91000000) hPDG=3122;
00545 #ifdef pdebug
00546 G4cout<<"->G4QFreeSc::GetElTotMean: P="<<pIU<<",pPDG="<<hPDG<<",Z="<<Z<<",N="<<N<<G4endl;
00547 #endif
00548 if(Z<1 && N<1)
00549 {
00550 G4cout<<"-Warning-G4QFreeScat::GetElTotMean: Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
00551 return std::make_pair(0.,0.);
00552 }
00553 std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
00554 std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
00555 #ifdef pdebug
00556 G4cout<<"-OUT->G4QFreeScat::GetElTotMean: hp("<<hp.first<<","<<hp.second<<"), hn("
00557 <<hn.first<<","<<hn.second<<")"<<G4endl;
00558 #endif
00559 G4double A=(Z+N)/millibarn;
00560 return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
00561 }
00562
00563
00564
00565 std::pair<G4LorentzVector,G4LorentzVector> G4QFreeScattering::Scatter(G4int NPDG,
00566 G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
00567 {
00568 static const G4double mNeut= G4QPDGCode(2112).GetMass();
00569 static const G4double mProt= G4QPDGCode(2212).GetMass();
00570 p4M/=megaelectronvolt;
00571 N4M/=megaelectronvolt;
00572 G4LorentzVector tot4M=N4M+p4M;
00573 if(pPDG==90001000) pPDG=2212;
00574 if(pPDG==90000001) pPDG=2112;
00575 if(pPDG==91000000) pPDG=3122;
00576 #ifdef ppdebug
00577 G4cout<<"->G4QFR::Scat:p4M="<<p4M<<",N4M="<<N4M<<",t4M="<<tot4M<<",NPDG="<<NPDG<<G4endl;
00578 #endif
00579 G4double mT=mNeut;
00580 #ifdef ppdebug
00581 G4int Z=0;
00582 G4int N=1;
00583 #endif
00584 if(NPDG==2212 || NPDG==90001000)
00585 {
00586 mT=mProt;
00587 #ifdef ppdebug
00588 Z=1;
00589 N=0;
00590 #endif
00591 }
00592 else if(NPDG!=2112 && NPDG!=90000001)
00593 {
00594 G4cout<<"Error:G4QFreeScattering::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
00595 G4Exception("G4QFreeScattering::Scatter:","21",FatalException,"CHIPScomplain");
00596
00597 }
00598 G4double mT2=mT*mT;
00599 G4double mP2=p4M.m2();
00600 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
00601 #ifdef pdebug
00602 G4cout<<"G4QFS::Scat:qM="<<mT<<",qM2="<<mT2<<",pM2="<<mP2<<",totM2="<<tot4M.m2()<<G4endl;
00603 #endif
00604 G4double E2=E*E;
00605 if( E < 0. || E2 < mP2)
00606 {
00607 #ifdef ppdebug
00608 G4cout<<"-Warning-G4QFS::Scat:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl;
00609 #endif
00610 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00611 }
00612 G4double pP2=E2-mP2;
00613 G4double P=std::sqrt(pP2);
00614 #ifdef ppdebug
00615 G4cout<<"G4QFreeS::Scatter: Before XS, P="<<P<<",Z="<<Z<<",N="<<N<<",PDG="<<pPDG<<G4endl;
00616 #endif
00617
00618 if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QFreeScat::Scatter: pPDG="<<pPDG<<G4endl;
00619
00620 G4double dmT=mT+mT;
00621 G4double s_value=dmT*E2+mP2+mT2;
00622 G4double maxt=dmT*dmT*pP2/s_value;
00623 G4double mint=0.;
00624 if(P < 14.) mint=maxt*G4UniformRand();
00625 else
00626 {
00627 G4double p4=pP2*pP2/1000000.;
00628 G4double theB1=(7.2+4.32/p4/(p4+12./P))/(1.+2.5/p4);
00629 mint=-std::log(G4UniformRand())/theB1;
00630 }
00631 #ifdef ppdebug
00632 G4cout<<"G4QFS::Scat:PDG="<<pPDG<<",P="<<P<<",-t="<<mint<<"<"<<maxt<<", Z="<<Z<<",N="<<N
00633 <<G4endl;
00634 #endif
00635 #ifdef nandebug
00636 if(mint>-.0000001);
00637 else G4cout<<"*Warning*G4QFreeScattering::Scatter: -t="<<mint<<G4endl;
00638 #endif
00639 G4double cost=1.-(mint+mint)/maxt;
00640 #ifdef ppdebug
00641 G4cout<<"G4QFS::Scat:-t="<<mint<<"<"<<maxt<<", cost="<<cost<<", Z="<<Z<<",N="<<N<<G4endl;
00642 #endif
00643 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
00644 {
00645 if (cost > 1.) cost = 1.;
00646 else if(cost <-1.) cost =-1.;
00647 else
00648 {
00649 G4cerr<<"G4QFreeScatter::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<maxt<<G4endl;
00650 return std::make_pair(G4LorentzVector(0.,0.,0.,0.), p4M);
00651 }
00652 }
00653 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT);
00654 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
00655 if(!G4QHadron(tot4M).RelDecayIn2(p4M, reco4M, dir4M, cost, cost))
00656 {
00657 G4cerr<<"G4QFS::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl;
00658
00659 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00660 }
00661 #ifdef ppdebug
00662 G4cout<<"G4QFS::Scat:p4M="<<p4M<<"+r4M="<<reco4M<<",dr="<<dir4M<<",t4M="<<tot4M<<G4endl;
00663 #endif
00664 return std::make_pair( reco4M*megaelectronvolt, p4M*megaelectronvolt );
00665 }
00666
00667 G4QHadronVector* G4QFreeScattering::InElF(G4int NPDG, G4LorentzVector N4M,
00668 G4int pPDG, G4LorentzVector p4M)
00669 {
00670 static const G4double mPi0 = G4QPDGCode(111).GetMass();
00671 static const G4double mPi = G4QPDGCode(211).GetMass();
00672 static const G4double mK = G4QPDGCode(321).GetMass();
00673 static const G4double mK0 = G4QPDGCode(311).GetMass();
00674 static const G4double mNeut= G4QPDGCode(2112).GetMass();
00675 static const G4double mProt= G4QPDGCode(2212).GetMass();
00676 static const G4double mLamb= G4QPDGCode(3122).GetMass();
00677
00678 static const G4double mSigM= G4QPDGCode(3112).GetMass();
00679 static const G4double mSigP= G4QPDGCode(3222).GetMass();
00680 static const G4double mXiM = G4QPDGCode(3312).GetMass();
00681 static const G4double mXiZ = G4QPDGCode(3322).GetMass();
00682
00683 static const G4double third =1./3.;
00684 static const G4double twothd =2./3.;
00685
00686 p4M/=megaelectronvolt;
00687 N4M/=megaelectronvolt;
00688 G4LorentzVector c4M=N4M+p4M;
00689 if(pPDG==90001000) pPDG=2212;
00690 if(pPDG==90000001) pPDG=2112;
00691 if(pPDG==91000000) pPDG=3122;
00692 #ifdef ppdebug
00693 G4cout<<"->G4QFR::InElF: p4M="<<p4M<<",N4M="<<N4M<<", c4M="<<c4M<<",NPDG="<<NPDG<<G4endl;
00694 #endif
00695 G4double mT=mNeut;
00696 G4int Z=0;
00697
00698 if(NPDG==2212 || NPDG==90001000)
00699 {
00700 NPDG=2212;
00701 mT=mProt;
00702 Z=1;
00703
00704 }
00705 else if(NPDG!=2112 && NPDG!=90000001)
00706 {
00707 G4cout<<"Error:G4QFreeScattering::InElF:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
00708 G4Exception("G4QFreeScattering::InElF:","21",FatalException,"CHIPS complain");
00709 }
00710 else NPDG=2112;
00711 G4double mC2=c4M.m2();
00712 G4double mC=0.;
00713 if( mC2 < -0.0000001 )
00714 {
00715 #ifdef ppdebug
00716 G4cout<<"-Warning-G4QFS::InElF: Negative compoundMass ="<<mC2<<", c4M="<<c4M<<G4endl;
00717 #endif
00718 return 0;
00719 }
00720 else if ( mC2 > 0.0000001) mC=std::sqrt(mC2);
00721 #ifdef ppdebug
00722 G4cout<<"G4QFS::InElF: mC ="<<mC<<G4endl;
00723 #endif
00724 G4double mP=0.;
00725 G4double mP2=p4M.m2();
00726 if( mP2 < -0.0000001 )
00727 {
00728 #ifdef ppdebug
00729 G4cout<<"-Warning-G4QFS::InElF: Negative projectileMass ="<<mC2<<", c4M="<<c4M<<G4endl;
00730 #endif
00731 return 0;
00732 }
00733 else if ( mP2 > 0.0000001) mP=std::sqrt(mP2);
00734 #ifdef pdebug
00735 G4cout<<"G4QFS::InElF:mT("<<mT<<")+mP("<<mP<<")+mPi0="<<mP+mT+mPi0<<"<? mC="<<mC<<G4endl;
00736 #endif
00737 if(pPDG > 3334 || pPDG < -321)
00738 {
00739 G4cout<<"-Warning-G4QFreeScat::InElF: pPDG="<<pPDG<<G4endl;
00740 return 0;
00741 }
00742 if(pPDG==130 || pPDG==310)
00743 {
00744 if(G4UniformRand()>.5) pPDG = 311;
00745 else pPDG =-311;
00746 }
00747 G4int mPDG=111;
00748 G4double mM=mPi0;
00749 G4double r=G4UniformRand();
00750 if (pPDG == 2212)
00751 {
00752 if(Z)
00753 {
00754 if(r<twothd)
00755 {
00756 pPDG=2112;
00757 mP=mNeut;
00758 mPDG= 211;
00759 mM=mPi;
00760 }
00761 }
00762 else
00763 {
00764 if(r<third)
00765 {
00766 pPDG=2112;
00767 mP=mNeut;
00768 mPDG= 211;
00769 mM=mPi;
00770 }
00771 else if(r<twothd)
00772 {
00773 mPDG=-211;
00774 mM=mPi;
00775 NPDG=2212;
00776 mT=mProt;
00777 }
00778 }
00779 }
00780 else if (pPDG == 2112)
00781 {
00782 if(Z)
00783 {
00784 if(r<third)
00785 {
00786 mPDG= 211;
00787 mM=mPi;
00788 NPDG=2112;
00789 mT=mNeut;
00790 }
00791 else if(r<twothd)
00792 {
00793 pPDG=2212;
00794 mP=mProt;
00795 mPDG=-211;
00796 mM=mPi;
00797 }
00798 }
00799 else
00800 {
00801 if(r<twothd)
00802 {
00803 pPDG=2212;
00804 mP=mProt;
00805 mPDG= -211;
00806 mM=mPi;
00807 }
00808 }
00809 }
00810 else if (pPDG == 211)
00811 {
00812 if(Z)
00813 {
00814 if(r<0.5)
00815 {
00816 mPDG= 211;
00817 mM=mPi;
00818 NPDG=2112;
00819 mT=mNeut;
00820 }
00821 }
00822 else
00823 {
00824 if(r<third)
00825 {
00826 pPDG= 111;
00827 mP=mPi0;
00828 NPDG=2212;
00829 mT=mProt;
00830 }
00831 else if(r<twothd)
00832 {
00833 mPDG=-211;
00834 mM=mPi;
00835 NPDG=2212;
00836 mT=mProt;
00837 }
00838 }
00839 }
00840 else if (pPDG == -211)
00841 {
00842 if(Z)
00843 {
00844 if(r<third)
00845 {
00846 pPDG= 111;
00847 mP=mPi0;
00848 NPDG=2112;
00849 mT=mNeut;
00850 }
00851 else if(r<twothd)
00852 {
00853 mPDG= 211;
00854 mM=mPi;
00855 NPDG=2112;
00856 mT=mNeut;
00857 }
00858 }
00859 else
00860 {
00861 if(r<0.5)
00862 {
00863 mPDG=-211;
00864 mM=mPi;
00865 NPDG=2212;
00866 mT=mProt;
00867 }
00868 }
00869 }
00870 else if (pPDG == -321)
00871 {
00872 if(Z)
00873 {
00874 if(r<0.25)
00875 {
00876 pPDG=-311;
00877 mP=mK0;
00878 NPDG=2112;
00879 mT=mNeut;
00880 }
00881 else if(r<0.5)
00882 {
00883 mPDG= 211;
00884 mM=mPi;
00885 NPDG=2112;
00886 mT=mNeut;
00887 }
00888 else if(r<0.75)
00889 {
00890 pPDG=-311;
00891 mP=mK0;
00892 mPDG=-211;
00893 mM=mPi;
00894 }
00895 }
00896 else
00897 {
00898 if(r<third)
00899 {
00900 mPDG=-211;
00901 mM=mPi;
00902 NPDG=2212;
00903 mT=mProt;
00904 }
00905 else if(r<twothd)
00906 {
00907 pPDG=-311;
00908 mP=mK0;
00909 mPDG=-211;
00910 mM=mPi;
00911 }
00912 }
00913 }
00914 else if (pPDG == -311)
00915 {
00916 if(Z)
00917 {
00918 if(r<third)
00919 {
00920 pPDG=-321;
00921 mP=mK;
00922 NPDG=2212;
00923 mT=mProt;
00924 }
00925 else if(r<twothd)
00926 {
00927 mPDG= 211;
00928 mM=mPi;
00929 NPDG=2112;
00930 mT=mNeut;
00931 }
00932 }
00933 else
00934 {
00935 if(r<0.25)
00936 {
00937 pPDG=-321;
00938 mP=mK;
00939 mPDG= 211;
00940 mM=mPi;
00941 }
00942 else if(r<0.5)
00943 {
00944 pPDG=-321;
00945 mP=mK;
00946 NPDG=2212;
00947 mT=mProt;
00948 }
00949 else if(r<0.75)
00950 {
00951 mPDG=-211;
00952 mM=mPi;
00953 NPDG=2212;
00954 mT=mProt;
00955 }
00956 }
00957 }
00958 else if (pPDG == 321)
00959 {
00960 if(Z)
00961 {
00962 if(r<third)
00963 {
00964 mPDG= 211;
00965 mM=mPi;
00966 NPDG=2112;
00967 mT=mNeut;
00968 }
00969 else if(r<twothd)
00970 {
00971 pPDG= 311;
00972 mP=mK0;
00973 mPDG= 211;
00974 mM=mPi;
00975 }
00976 }
00977 else
00978 {
00979 if(r<0.25)
00980 {
00981 pPDG= 311;
00982 mP=mK0;
00983 mPDG= 211;
00984 mM=mPi;
00985 }
00986 else if(r<0.5)
00987 {
00988 mPDG=-211;
00989 mM=mPi;
00990 NPDG=2212;
00991 mT=mProt;
00992 }
00993 else if(r<0.75)
00994 {
00995 pPDG= 311;
00996 mP=mK0;
00997 NPDG=2212;
00998 mT=mProt;
00999 }
01000 }
01001 }
01002 else if (pPDG == 311)
01003 {
01004 if(Z)
01005 {
01006 if(r<0.25)
01007 {
01008 pPDG= 321;
01009 mP=mK;
01010 mPDG=-211;
01011 mM=mPi;
01012 }
01013 else if(r<0.5)
01014 {
01015 mPDG= 211;
01016 mM=mPi;
01017 NPDG=2112;
01018 mT=mNeut;
01019 }
01020 else if(r<0.75)
01021 {
01022 pPDG= 321;
01023 mP=mK;
01024 NPDG=2112;
01025 mT=mNeut;
01026 }
01027 }
01028 else
01029 {
01030 if(r<third)
01031 {
01032 mPDG=-211;
01033 mM=mPi;
01034 NPDG=2212;
01035 mT=mProt;
01036 }
01037 else if(r<twothd)
01038 {
01039 pPDG= 321;
01040 mP=mK;
01041 mPDG=-211;
01042 mM=mPi;
01043 }
01044 }
01045 }
01046 else if (pPDG == 3122 || pPDG== 3212)
01047 {
01048 if(pPDG == 3212)
01049 {
01050 pPDG=3122;
01051 mP=mLamb;
01052 }
01053 if(Z)
01054 {
01055 if(r<0.2)
01056 {
01057 pPDG=3222;
01058 mP=mSigP;
01059 NPDG=2112;
01060 mT=mNeut;
01061 }
01062 else if(r<0.4)
01063 {
01064 pPDG=3222;
01065 mP=mSigP;
01066 mPDG=-211;
01067 mM=mPi;
01068 }
01069 else if(r<0.6)
01070 {
01071 pPDG=3112;
01072 mP=mSigM;
01073 mPDG= 211;
01074 mM=mPi;
01075 }
01076 else if(r<0.8)
01077 {
01078 mPDG= 211;
01079 mM=mPi;
01080 NPDG=2112;
01081 mT=mNeut;
01082 }
01083 }
01084 else
01085 {
01086 if(r<0.2)
01087 {
01088 pPDG=3112;
01089 mP=mSigM;
01090 NPDG=2212;
01091 mT=mProt;
01092 }
01093 else if(r<0.4)
01094 {
01095 pPDG=3222;
01096 mP=mSigP;
01097 mPDG=-211;
01098 mM=mPi;
01099 }
01100 else if(r<0.6)
01101 {
01102 pPDG=3112;
01103 mP=mSigM;
01104 mPDG= 211;
01105 mM=mPi;
01106 }
01107 else if(r<0.8)
01108 {
01109 mPDG=-211;
01110 mM=mPi;
01111 NPDG=2212;
01112 mT=mProt;
01113 }
01114 }
01115 }
01116 else if (pPDG == 3112)
01117 {
01118 if(Z)
01119 {
01120 if(r<0.25)
01121 {
01122 pPDG=3122;
01123 mP=mLamb;
01124 NPDG=2112;
01125 mT=mNeut;
01126 }
01127 else if(r<0.5)
01128 {
01129 pPDG=3122;
01130 mP=mLamb;
01131 mPDG=-211;
01132 mM=mPi;
01133 }
01134 else if(r<0.75)
01135 {
01136 mPDG= 211;
01137 mM=mPi;
01138 NPDG=2112;
01139 mT=mNeut;
01140 }
01141 }
01142 else
01143 {
01144 if(r<third)
01145 {
01146 pPDG=3122;
01147 mP=mLamb;
01148 mPDG=-211;
01149 mM=mPi;
01150 }
01151 else if(r<twothd)
01152 {
01153 mPDG=-211;
01154 mM=mPi;
01155 NPDG=2212;
01156 mT=mProt;
01157 }
01158 }
01159 }
01160 else if (pPDG == 3222)
01161 {
01162 if(Z)
01163 {
01164 if(r<third)
01165 {
01166 pPDG=3122;
01167 mP=mLamb;
01168 mPDG= 211;
01169 mM=mPi;
01170 }
01171 else if(r<twothd)
01172 {
01173 mPDG= 211;
01174 mM=mPi;
01175 NPDG=2112;
01176 mT=mNeut;
01177 }
01178 }
01179 else
01180 {
01181 if(r<0.25)
01182 {
01183 pPDG=3122;
01184 mP=mLamb;
01185 NPDG=2212;
01186 mT=mProt;
01187 }
01188 else if(r<0.5)
01189 {
01190 pPDG=3122;
01191 mP=mLamb;
01192 mPDG= 211;
01193 mM=mPi;
01194 }
01195 else if(r<0.75)
01196 {
01197 mPDG=-211;
01198 mM=mPi;
01199 NPDG=2212;
01200 mT=mProt;
01201 }
01202 }
01203 }
01204 else if (pPDG == 3312)
01205 {
01206 if(Z)
01207 {
01208 if(r<0.25)
01209 {
01210 pPDG=3322;
01211 mP=mXiZ;
01212 NPDG=2112;
01213 mT=mNeut;
01214 }
01215 else if(r<0.5)
01216 {
01217 pPDG=3322;
01218 mP=mXiZ;
01219 mPDG=-211;
01220 mM=mPi;
01221 }
01222 else if(r<0.75)
01223 {
01224 mPDG= 211;
01225 mM=mPi;
01226 NPDG=2112;
01227 mT=mNeut;
01228 }
01229 }
01230 else
01231 {
01232 if(r<third)
01233 {
01234 pPDG=3322;
01235 mP=mXiZ;
01236 mPDG=-211;
01237 mM=mPi;
01238 }
01239 else if(r<twothd)
01240 {
01241 mPDG=-211;
01242 mM=mPi;
01243 NPDG=2212;
01244 mT=mProt;
01245 }
01246 }
01247 }
01248 else if (pPDG == 3322)
01249 {
01250 if(Z)
01251 {
01252 if(r<third)
01253 {
01254 pPDG=3312;
01255 mP=mXiM;
01256 mPDG= 211;
01257 mM=mPi;
01258 }
01259 else if(r<twothd)
01260 {
01261 mPDG= 211;
01262 mM=mPi;
01263 NPDG=2112;
01264 mT=mNeut;
01265 }
01266 }
01267 else
01268 {
01269 if(r<0.25)
01270 {
01271 pPDG=3312;
01272 mP=mXiM;
01273 NPDG=2212;
01274 mT=mProt;
01275 }
01276 else if(r<0.5)
01277 {
01278 pPDG=3312;
01279 mP=mXiM;
01280 mPDG= 211;
01281 mM=mPi;
01282 }
01283 else if(r<0.75)
01284 {
01285 mPDG=-211;
01286 mM=mPi;
01287 NPDG=2212;
01288 mT=mProt;
01289 }
01290 }
01291 }
01292 else if (pPDG == 3334)
01293 {
01294 if(Z)
01295 {
01296 if(r<0.5)
01297 {
01298 mPDG= 211;
01299 mM=mPi;
01300 NPDG=2112;
01301 mT=mNeut;
01302 }
01303 }
01304 else
01305 {
01306 if(r<0.5)
01307 {
01308 mPDG=-211;
01309 mM=mPi;
01310 NPDG=2212;
01311 mT=mProt;
01312 }
01313 }
01314 }
01315 else
01316 {
01317 G4cout<<"*Error*G4QFreeScattering::InElF: PDG="<<pPDG
01318 <<", while it is defined only for p,n,hyperons(not Omega),pi,K/antiK"<<G4endl;
01319 G4Exception("G4QFreeScattering::InElF:","22",FatalException,"CHIPS_crash");
01320 }
01321 if (mC-mP-mM-mT <-0.000001) return 0;
01322 G4QHadronVector* TripQH = new G4QHadronVector;
01323 G4LorentzVector m4M(0.,0.,0.,mM);
01324 if(mC-mP-mM-mT < 0.000001)
01325 {
01326 p4M=(mP/mC)*c4M;
01327 m4M=(mM/mC)*c4M;
01328 N4M=(mT/mC)*c4M;
01329 }
01330 else
01331 {
01332 p4M=G4LorentzVector(0.,0.,0.,mP);
01333 N4M=G4LorentzVector(0.,0.,0.,mT);
01334 if(!G4QHadron(c4M).DecayIn3(p4M,m4M,N4M))
01335 {
01336 G4ExceptionDescription ed;
01337 ed << "DecayIn3, TotM=" << mC << G4endl;
01338 G4Exception("G4QFreeScattering::InElF()","HAD_CHPS_0027", FatalException, ed);
01339 }
01340 G4QHadron* h1 = new G4QHadron(pPDG, p4M);
01341 TripQH->push_back(h1);
01342 #ifdef debug
01343 G4cout << "G4QFreeScat::InElF: H1=" << pPDG << p4M << G4endl;
01344 #endif
01345 G4QHadron* h2 = new G4QHadron(mPDG, m4M);
01346 TripQH->push_back(h2);
01347 #ifdef debug
01348 G4cout << "G4QFreeScat::InElF: H2=" << mPDG << m4M << G4endl;
01349 #endif
01350 G4QHadron* h3 = new G4QHadron(NPDG, N4M);
01351 TripQH->push_back(h3);
01352 #ifdef debug
01353 G4cout << "G4QFreeScat::InElF: H3=" << NPDG << N4M << G4endl;
01354 #endif
01355 }
01356 return TripQH;
01357 }