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 #include "G4QuasiFreeRatios.hh"
00045 #include "G4SystemOfUnits.hh"
00046
00047
00048 std::vector<G4double*> G4QuasiFreeRatios::vT;
00049 std::vector<G4double*> G4QuasiFreeRatios::vL;
00050
00051 G4QuasiFreeRatios::G4QuasiFreeRatios()
00052 {
00053 #ifdef pdebug
00054 G4cout<<"***^^^*** G4QuasiFreeRatios singletone is created ***^^^***"<<G4endl;
00055 #endif
00056 QFreeScat=G4QFreeScattering::GetPointer();
00057 }
00058
00059 G4QuasiFreeRatios::~G4QuasiFreeRatios()
00060 {
00061 std::vector<G4double*>::iterator pos;
00062 for(pos=vT.begin(); pos<vT.end(); ++pos) delete [] *pos;
00063 vT.clear();
00064 for(pos=vL.begin(); pos<vL.end(); ++pos) delete [] *pos;
00065 vL.clear();
00066 }
00067
00068
00069 G4QuasiFreeRatios* G4QuasiFreeRatios::GetPointer()
00070 {
00071 static G4QuasiFreeRatios theRatios;
00072 return &theRatios;
00073 }
00074
00075
00076 std::pair<G4double,G4double> G4QuasiFreeRatios::GetRatios(G4double pIU, G4int pPDG,
00077 G4int tgZ, G4int tgN)
00078 {
00079 #ifdef pdebug
00080 G4cout<<">>>IN>>>G4QFRat::GetQF:P="<<pIU<<",pPDG="<<pPDG<<",Z="<<tgZ<<",N="<<tgN<<G4endl;
00081 #endif
00082 G4double R=0.;
00083 G4double QF2In=1.;
00084 G4int tgA=tgZ+tgN;
00085 if(tgA<2) return std::make_pair(QF2In,R);
00086 std::pair<G4double,G4double> ElTot=GetElTot(pIU, pPDG, tgZ, tgN);
00087
00088 if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.;
00089 else if(ElTot.second>0.)
00090 {
00091 R=ElTot.first/ElTot.second;
00092 QF2In=GetQF2IN_Ratio(ElTot.second/millibarn, tgZ+tgN);
00093 }
00094 #ifdef pdebug
00095 G4cout<<">>>OUT>>>G4QuasiFreeRatio::GetQF2IN_Ratio: QF2In="<<QF2In<<", R="<<R<<G4endl;
00096 #endif
00097 return std::make_pair(QF2In,R);
00098 }
00099
00100
00101 G4double G4QuasiFreeRatios::GetQF2IN_Ratio(G4double s_value, G4int A)
00102 {
00103 static const G4int nps=150;
00104 static const G4int mps=nps+1;
00105 static const G4double sma=150.;
00106 static const G4double ds=sma/nps;
00107 static const G4int nls=100;
00108 static const G4int mls=nls+1;
00109 static const G4double lsi=5.;
00110 static const G4double lsa=9.;
00111 static const G4double mi=std::exp(lsi);
00112 static const G4double max_s=std::exp(lsa);
00113 static const G4double dl=(lsa-lsi)/nls;
00114 static const G4double edl=std::exp(dl);
00115 static const G4double toler=.01;
00116 static G4double lastS=0.;
00117 static G4double lastR=0.;
00118
00119 static std::vector<G4int> vA;
00120 static std::vector<G4double> vH;
00121 static std::vector<G4int> vN;
00122 static std::vector<G4double> vM;
00123 static std::vector<G4int> vK;
00124
00125 static G4int lastA=0;
00126 static G4double lastH=0.;
00127 static G4int lastN=0;
00128 static G4double lastM=0.;
00129 static G4int lastK=0;
00130 static G4double* lastT=0;
00131 static G4double* lastL=0;
00132
00133 #ifdef pdebug
00134 G4cout<<"+++G4QuasiFreeRatio::GetQF2IN_Ratio:A="<<A<<", s="<<s_value<<G4endl;
00135 #endif
00136 if(s_value<toler || A<2) return 1.;
00137 if(s_value>max_s) return 0.;
00138 if(A>238)
00139 {
00140 G4cout<<"-Warning-G4QuasiFreeRatio::GetQF2IN_Ratio:A="<<A<<">238, return zero"<<G4endl;
00141 return 0.;
00142 }
00143 G4int nDB=vA.size();
00144
00145 if(nDB && lastA==A && s_value==lastS) return lastR;
00146 G4bool found=false;
00147 G4int i=-1;
00148 if(nDB) for (i=0; i<nDB; i++) if(A==vA[i])
00149 {
00150 found=true;
00151 break;
00152 }
00153 #ifdef pdebug
00154 G4cout<<"+++G4QuasiFreeRatio::GetQF2IN_Ratio: nDB="<<nDB<<", found="<<found<<G4endl;
00155 #endif
00156 if(!nDB || !found)
00157 {
00158 lastA = A;
00159 #ifdef pdebug
00160 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: NewT, A="<<A<<", nDB="<<nDB<<G4endl;
00161 #endif
00162 lastT = new G4double[mps];
00163 lastN = static_cast<int>(s_value/ds)+1;
00164 if(lastN>nps)
00165 {
00166 lastN=nps;
00167 lastH=sma;
00168 }
00169 else lastH = lastN*ds;
00170 G4double sv=0;
00171 lastT[0]=1.;
00172 for(G4int j=1; j<=lastN; j++)
00173 {
00174 sv+=ds;
00175 lastT[j]=CalcQF2IN_Ratio(sv,A);
00176 }
00177 lastL=new G4double[mls];
00178 if(s_value>sma)
00179 {
00180 #ifdef pdebug
00181 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: NewL, A="<<A<<", nDB="<<nDB<<G4endl;
00182 #endif
00183 G4double ls=std::log(s_value);
00184 lastK = static_cast<int>((ls-lsi)/dl)+1;
00185 if(lastK>nls)
00186 {
00187 lastK=nls;
00188 lastM=lsa-lsi;
00189 }
00190 else lastM = lastK*dl;
00191 sv=mi;
00192 for(G4int j=0; j<=lastK; j++)
00193 {
00194 lastL[j]=CalcQF2IN_Ratio(sv,A);
00195 if(j!=lastK) sv*=edl;
00196 }
00197 }
00198 else
00199 {
00200 lastK = 0;
00201 lastM = 0.;
00202 }
00203 i++;
00204 vA.push_back(lastA);
00205 vH.push_back(lastH);
00206 vN.push_back(lastN);
00207 vM.push_back(lastM);
00208 vK.push_back(lastK);
00209 vT.push_back(lastT);
00210 vL.push_back(lastL);
00211 }
00212 else
00213 {
00214 lastA=vA[i];
00215 lastH=vH[i];
00216 lastN=vN[i];
00217 lastM=vM[i];
00218 lastK=vK[i];
00219 lastT=vT[i];
00220 lastL=vL[i];
00221 #ifdef pdebug
00222 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: Found, s="<<s_value<<", lastM="<<lastM<<G4endl;
00223 #endif
00224 if(s_value>lastH)
00225 {
00226 G4int nextN=lastN+1;
00227 #ifdef pdebug
00228 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: lastN="<<lastN<<" ?< nps="<<nps<<G4endl;
00229 #endif
00230 if(lastN < nps)
00231 {
00232 lastN = static_cast<int>(s_value/ds)+1;
00233 G4double sv=lastH;
00234 if(lastN>nps)
00235 {
00236 lastN=nps;
00237 lastH=sma;
00238 }
00239 else lastH = lastN*ds;
00240 for(G4int j=nextN; j<=lastN; j++)
00241 {
00242 sv+=ds;
00243 lastT[j]=CalcQF2IN_Ratio(sv,A);
00244 }
00245 #ifdef pdebug
00246 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: End of LinTab update"<<G4endl;
00247 #endif
00248 }
00249 #ifdef pdebug
00250 G4cout<<"G4QFRatios::GetQF2IN_Ratio: lN="<<lastN<<", nN="<<nextN<<", i="<<i<<G4endl;
00251 #endif
00252 if(lastN >= nextN)
00253 {
00254 vH[i]=lastH;
00255 vN[i]=lastN;
00256 }
00257 G4int nextK=lastK+1;
00258 if(!lastK) nextK=0;
00259 #ifdef pdebug
00260 G4cout<<"G4QFRat::GetQF2IN_Ratio: sma="<<sma<<", lastK="<<lastK<<" < "<<nls<<G4endl;
00261 #endif
00262 if(s_value>sma && lastK<nls)
00263 {
00264 G4double sv=std::exp(lastM+lsi);
00265 G4double ls=std::log(s_value);
00266 lastK = static_cast<int>((ls-lsi)/dl)+1;
00267 if(lastK>nls)
00268 {
00269 lastK=nls;
00270 lastM=lsa-lsi;
00271 }
00272 else lastM = lastK*dl;
00273 #ifdef pdebug
00274 G4cout<<"G4QFRat::GetQF2IN_Ratio: nK="<<nextK<<", lK="<<lastK<<", sv="<<sv<<G4endl;
00275 #endif
00276 for(G4int j=nextK; j<=lastK; j++)
00277 {
00278 sv*=edl;
00279 #ifdef pdebug
00280 G4cout<<"G4QFRat::GetQF2IN_Ratio: j="<<j<<", sv="<<sv<<", A="<<A<<G4endl;
00281 #endif
00282 lastL[j]=CalcQF2IN_Ratio(sv,A);
00283 }
00284 #ifdef pdebug
00285 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: End of LinTab update"<<G4endl;
00286 #endif
00287 }
00288 #ifdef pdebug
00289 G4cout<<"G4QFRatios::GetQF2IN_Ratio: lK="<<lastK<<", nK="<<nextK<<", i="<<i<<G4endl;
00290 #endif
00291 if(lastK >= nextK)
00292 {
00293 vM[i]=lastM;
00294 vK[i]=lastK;
00295 }
00296 }
00297 }
00298 #ifdef pdebug
00299 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: BeforeTab s="<<s_value<<", sma="<<sma<<G4endl;
00300 #endif
00301
00302 if(s_value<sma)
00303 {
00304 G4int n=static_cast<int>(s_value/ds);
00305 G4double d=s_value-n*ds;
00306 G4double v=lastT[n];
00307 lastR=v+d*(lastT[n+1]-v)/ds;
00308 }
00309 else
00310 {
00311 G4double ls=std::log(s_value)-lsi;
00312 G4int n=static_cast<int>(ls/dl);
00313 G4double d=ls-n*dl;
00314 G4double v=lastL[n];
00315 lastR=v+d*(lastL[n+1]-v)/dl;
00316 }
00317 if(lastR<0.) lastR=0.;
00318 if(lastR>1.) lastR=1.;
00319 #ifdef pdebug
00320 G4cout<<"G4QuasiFreeRatios::GetQF2IN_Ratio: BeforeRet lastR="<<lastR<<G4endl;
00321 #endif
00322 return lastR;
00323 }
00324
00325
00326 G4double G4QuasiFreeRatios::CalcQF2IN_Ratio(G4double s_value, G4int A)
00327 {
00328 static const G4double C=1.246;
00329 G4double s2=s_value*s_value;
00330 G4double s4=s2*s2;
00331 G4double ss=std::sqrt(std::sqrt(s_value));
00332 G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2);
00333 G4double E=.2644+.016/(1.+std::exp((29.54-s_value)/2.49));
00334 G4double F=ss*.1526*std::exp(-s2*ss*.0000859);
00335 return C*std::exp(-E*std::pow(G4double(A-1.),F))/std::pow(G4double(A),P);
00336 }
00337
00338
00339 std::pair<G4double,G4double> G4QuasiFreeRatios::GetElTotXS(G4double p, G4int PDG, G4bool F)
00340 {
00341 G4int ind=0;
00342 G4bool kfl=true;
00343 G4bool kf=false;
00344 if(PDG==130||PDG==310)
00345 {
00346 kf=true;
00347 if(G4UniformRand()>.5) kfl=false;
00348 }
00349 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0;
00350 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1;
00351 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2;
00352 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3;
00353 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4;
00354 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5;
00355 else if ( PDG > 3000 && PDG < 3335) ind=6;
00356 else if ( PDG > -3335 && PDG < -2000) ind=7;
00357 else {
00358 G4cout<<"*Error*G4QuasiFreeRatios::CalcElTotXS: PDG="<<PDG
00359 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
00360 G4Exception("G4QuasiFreeRatio::CalcElTotXS:","22",FatalException,"CHIPScrash");
00361 }
00362 return QFreeScat->CalcElTot(p,ind);
00363 }
00364
00365
00366 std::pair<G4double,G4double> G4QuasiFreeRatios::GetElTot(G4double pIU, G4int hPDG,
00367 G4int Z, G4int N)
00368 {
00369 G4double pGeV=pIU/gigaelectronvolt;
00370 #ifdef pdebug
00371 G4cout<<"-->G4QuasiFreeR::GetElTot: P="<<pIU<<",pPDG="<<hPDG<<",Z="<<Z<<",N="<<N<<G4endl;
00372 #endif
00373 if(Z<1 && N<1)
00374 {
00375 G4cout<<"-Warning-G4QuasiFreeRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
00376 return std::make_pair(0.,0.);
00377 }
00378 std::pair<G4double,G4double> hp=QFreeScat->FetchElTot(pGeV, hPDG, true);
00379 std::pair<G4double,G4double> hn=QFreeScat->FetchElTot(pGeV, hPDG, false);
00380 #ifdef pdebug
00381 G4cout<<"-OUT->G4QFRat::GetElTot: hp("<<hp.first<<","<<hp.second<<"), hn("<<hn.first<<","
00382 <<hn.second<<")"<<G4endl;
00383 #endif
00384 G4double A=(Z+N)/millibarn;
00385 return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
00386 }
00387
00388
00389 std::pair<G4double,G4double> G4QuasiFreeRatios::GetChExFactor(G4double pIU, G4int hPDG,
00390 G4int Z, G4int N)
00391 {
00392 G4double pGeV=pIU/gigaelectronvolt;
00393 G4double resP=0.;
00394 G4double resN=0.;
00395 if(Z<1 && N<1)
00396 {
00397 G4cout<<"-Warning-G4QuasiFreeRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
00398 return std::make_pair(resP,resN);
00399 }
00400 G4double A=Z+N;
00401 G4double pf=0.;
00402 G4double nf=0.;
00403 if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312||hPDG==3334) pf=Z/(A+N);
00404 else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z);
00405 else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310)
00406 {
00407 G4double dA=A+A;
00408 pf=Z/(dA+N+N);
00409 nf=N/(dA+Z+Z);
00410 }
00411 G4double mult=1.;
00412 if(pGeV>.5)
00413 {
00414 mult=1./(1.+std::log(pGeV+pGeV))/pGeV;
00415 if(mult>1.) mult=1.;
00416 }
00417 if(pf)
00418 {
00419 std::pair<G4double,G4double> hp=QFreeScat->FetchElTot(pGeV, hPDG, true);
00420 resP=pf*(hp.second/hp.first-1.)*mult;
00421 }
00422 if(nf)
00423 {
00424 std::pair<G4double,G4double> hn=QFreeScat->FetchElTot(pGeV, hPDG, false);
00425 resN=nf*(hn.second/hn.first-1.)*mult;
00426 }
00427 return std::make_pair(resP,resN);
00428 }
00429
00430
00431
00432 std::pair<G4LorentzVector,G4LorentzVector> G4QuasiFreeRatios::Scatter(G4int NPDG,
00433 G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
00434 {
00435 static const G4double mNeut= G4QPDGCode(2112).GetMass();
00436 static const G4double mProt= G4QPDGCode(2212).GetMass();
00437 static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);
00438 static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);
00439 static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);
00440 static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);
00441 static const G4double two3d= 2./3.;
00442 static const G4double two3d2= std::pow(2.,two3d);
00443 static const G4double two3d3= std::pow(3.,two3d);
00444 static const G4double two3d4= std::pow(4.,two3d);
00445 G4LorentzVector pr4M=p4M/megaelectronvolt;
00446 N4M/=megaelectronvolt;
00447 G4LorentzVector tot4M=N4M+p4M;
00448 #ifdef ppdebug
00449 G4cerr<<"->G4QFR::Scat:p4M="<<pr4M<<",N4M="<<N4M<<",t4M="<<tot4M<<",NPDG="<<NPDG<<G4endl;
00450 #endif
00451 G4double mT=mNeut;
00452 G4int Z=0;
00453 G4int N=1;
00454 if(NPDG==2212||NPDG==90001000)
00455 {
00456 mT=mProt;
00457 Z=1;
00458 N=0;
00459 }
00460 else if(NPDG==90001001)
00461 {
00462 mT=mDeut;
00463 Z=1;
00464 N=1;
00465 }
00466 else if(NPDG==90002001)
00467 {
00468 mT=mHel3;
00469 Z=2;
00470 N=1;
00471 }
00472 else if(NPDG==90001002)
00473 {
00474 mT=mTrit;
00475 Z=1;
00476 N=2;
00477 }
00478 else if(NPDG==90002002)
00479 {
00480 mT=mAlph;
00481 Z=2;
00482 N=2;
00483 }
00484 else if(NPDG!=2112&&NPDG!=90000001)
00485 {
00486 G4cout<<"Error:G4QuasiFreeRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
00487 G4Exception("G4QuasiFreeRatios::Scatter:","21",FatalException,"CHIPScomplain");
00488
00489 }
00490 G4double mT2=mT*mT;
00491 G4double mP2=pr4M.m2();
00492 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
00493 #ifdef pdebug
00494 G4cerr<<"G4QFR::Scat:qM="<<mT<<",qM2="<<mT2<<",pM2="<<mP2<<",totM2="<<tot4M.m2()<<G4endl;
00495 #endif
00496 G4double E2=E*E;
00497 if(E<0. || E2<mP2)
00498 {
00499 #ifdef ppdebug
00500 G4cerr<<"-Warning-G4QFR::Scat:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl;
00501 #endif
00502 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00503 }
00504 G4double P=std::sqrt(E2-mP2);
00505 G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
00506 G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer();
00507 #ifdef ppdebug
00508 G4cout<<"G4QFR::Scatter: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl;
00509 #endif
00510
00511 if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QuasiFree::Scatter: pPDG="<<pPDG<<G4endl;
00512 G4int PDG=2212;
00513 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;
00514 if(!Z && N==1)
00515 {
00516 Z=1;
00517 N=0;
00518 if (PDG==2212) PDG=2112;
00519 else if(PDG==2112) PDG=2212;
00520 }
00521 G4double xSec=0.;
00522 if(PDG==2212) xSec=PCSmanager->GetCrossSection(false, P, Z, N, PDG);
00523 else xSec=NCSmanager->GetCrossSection(false, P, Z, N, PDG);
00524 #ifdef ppdebug
00525 G4cout<<"G4QuasiFreeRat::Scatter: pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl;
00526 #endif
00527 #ifdef nandebug
00528 if(xSec>0. || xSec<0. || xSec==0);
00529 else G4cout<<"*Warning*G4QuasiFreeRatios::Scatter: xSec="<<xSec/millibarn<<G4endl;
00530 #endif
00531
00532 if(xSec <= 0.)
00533 {
00534 #ifdef ppdebug
00535 G4cerr<<"-Warning-G4QFR::Scat:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl;
00536 #endif
00537 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00538 }
00539 G4double mint=0.;
00540 if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);
00541 else mint=NCSmanager->GetExchangeT(Z,N,PDG);
00542 if (mT==mDeut) mint/=two3d2;
00543 else if(mT==mTrit || mT==mHel3) mint/=two3d3;
00544 else if(mT==mAlph) mint/=two3d4;
00545 G4double maxt=0.;
00546 if(PDG==2212) maxt=PCSmanager->GetHMaxT();
00547 else maxt=NCSmanager->GetHMaxT();
00548 #ifdef ppdebug
00549 G4cout<<"G4QFR::Scat:PDG="<<PDG<<",P="<<P<<",X="<<xSec<<",-t="<<mint<<"<"<<maxt<<", Z="
00550 <<Z<<",N="<<N<<G4endl;
00551 #endif
00552 #ifdef nandebug
00553 if(mint>-.0000001);
00554 else G4cout<<"*Warning*G4QFR::Scat: -t="<<mint<<G4endl;
00555 #endif
00556 G4double cost=1.-(mint+mint)/maxt;
00557 #ifdef ppdebug
00558 G4cout<<"G4QFR::Scat:-t="<<mint<<"<"<<maxt<<", cost="<<cost<<", Z="<<Z<<",N="<<N<<G4endl;
00559 #endif
00560 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
00561 {
00562 if (cost>1.) cost=1.;
00563 else if(cost<-1.) cost=-1.;
00564 else
00565 {
00566 G4double tm=0.;
00567 if(PDG==2212) tm=PCSmanager->GetHMaxT();
00568 else tm=NCSmanager->GetHMaxT();
00569 G4cerr<<"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<tm<<G4endl;
00570 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00571 }
00572 }
00573 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT);
00574 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
00575 if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost))
00576 {
00577 G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl;
00578
00579 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00580 }
00581 #ifdef ppdebug
00582 G4cout<<"G4QFR::Scat:p4M="<<pr4M<<"+r4M="<<reco4M<<",dr="<<dir4M<<",t4M="<<tot4M<<G4endl;
00583 #endif
00584 return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt);
00585 }
00586
00587
00588
00589
00590 std::pair<G4LorentzVector,G4LorentzVector> G4QuasiFreeRatios::ChExer(G4int NPDG,
00591 G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
00592 {
00593 static const G4double mNeut= G4QPDGCode(2112).GetMass();
00594 static const G4double mProt= G4QPDGCode(2212).GetMass();
00595 G4LorentzVector pr4M=p4M/megaelectronvolt;
00596 N4M/=megaelectronvolt;
00597 G4LorentzVector tot4M=N4M+p4M;
00598 G4int Z=0;
00599 G4int N=1;
00600 G4int sPDG=0;
00601 G4double mS=0.;
00602 G4double mT=mProt;
00603 G4cout<<"-Warning-G4QFR::ChEx: Impossible for Omega-"<<G4endl;
00604 if(NPDG==2212)
00605 {
00606 mT=mNeut;
00607 Z=1;
00608 N=0;
00609 if(pPDG==-211) sPDG=111;
00610 else if(pPDG==-321)
00611 {
00612 sPDG=310;
00613 if(G4UniformRand()>.5) sPDG=130;
00614 }
00615 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321;
00616 else if(pPDG==3112) sPDG=3212;
00617 else if(pPDG==3212) sPDG=3222;
00618 else if(pPDG==3312) sPDG=3322;
00619 }
00620 else if(NPDG==2112)
00621 {
00622 if(pPDG==211) sPDG=111;
00623 else if(pPDG==321)
00624 {
00625 sPDG=310;
00626 if(G4UniformRand()>.5) sPDG=130;
00627 }
00628 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321;
00629 else if(pPDG==3222) sPDG=3212;
00630 else if(pPDG==3212) sPDG=3112;
00631 else if(pPDG==3322) sPDG=3312;
00632 }
00633 else
00634 {
00635 G4cout<<"Error:G4QuasiFreeRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
00636 G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain");
00637
00638 }
00639 if(sPDG) mS=mNeut;
00640 else
00641 {
00642 G4cout<<"Error:G4QuasiFreeRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl;
00643 G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain");
00644
00645 }
00646 G4double mT2=mT*mT;
00647 G4double mS2=mS*mS;
00648 G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);
00649 G4double E2=E*E;
00650 if(E<0. || E2<mS2)
00651 {
00652 #ifdef pdebug
00653 G4cerr<<"-Warning-G4QFR::ChEx:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mS2<<G4endl;
00654 #endif
00655 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00656 }
00657 G4double P=std::sqrt(E2-mS2);
00658 G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
00659 G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer();
00660 #ifdef debug
00661 G4cout<<"G4QFR::ChExer: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl;
00662 #endif
00663
00664 G4int PDG=2212;
00665 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;
00666 if(!Z && N==1)
00667 {
00668 Z=1;
00669 N=0;
00670 if (PDG==2212) PDG=2112;
00671 else if(PDG==2112) PDG=2212;
00672 }
00673 G4double xSec=0.;
00674 if(PDG==2212) xSec=PCSmanager->GetCrossSection(false, P, Z, N, PDG);
00675 else xSec=NCSmanager->GetCrossSection(false, P, Z, N, PDG);
00676 #ifdef debug
00677 G4cout<<"G4QuasiFreeRat::ChExer: pPDG="<<pPDG<<",P="<<P<<", CS="<<xSec/millibarn<<G4endl;
00678 #endif
00679 #ifdef nandebug
00680 if(xSec>0. || xSec<0. || xSec==0);
00681 else G4cout<<"*Warning*G4QuasiFreeRatios::ChExer: xSec="<<xSec/millibarn<<G4endl;
00682 #endif
00683
00684 if(xSec <= 0.)
00685 {
00686 #ifdef pdebug
00687 G4cerr<<"-Warning-G4QFR::ChEx:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl;
00688 #endif
00689 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00690 }
00691 G4double mint=0.;
00692 if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);
00693 else mint=NCSmanager->GetExchangeT(Z,N,PDG);
00694 #ifdef pdebug
00695 G4cout<<"G4QFR::ChEx:PDG="<<pPDG<<", P="<<P<<", CS="<<xSec<<", -t="<<mint<<G4endl;
00696 #endif
00697 #ifdef nandebug
00698 if(mint>-.0000001);
00699 else G4cout<<"*Warning*G4QFR::ChExer: -t="<<mint<<G4endl;
00700 #endif
00701 G4double maxt=0.;
00702 if(PDG==2212) maxt=PCSmanager->GetHMaxT();
00703 else maxt=NCSmanager->GetHMaxT();
00704 G4double cost=1.-mint/maxt;
00705 #ifdef pdebug
00706 G4cout<<"G4QuasiFfreeRatio::ChExer: -t="<<mint<<", maxt="<<maxt<<", cost="<<cost<<G4endl;
00707 #endif
00708 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
00709 {
00710 if (cost>1.) cost=1.;
00711 else if(cost<-1.) cost=-1.;
00712 else
00713 {
00714 G4cerr<<"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<",t="<<mint<<",tm="<<maxt<<G4endl;
00715 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00716 }
00717 }
00718 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT);
00719 pr4M=G4LorentzVector(0.,0.,0.,mS);
00720 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
00721 if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost))
00722 {
00723 G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl;
00724
00725 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);
00726 }
00727 #ifdef debug
00728 G4cout<<"G4QFR::ChEx:p4M="<<p4M<<"+r4M="<<reco4M<<"="<<p4M+reco4M<<"="<<tot4M<<G4endl;
00729 #endif
00730 return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt);
00731 }
00732
00733
00734 G4double G4QuasiFreeRatios::ChExElCoef(G4double p, G4int Z, G4int N, G4int pPDG)
00735 {
00736 p/=MeV;
00737 G4double A=Z+N;
00738 if(A<1.5) return 0.;
00739 G4double C=0.;
00740 if (pPDG==2212) C=N/(A+Z);
00741 else if(pPDG==2112) C=Z/(A+N);
00742 else G4cout<<"*Warning*G4QCohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl;
00743 C*=C;
00744
00745 G4double sp=std::sqrt(p);
00746 G4double p2=p*p;
00747 G4double p4=p2*p2;
00748 G4double dl1=std::log(p)-5.;
00749 G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013);
00750 G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p;
00751 G4double R=U/T;
00752 return C*R*R;
00753 }