G4QFreeScattering Class Reference

#include <G4QFreeScattering.hh>


Public Member Functions

 ~G4QFreeScattering ()
std::pair< G4LorentzVector,
G4LorentzVector
Scatter (G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
G4QHadronVectorInElF (G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
std::pair< G4double, G4doubleGetElTotMean (G4double pIU, G4int hPDG, G4int Z, G4int N)
std::pair< G4double, G4doubleGetElTotXS (G4double Mom, G4int PDG, G4bool F)
std::pair< G4double, G4doubleFetchElTot (G4double pGeV, G4int PDG, G4bool F)
std::pair< G4double, G4doubleCalcElTot (G4double pGeV, G4int Index)

Static Public Member Functions

static G4QFreeScatteringGetPointer ()

Protected Member Functions

 G4QFreeScattering ()


Detailed Description

Definition at line 49 of file G4QFreeScattering.hh.


Constructor & Destructor Documentation

G4QFreeScattering::G4QFreeScattering (  )  [protected]

Definition at line 49 of file G4QFreeScattering.cc.

References G4cout, and G4endl.

00050 {
00051 #ifdef pdebug
00052   G4cout<<"***^^^*** G4QFreeScattering singletone is created ***^^^***"<<G4endl;
00053 #endif
00054 }

G4QFreeScattering::~G4QFreeScattering (  ) 

Definition at line 56 of file G4QFreeScattering.cc.

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 }


Member Function Documentation

std::pair< G4double, G4double > G4QFreeScattering::CalcElTot ( G4double  pGeV,
G4int  Index 
)

Definition at line 71 of file G4QFreeScattering.cc.

References FatalException, G4cout, G4endl, G4Exception(), LE, and G4InuclParticleNames::sp.

Referenced by FetchElTot(), G4QuasiFreeRatios::GetElTotXS(), and GetElTotXS().

00072 {
00073   // ---------> Each parameter set can have not more than nPoints=128 parameters
00074   static const G4double lmi=3.5;       // min of (lnP-lmi)^2 parabola
00075   static const G4double pbe=.0557;     // elastic (lnP-lmi)^2 parabola coefficient
00076   static const G4double pbt=.3;        // total (lnP-lmi)^2 parabola coefficient
00077   static const G4double pmi=.1;        // Below that fast LE calculation is made
00078   static const G4double pma=1000.;     // Above that fast HE calculation is made
00079   G4double El=0.;                      // prototype of the elastic hN cross-section
00080   G4double To=0.;                      // prototype of the total hN cross-section
00081   if(p<=0.)
00082   {
00083     //G4cout<<"-Warning-G4QFreeScattering::CalcElTot: p="<<p<<" is 0 or negative"<<G4endl;
00084     return std::make_pair(El,To);
00085   }
00086   if     (!I)                          // pp/nn
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)                        // np/pn
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)                        // pimp/pipn
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;                    // p1
00170       G4double LE=1.53/(lr*lr+.0676);         // p2, p3        
00171       G4double ld=lp-lmi;                     // p4 (lmi=3.5)
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;                     // p5
00177       G4double md=lm*lm+.04;                  // p6
00178       G4double lh=lp-.017;                    // p7
00179       G4double hd=lh*lh+.0025;                // p8
00180       El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd;//p9(pbe=.0557),p10,p11,p12,p13,p14
00181       To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+.4/p4)+1./md+.06/hd;
00182     }
00183   }
00184   else if(I==3)                        // pipp/pimn
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;                   // p1
00205       G4double lr2=lr*lr;
00206       G4double LE=13./(lr2+lr2*lr2+.0676);   // p2, p3
00207       G4double ld=lp-lmi;                    // p4 (lmi=3.5)
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;                    // p5
00213       G4double md=lm*lm+.0576;               // p6
00214       El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md; // p7(pbe=.0557), p8, p9, p10, p11
00215       To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md;
00216     }
00217   }
00218   else if(I==4)                        // Kmp/Kmn/K0p/K0n
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)                        // Kpp/Kpn/aKp/aKn
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)                        // hyperon-N
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)                        // antibaryon-N
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 } // End of CalcElTot

std::pair< G4double, G4double > G4QFreeScattering::FetchElTot ( G4double  pGeV,
G4int  PDG,
G4bool  F 
)

Definition at line 373 of file G4QFreeScattering.cc.

References CalcElTot(), G4cout, G4endl, G4Exception(), G4UniformRand, JustWarning, and CLHEP::detail::n.

Referenced by G4QuasiFreeRatios::GetChExFactor(), G4QuasiFreeRatios::GetElTot(), and GetElTotMean().

00374 {
00375   static const G4int    nlp=300;         // Number of steps in the S(lnp) logTable(5% step)
00376   static const G4int    mlp=nlp+1;       // Number of elements in the S(lnp) logTable
00377   static const G4double lpi=-5.;         // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
00378   static const G4double lpa=10.;         // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
00379   static const G4double mi=std::exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c)
00380   static const G4double ma=std::exp(lpa);// The max p of logTabEl(~ 22. TeV)
00381   static const G4double dl=(lpa-lpi)/nlp;// Step of the logarithmic Table
00382   static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table
00383   //static const G4double toler=.001;      // Relative Tolarence defining "theSameMomentum"
00384   static G4double lastP=0.;              // The last momentum for which XS was calculated
00385   static G4int    lastH=0;               // The last projPDG for which XS was calculated
00386   static G4bool   lastF=true;            // The last nucleon for which XS was calculated
00387   static std::pair<G4double,G4double> lastR(0.,0.); // The last result
00388   // Local Associative Data Base:
00389   static std::vector<G4int>     vI;      // Vector of index for which XS was calculated
00390   static std::vector<G4double>  vM;      // Vector of rel max ln(p) initialized in LogTable
00391   static std::vector<G4int>     vK;      // Vector of topBin number initialized in LogTable
00392   // Last values of the Associative Data Base:
00393   static G4int     lastI=0;              // The Last index for which XS was calculated
00394   static G4double  lastM=0.;             // The Last rel max ln(p) initialized in LogTable
00395   static G4int     lastK=0;              // The Last topBin number initialized in LogTable
00396   static std::pair<G4double,G4double>* lastX=0; // The Last ETPointers to LogTable in heap
00397   // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei
00398   G4int nDB=vI.size();                   // A number of hadrons already initialized in AMDB
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;// VI don't use toler.
00406   //  if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
00407   lastH=PDG;
00408   lastF=F;
00409   G4int ind=-1;                          // Prototipe of the index of the PDG/F combination
00410   // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p),
00411   // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp) 
00412   G4bool kfl=true;                             // Flag of K0/aK0 oscillation
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; // pp/nn
00420   else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
00421   else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
00422   else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
00423   else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
00424   else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
00425   else if ( PDG >  3000 && PDG <  3335) ind=6; // @@ for all hyperons - take Lambda
00426   else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
00427   else
00428   {
00429     // VI changed Fatal to Warning, because this method cannot do better
00430     //    source of problem in classes which make this call
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;  // VI do not use toler
00437   //  if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
00438   if(p<=mi || p>=ma) return CalcElTot(p,ind);   // @@ Slow calculation ! (Warning?)
00439   G4bool found=false;
00440   G4int i=-1;
00441   if(nDB) for (i=0; i<nDB; ++i) if(ind==vI[i])  // Sirch for this index in AMDB
00442   {
00443     found=true;                                 // The index is found
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)                            // Create new line in the AMDB
00451   {
00452 #ifdef pdebug
00453     G4cout<<"G4QFreeScattering::FetchElTot: NewX, ind="<<ind<<", nDB="<<nDB<<G4endl;
00454 #endif
00455     lastX = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot
00456     lastI = ind;                                // Remember the initialized inex
00457     lastK = static_cast<int>((lp-lpi)/dl)+1;    // MaxBin to be initialized in LogTaB
00458     if(lastK>nlp)
00459     {
00460       lastK=nlp;
00461       lastM=lpa-lpi;
00462     }
00463     else lastM = lastK*dl;               // Calculate max initialized ln(p)-lpi for LogTab
00464     G4double pv=mi;
00465     for(G4int j=0; j<=lastK; ++j)        // Calculate LogTab values
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++;                                 // Make a new record to AMDB and position on it
00475     vI.push_back(lastI);
00476     vM.push_back(lastM);
00477     vK.push_back(lastK);
00478     vX.push_back(lastX);
00479   }
00480   else                                   // The A value was found in AMDB
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)              // LogTab must be updated
00492     {
00493       lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTab
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;             // Calculate max initialized ln(p)-lpi for LogTab
00503       G4double pv=std::exp(lpM);         // momentum of the last calculated beam
00504       for(G4int j=nextK; j<=lastK; ++j)  // Calculate LogTab values
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     } // End of LogTab update
00514     if(lastK >= nextK)                   // The AMDB was apdated
00515     {
00516       vM[i]=lastM;
00517       vK[i]=lastK;
00518     }
00519   }
00520   // Now one can use tabeles to calculate the value
00521   G4double dlp=lp-lpi;                       // Shifted log(p) value
00522   G4int n=static_cast<int>(dlp/dl);          // Low edge number of the bin
00523   G4double d=dlp-n*dl;                       // Log shift
00524   G4double e=lastX[n].first;                 // E-Base
00525   lastR.first=e+d*(lastX[n+1].first-e)/dl;   // E-Result
00526   if(lastR.first<0.)  lastR.first = 0.;
00527   G4double t=lastX[n].second;                // T-Base
00528   lastR.second=t+d*(lastX[n+1].second-t)/dl; // T-Result
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 } // End of FetchElTot

std::pair< G4double, G4double > G4QFreeScattering::GetElTotMean ( G4double  pIU,
G4int  hPDG,
G4int  Z,
G4int  N 
)

Definition at line 538 of file G4QFreeScattering.cc.

References FetchElTot(), G4cout, and G4endl.

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;                // To make the result in independent units(IU)
00560   return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
00561 } // End of GetElTotMean

std::pair< G4double, G4double > G4QFreeScattering::GetElTotXS ( G4double  Mom,
G4int  PDG,
G4bool  F 
)

Definition at line 340 of file G4QFreeScattering.cc.

References CalcElTot(), G4cout, G4endl, G4Exception(), G4UniformRand, and JustWarning.

00341 {
00342   G4int ind=0;                                 // Prototype of the reaction index
00343   G4bool kfl=true;                             // Flag of K0/aK0 oscillation
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; // pp/nn
00354   else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
00355   else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
00356   else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
00357   else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
00358   else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
00359   else if ( PDG >  3000 && PDG <  3335) ind=6; // @@ for all hyperons - take Lambda
00360   else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
00361   else {
00362     // VI changed Fatal to Warning, because this method cannot do better
00363     //    source of problem in classes which make this call
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); // This is a slow direct method, better use FetchElTot
00370 }

G4QFreeScattering * G4QFreeScattering::GetPointer (  )  [static]

Definition at line 64 of file G4QFreeScattering.cc.

Referenced by G4QEnvironment::G4QEnvironment(), G4QuasiFreeRatios::G4QuasiFreeRatios(), and G4QEnvironment::operator=().

00065 {
00066   static G4QFreeScattering theQFS;   // *** Static body of the Quasi-Free Scattering ***
00067   return &theQFS;
00068 }

G4QHadronVector * G4QFreeScattering::InElF ( G4int  NPDG,
G4LorentzVector  N4M,
G4int  pPDG,
G4LorentzVector  p4M 
)

Definition at line 667 of file G4QFreeScattering.cc.

References FatalException, G4cout, G4endl, G4Exception(), and G4UniformRand.

00669 {
00670   static const G4double mPi0 = G4QPDGCode(111).GetMass();
00671   static const G4double mPi  = G4QPDGCode(211).GetMass(); // Pi+ (Pi-: -211)
00672   static const G4double mK   = G4QPDGCode(321).GetMass(); // K+  (K- : -321)
00673   static const G4double mK0  = G4QPDGCode(311).GetMass(); // aK0 (K0 : -311)
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   //static const G4double mSigZ= G4QPDGCode(3212).GetMass();
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   //static const G4double mOmM = G4QPDGCode(3334).GetMass();
00683   static const G4double third =1./3.;
00684   static const G4double twothd =2./3.;
00685 
00686   p4M/=megaelectronvolt;   // Convert 4-momenta in MeV (CHIPS units)
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   //G4int N=1;
00698   if(NPDG==2212 || NPDG==90001000)
00699   {
00700     NPDG=2212;
00701     mT=mProt;
00702     Z=1;
00703     //N=0;
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; // Do Nothing Action
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(); // a projectile squared mass
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; // Do Nothing Action
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) // Annihilation/Inelastic of anti-barions not implemented
00738   {
00739     G4cout<<"-Warning-G4QFreeScat::InElF: pPDG="<<pPDG<<G4endl;
00740     return 0; // Do Nothing Action
00741   }
00742   if(pPDG==130 || pPDG==310)
00743   {
00744     if(G4UniformRand()>.5) pPDG = 311; //  K0
00745     else                   pPDG =-311; // aK0
00746   }
00747   G4int mPDG=111;                     // Additional newtral pion by default
00748   G4double mM=mPi0;                   // Default Pi0 mass
00749   G4double r=G4UniformRand();         // A random number to split pi0/pi
00750   if      (pPDG == 2212) // p-N
00751   {
00752     if(Z) // p-p
00753     {
00754       if(r<twothd) // -> n + Pi+ + p
00755       {
00756         pPDG=2112; // n
00757         mP=mNeut;
00758         mPDG= 211; // Pi+
00759         mM=mPi;
00760       }
00761     }
00762     else  // p-n
00763     {
00764       if(r<third) // -> n + Pi+ + n
00765       {
00766         pPDG=2112; // n
00767         mP=mNeut;
00768         mPDG= 211; // Pi+
00769         mM=mPi;
00770       }
00771       else if(r<twothd) // -> p + Pi- + p
00772       {
00773         mPDG=-211; // Pi-
00774         mM=mPi;
00775         NPDG=2212; // p
00776         mT=mProt;
00777       }
00778     }
00779   }
00780   else if (pPDG == 2112) // n-N
00781   {
00782     if(Z) // n-p
00783     {
00784       if(r<third) // -> n + Pi+ + n
00785       {
00786         mPDG= 211; // Pi+
00787         mM=mPi;
00788         NPDG=2112; // n
00789         mT=mNeut;
00790       }
00791       else if(r<twothd) // -> p + Pi- + p
00792       {
00793         pPDG=2212; // p
00794         mP=mProt;
00795         mPDG=-211; // Pi-
00796         mM=mPi;
00797       }
00798     }
00799     else  // n-n
00800     {
00801       if(r<twothd) // -> p + Pi- + n
00802       {
00803         pPDG=2212; // p
00804         mP=mProt;
00805         mPDG= -211; // Pi-
00806         mM=mPi;
00807       }
00808     }
00809   }
00810   else if (pPDG ==  211) // pip-N
00811   {
00812     if(Z) // pip-p
00813     {
00814       if(r<0.5) // -> Pi+ + Pi+ + n
00815       {
00816         mPDG= 211; // Pi+
00817         mM=mPi;
00818         NPDG=2112; // n
00819         mT=mNeut;
00820       }
00821     }
00822     else  // pip-n
00823     {
00824       if(r<third) // -> Pi0 + Pi0 + p
00825       {
00826         pPDG= 111; // Pi0
00827         mP=mPi0;
00828         NPDG=2212; // p
00829         mT=mProt;
00830       }
00831       else if(r<twothd) // -> Pi+ + Pi- + p
00832       {
00833         mPDG=-211; // Pi-
00834         mM=mPi;
00835         NPDG=2212; // p
00836         mT=mProt;
00837       }
00838     }
00839   }
00840   else if (pPDG == -211) // pim-N
00841   {
00842     if(Z) // pim-p
00843     {
00844       if(r<third) // -> Pi0 + Pi0 + n
00845       {
00846         pPDG= 111; // Pi0
00847         mP=mPi0;
00848         NPDG=2112; // n
00849         mT=mNeut;
00850       }
00851       else if(r<twothd) // -> Pi- + Pi+ + n
00852       {
00853         mPDG= 211; // Pi+
00854         mM=mPi;
00855         NPDG=2112; // n
00856         mT=mNeut;
00857       }
00858     }
00859     else  // pim-n
00860     {
00861       if(r<0.5) // -> Pi- + Pi- + p
00862       {
00863         mPDG=-211; // Pi-
00864         mM=mPi;
00865         NPDG=2212; // p
00866         mT=mProt;
00867       }
00868     }
00869   }
00870   else if (pPDG == -321) // Km-N
00871   {
00872     if(Z) // Km-p
00873     {
00874       if(r<0.25) // K0 + Pi0 + n
00875       {
00876         pPDG=-311;  // K0
00877         mP=mK0;
00878         NPDG=2112; // n
00879         mT=mNeut;
00880       }
00881       else if(r<0.5) // -> K- + Pi+ + n
00882       {
00883         mPDG= 211; // Pi+
00884         mM=mPi;
00885         NPDG=2112; // n
00886         mT=mNeut;
00887       }
00888       else if(r<0.75) // -> K0 + Pi- + p
00889       {
00890         pPDG=-311;  // K0
00891         mP=mK0;
00892         mPDG=-211; // Pi-
00893         mM=mPi;
00894       }
00895     }
00896     else  // Km-n
00897     {
00898       if(r<third) // -> K- + Pi- + p
00899       {
00900         mPDG=-211; // Pi-
00901         mM=mPi;
00902         NPDG=2212; // p
00903         mT=mProt;
00904       }
00905       else if(r<twothd) // -> K0 + Pi- + n
00906       {
00907         pPDG=-311;  // K0
00908         mP=mK0;
00909         mPDG=-211; // Pi-
00910         mM=mPi;
00911       }
00912     }
00913   }
00914   else if (pPDG == -311) // K0-N
00915   {
00916     if(Z) // K0-p
00917     {
00918       if(r<third) // K- + Pi+ + p
00919       {
00920         pPDG=-321;  // K-
00921         mP=mK;
00922         NPDG=2212; // p
00923         mT=mProt;
00924       }
00925       else if(r<twothd) // -> K0 + Pi+ + n
00926       {
00927         mPDG= 211; // Pi+
00928         mM=mPi;
00929         NPDG=2112; // n
00930         mT=mNeut;
00931       }
00932     }
00933     else  // K0-n
00934     {
00935       if(r<0.25) // -> K- + Pi+ + n
00936       {
00937         pPDG=-321; // K-
00938         mP=mK;
00939         mPDG= 211; // Pi+
00940         mM=mPi;
00941       }
00942       else if(r<0.5) // -> K- + Pi0 + p
00943       {
00944         pPDG=-321; // K-
00945         mP=mK;
00946         NPDG=2212; // p
00947         mT=mProt;
00948       }
00949       else if(r<0.75) // -> K0 + Pi- + p
00950       {
00951         mPDG=-211; // Pi-
00952         mM=mPi;
00953         NPDG=2212; // p
00954         mT=mProt;
00955       }
00956     }
00957   }
00958   else if (pPDG ==  321) // Kp-N
00959   {
00960     if(Z) // Kp-p
00961     {
00962       if(r<third) // -> K+ + Pi+ + n
00963       {
00964         mPDG= 211; // Pi+
00965         mM=mPi;
00966         NPDG=2112; // n
00967         mT=mNeut;
00968       }
00969       else if(r<twothd) // -> aK0 + Pi+ + p
00970       {
00971         pPDG= 311; // aK0
00972         mP=mK0;
00973         mPDG= 211; // Pi+
00974         mM=mPi;
00975       }
00976     }
00977     else  // Kp-n
00978     {
00979       if(r<0.25) // -> aK0 + Pi+ + n
00980       {
00981         pPDG= 311; // aK0
00982         mP=mK0;
00983         mPDG= 211; // Pi+
00984         mM=mPi;
00985       }
00986       else if(r<0.5) // -> Kp + Pi- + p
00987       {
00988         mPDG=-211; // Pi-
00989         mM=mPi;
00990         NPDG=2212; // p
00991         mT=mProt;
00992       }
00993       else if(r<0.75) // -> aK0 + Pi0 + p
00994       {
00995         pPDG= 311; // aK0
00996         mP=mK0;
00997         NPDG=2212; // p
00998         mT=mProt;
00999       }
01000     }
01001   }
01002   else if (pPDG ==  311) // aK0-N
01003   {
01004     if(Z) // aK0-p
01005     {
01006       if(r<0.25) // -> K+ + Pi- + p
01007       {
01008         pPDG= 321; // K+
01009         mP=mK;
01010         mPDG=-211; // Pi-
01011         mM=mPi;
01012       }
01013       else if(r<0.5) // -> aK0 + Pi+ + n
01014       {
01015         mPDG= 211; // Pi+
01016         mM=mPi;
01017         NPDG=2112; // n
01018         mT=mNeut;
01019       }
01020       else if(r<0.75) // -> K+ + Pi0 + n
01021       {
01022         pPDG= 321; // K+
01023         mP=mK;
01024         NPDG=2112; // n
01025         mT=mNeut;
01026       }
01027     }
01028     else  // aK0-n
01029     {
01030       if(r<third) // -> aK0 + Pi- + p
01031       {
01032         mPDG=-211; // Pi-
01033         mM=mPi;
01034         NPDG=2212; // p
01035         mT=mProt;
01036       }
01037       else if(r<twothd) // -> K+ + Pi- + n
01038       {
01039         pPDG= 321; // K+
01040         mP=mK;
01041         mPDG=-211; // Pi-
01042         mM=mPi;
01043       }
01044     }
01045   }
01046   else if (pPDG == 3122 || pPDG== 3212) // Lambda/Sigma0-N
01047   {
01048     if(pPDG == 3212)
01049     {
01050       pPDG=3122;
01051       mP=mLamb;
01052     }
01053     if(Z) // L/S0-p
01054     {
01055       if(r<0.2) // -> SigP + Pi0 + n
01056       {
01057         pPDG=3222; // SigP
01058         mP=mSigP;
01059         NPDG=2112; // n
01060         mT=mNeut;
01061       }
01062       else if(r<0.4) // -> SigP + Pi- + p
01063       {
01064         pPDG=3222; // SigP
01065         mP=mSigP;
01066         mPDG=-211; // Pi-
01067         mM=mPi;
01068       }
01069       else if(r<0.6) // -> SigM + Pi+ + p
01070       {
01071         pPDG=3112; // SigM
01072         mP=mSigM;
01073         mPDG= 211; // Pi+
01074         mM=mPi;
01075       }
01076       else if(r<0.8) // -> Lamb + Pi+ + n
01077       {
01078         mPDG= 211; // Pi+
01079         mM=mPi;
01080         NPDG=2112; // n
01081         mT=mNeut;
01082       }
01083     }
01084     else  // L/S0-n
01085     {
01086       if(r<0.2) // -> SigM + Pi0 + p
01087       {
01088         pPDG=3112; // SigM
01089         mP=mSigM;
01090         NPDG=2212; // p
01091         mT=mProt;
01092       }
01093       else if(r<0.4) // -> SigP + Pi- + n
01094       {
01095         pPDG=3222; // SigP
01096         mP=mSigP;
01097         mPDG=-211; // Pi-
01098         mM=mPi;
01099       }
01100       else if(r<0.6) // -> SigM + Pi+ + n
01101       {
01102         pPDG=3112; // SigM
01103         mP=mSigM;
01104         mPDG= 211; // Pi+
01105         mM=mPi;
01106       }
01107       else if(r<0.8) // -> Lamb + Pi- + p
01108       {
01109         mPDG=-211; // Pi-
01110         mM=mPi;
01111         NPDG=2212; // p
01112         mT=mProt;
01113       }
01114     }
01115   }
01116   else if (pPDG == 3112) // Sigma- -N
01117   {
01118     if(Z) // SigM-p
01119     {
01120       if(r<0.25) // -> Lamb + Pi0 + n
01121       {
01122         pPDG=3122; // Lamb
01123         mP=mLamb;
01124         NPDG=2112; // n
01125         mT=mNeut;
01126       }
01127       else if(r<0.5) // -> Lamb + Pi- + p
01128       {
01129         pPDG=3122; // Lamb
01130         mP=mLamb;
01131         mPDG=-211; // Pi-
01132         mM=mPi;
01133       }
01134       else if(r<0.75) // -> SigM + Pi+ + n
01135       {
01136         mPDG= 211; // Pi+
01137         mM=mPi;
01138         NPDG=2112; // n
01139         mT=mNeut;
01140       }
01141     }
01142     else  // SigM-n
01143     {
01144       if(r<third) // -> Lamb + Pi- + n
01145       {
01146         pPDG=3122; // Lamb
01147         mP=mLamb;
01148         mPDG=-211; // Pi-
01149         mM=mPi;
01150       }
01151       else if(r<twothd) // -> SigM + Pi- + p
01152       {
01153         mPDG=-211; // Pi-
01154         mM=mPi;
01155         NPDG=2212; // p
01156         mT=mProt;
01157       }
01158     }
01159   }
01160   else if (pPDG == 3222) // Sigma+ -N
01161   {
01162     if(Z) // SigP-p
01163     {
01164       if(r<third) // -> Lamb + Pi+ + p
01165       {
01166         pPDG=3122; // Lamb
01167         mP=mLamb;
01168         mPDG= 211; // Pi+
01169         mM=mPi;
01170       }
01171       else if(r<twothd) // -> SigP + Pi+ + n
01172       {
01173         mPDG= 211; // Pi+
01174         mM=mPi;
01175         NPDG=2112; // n
01176         mT=mNeut;
01177       }
01178     }
01179     else  // SigP-n
01180     {
01181       if(r<0.25) // -> Lamb + Pi0 + p
01182       {
01183         pPDG=3122; // Lamb
01184         mP=mLamb;
01185         NPDG=2212; // p
01186         mT=mProt;
01187       }
01188       else if(r<0.5) // -> Lamb + Pi+ + n
01189       {
01190         pPDG=3122; // Lamb
01191         mP=mLamb;
01192         mPDG= 211; // Pi+
01193         mM=mPi;
01194       }
01195       else if(r<0.75) // -> SigP + Pi- + p
01196       {
01197         mPDG=-211; // Pi-
01198         mM=mPi;
01199         NPDG=2212; // p
01200         mT=mProt;
01201       }
01202     }
01203   }
01204   else if (pPDG == 3312) // Xi- -N
01205   {
01206     if(Z) // XiM-p
01207     {
01208       if(r<0.25) // -> Xi0 + Pi0 + n
01209       {
01210         pPDG=3322; // Xi0
01211         mP=mXiZ;
01212         NPDG=2112; // n
01213         mT=mNeut;
01214       }
01215       else if(r<0.5) // -> Xi0 + Pi- + p
01216       {
01217         pPDG=3322; // Xi0
01218         mP=mXiZ;
01219         mPDG=-211; // Pi-
01220         mM=mPi;
01221       }
01222       else if(r<0.75) // -> XiM + Pi+ + n
01223       {
01224         mPDG= 211; // Pi+
01225         mM=mPi;
01226         NPDG=2112; // n
01227         mT=mNeut;
01228       }
01229     }
01230     else  // XiM-n
01231     {
01232       if(r<third) // -> Xi0 + Pi- + n
01233       {
01234         pPDG=3322; // Xi0
01235         mP=mXiZ;
01236         mPDG=-211; // Pi-
01237         mM=mPi;
01238       }
01239       else if(r<twothd) // -> XiM + Pi- + p
01240       {
01241         mPDG=-211; // Pi-
01242         mM=mPi;
01243         NPDG=2212; // p
01244         mT=mProt;
01245       }
01246     }
01247   }
01248   else if (pPDG == 3322) // Xi0 -N
01249   {
01250     if(Z) // Xi0-p
01251     {
01252       if(r<third) // -> Xi- + Pi+ + p
01253       {
01254         pPDG=3312; // Xi-
01255         mP=mXiM;
01256         mPDG= 211; // Pi+
01257         mM=mPi;
01258       }
01259       else if(r<twothd) // -> Xi0 + Pi+ + n
01260       {
01261         mPDG= 211; // Pi+
01262         mM=mPi;
01263         NPDG=2112; // n
01264         mT=mNeut;
01265       }
01266     }
01267     else  // Xi0-n
01268     {
01269       if(r<0.25) // -> Xi- + Pi0 + p
01270       {
01271         pPDG=3312; // Xi-
01272         mP=mXiM;
01273         NPDG=2212; // p
01274         mT=mProt;
01275       }
01276       else if(r<0.5) // -> Xi- + Pi+ + n
01277       {
01278         pPDG=3312; // Xi-
01279         mP=mXiM;
01280         mPDG= 211; // Pi+
01281         mM=mPi;
01282       }
01283       else if(r<0.75) // -> Xi0 + Pi- + p
01284       {
01285         mPDG=-211; // Pi-
01286         mM=mPi;
01287         NPDG=2212; // p
01288         mT=mProt;
01289       }
01290     }
01291   }
01292   else if (pPDG == 3334) // Om- -N
01293   {
01294     if(Z) // OmM-p
01295     {
01296       if(r<0.5)    // -> Om- + Pi+ + n
01297       {
01298         mPDG= 211; // Pi+
01299         mM=mPi;
01300         NPDG=2112; // n
01301         mT=mNeut;
01302       } 
01303     }
01304     else  // OmM-n
01305     {
01306       if(r<0.5)    // -> Om- + Pi- + p
01307       {
01308         mPDG=-211; // Pi-
01309         mM=mPi;
01310         NPDG=2212; // p
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; // Proto of the Result
01323   G4LorentzVector m4M(0.,0.,0.,mM);
01324   if(mC-mP-mM-mT < 0.000001) // Equal share
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 /* << " <? " << mT + mP + mM */ << G4endl;
01338       G4Exception("G4QFreeScattering::InElF()","HAD_CHPS_0027", FatalException, ed);
01339     }
01340     G4QHadron* h1 = new G4QHadron(pPDG, p4M);
01341     TripQH->push_back(h1); // (delete equivalent, responsibility of users)
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); // (delete equivalent, responsibility of users)
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); // (delete equivalent, responsibility of users)
01352 #ifdef debug
01353     G4cout << "G4QFreeScat::InElF: H3=" << NPDG << N4M << G4endl;
01354 #endif
01355   }
01356   return TripQH; // The Result
01357 } // End of Scatter

std::pair< G4LorentzVector, G4LorentzVector > G4QFreeScattering::Scatter ( G4int  NPDG,
G4LorentzVector  N4M,
G4int  pPDG,
G4LorentzVector  p4M 
)

Definition at line 565 of file G4QFreeScattering.cc.

References FatalException, G4cerr, G4cout, G4endl, G4Exception(), and G4UniformRand.

00567 {
00568   static const G4double mNeut= G4QPDGCode(2112).GetMass();
00569   static const G4double mProt= G4QPDGCode(2212).GetMass();
00570   p4M/=megaelectronvolt;   // Convert 4-momenta in MeV (CHIPS units)
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     //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Use this if not exception
00597   }
00598   G4double mT2=mT*mT;    // a squared mass of the free scattered nuclead cluster (FSNC)
00599   G4double mP2=p4M.m2(); // a projectile squared mass
00600   G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT); // a projectile energy in the CMS of FSNC
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); // Do Nothing Action
00611   }
00612   G4double pP2=E2-mP2; // Squared Momentum in pseudo laboratory system (final particles)
00613   G4double P=std::sqrt(pP2); // Momentum in pseudo laboratory system (final particles) ?
00614 #ifdef ppdebug
00615   G4cout<<"G4QFreeS::Scatter: Before XS, P="<<P<<",Z="<<Z<<",N="<<N<<",PDG="<<pPDG<<G4endl;
00616 #endif
00617   // @@ Temporary NN t-dependence for all hadrons
00618   if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QFreeScat::Scatter: pPDG="<<pPDG<<G4endl;
00619   // @@ check a possibility to separate p, n, or alpha (!)
00620   G4double dmT=mT+mT;
00621   G4double s_value=dmT*E2+mP2+mT2;         // Mondelstam s (?)
00622   G4double maxt=dmT*dmT*pP2/s_value;       // max possible |t|
00623   G4double mint=0.;                        // Prototype of functional rand -t (MeV^2)
00624   if(P < 14.) mint=maxt*G4UniformRand();   // S-wave
00625   else                                     // Calculate slopes (no cash !)
00626   {
00627     G4double p4=pP2*pP2/1000000.;          // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00628     G4double theB1=(7.2+4.32/p4/(p4+12./P))/(1.+2.5/p4); // p4 in GeV, P in MeV
00629     mint=-std::log(G4UniformRand())/theB1; // t-chan only
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;        // cos(theta) in CMS
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); // Do Nothing Action
00651     }
00652   }
00653   G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT);      // 4mom of the recoil nucleon
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     //G4Exception("G4QFS::Scat:","009",FatalException,"Decay of ElasticComp");
00659     return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
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 ); // The Result
00665 } // End of Scatter


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:07 2013 for Geant4 by  doxygen 1.4.7