Geant4: G4QuasiElRatios Class Reference

G4QuasiElRatios::G4QuasiElRatios ( ) [protected]

Definition at line 60 of file G4QuasiElRatios.cc.

References G4ChipsNeutronElasticXS::Default_Name(), G4ChipsProtonElasticXS::Default_Name(), G4CrossSectionDataSetRegistry::GetCrossSectionDataSet(), and G4CrossSectionDataSetRegistry::Instance().

00061 {
00062     
00063     PCSmanager=(G4ChipsProtonElasticXS*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsProtonElasticXS::Default_Name());
00064     
00065     NCSmanager=(G4ChipsNeutronElasticXS*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsNeutronElasticXS::Default_Name());
00066 }

G4QuasiElRatios::~G4QuasiElRatios ( )

Definition at line 68 of file G4QuasiElRatios.cc.

00069 {
00070     std::vector<G4double*>::iterator pos;
00071     for(pos=vT.begin(); pos<vT.end(); pos++)
00072     { delete [] *pos; }
00073     vT.clear();
00074     for(pos=vL.begin(); pos<vL.end(); pos++)
00075     { delete [] *pos; }
00076     vL.clear();
00077     
00078     std::vector<std::pair<G4double,G4double>*>::iterator pos2;
00079     for(pos2=vX.begin(); pos2<vX.end(); pos2++)
00080     { delete [] *pos2; }
00081     vX.clear();
00082 }

G4double G4QuasiElRatios::ChExElCoef	(	G4double	p,
		G4int	Z,
		G4int	N,
		G4int	pPDG
	)

Definition at line 1005 of file G4QuasiElRatios.cc.

References G4cout, G4endl, and G4InuclParticleNames::sp.

01006 {
01007     p/=MeV;                                // Converted from independent units
01008     G4double A=Z+N;
01009     if(A<1.5) return 0.;
01010     G4double C=0.;
01011     if     (pPDG==2212) C=N/(A+Z);
01012     else if(pPDG==2112) C=Z/(A+N);
01013     else G4cout<<"*Warning*G4CohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl;
01014     C*=C;                         // Coherent processes squares the amplitude
01015     // @@ This is true only for nucleons: other projectiles must be treated differently
01016     G4double sp=std::sqrt(p);
01017     G4double p2=p*p;            
01018     G4double p4=p2*p2;
01019     G4double dl1=std::log(p)-5.;
01020     G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013);
01021     G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p; 
01022     G4double R=U/T;
01023     return C*R*R;
01024 }

std::pair< G4LorentzVector, G4LorentzVector > G4QuasiElRatios::ChExer	(	G4int	NPDG,
		G4LorentzVector	N4M,
		G4int	pPDG,
		G4LorentzVector	p4M
	)

Definition at line 893 of file G4QuasiElRatios.cc.

References FatalException, G4cerr, G4cout, G4endl, G4Exception(), G4UniformRand, G4ChipsNeutronElasticXS::GetChipsCrossSection(), G4ChipsProtonElasticXS::GetChipsCrossSection(), G4ChipsNeutronElasticXS::GetExchangeT(), G4ChipsProtonElasticXS::GetExchangeT(), G4ChipsNeutronElasticXS::GetHMaxT(), G4ChipsProtonElasticXS::GetHMaxT(), G4ParticleDefinition::GetPDGMass(), G4Neutron::Neutron(), G4Proton::Proton(), and RelDecayIn2().

00895 {
00896     static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass();
00897     static const G4double mProt= G4Proton::Proton()->GetPDGMass();
00898     G4LorentzVector pr4M=p4M/megaelectronvolt;          // Convert 4-momenta in MeV(keep p4M)
00899     N4M/=megaelectronvolt;
00900     G4LorentzVector tot4M=N4M+p4M;
00901     G4int Z=0;
00902     G4int N=1;
00903     G4int sPDG=0;                                        // PDG code of the scattered hadron
00904     G4double mS=0.;                                      // proto of mass of scattered hadron
00905     G4double mT=mProt;                                   // mass of the recoil nucleon
00906     if(NPDG==2212)
00907     {
00908         mT=mNeut;
00909         Z=1;
00910         N=0;
00911         if(pPDG==-211) sPDG=111;                           // pi+    -> pi0
00912         else if(pPDG==-321)
00913         {
00914             sPDG=310;                                        // K+     -> K0S
00915             if(G4UniformRand()>.5) sPDG=130;                 // K+     -> K0L
00916         }
00917         else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321;  // K0     -> K+ (?)
00918         else if(pPDG==3112) sPDG=3212;                     // Sigma- -> Sigma0
00919         else if(pPDG==3212) sPDG=3222;                     // Sigma0 -> Sigma+
00920         else if(pPDG==3312) sPDG=3322;                     // Xi-    -> Xi0
00921     }
00922     else if(NPDG==2112) // Default
00923     {
00924         if(pPDG==211)  sPDG=111;                           // pi+    -> pi0
00925         else if(pPDG==321)
00926         {
00927             sPDG=310;                                        // K+     -> K0S
00928             if(G4UniformRand()>.5) sPDG=130;                 // K+     -> K0L
00929         }
00930         else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321; // K0     -> K- (?)
00931         else if(pPDG==3222) sPDG=3212;                     // Sigma+ -> Sigma0
00932         else if(pPDG==3212) sPDG=3112;                     // Sigma0 -> Sigma-
00933         else if(pPDG==3322) sPDG=3312;                     // Xi0    -> Xi-
00934     }
00935     else
00936     {
00937         G4cout<<"Error:G4QuasiElRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
00938         G4Exception("G4QuasiElRatios::ChExer:","21",FatalException,"QE complain");
00939         //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
00940     }
00941     if(sPDG) mS=mNeut;
00942     else
00943     {
00944         G4cout<<"Error:G4QuasiElRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl;
00945         G4Exception("G4QuasiElRatios::ChExer:","21",FatalException,"QE complain");
00946         //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
00947     }
00948     G4double mT2=mT*mT;
00949     G4double mS2=mS*mS;
00950     G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);
00951     G4double E2=E*E;
00952     if(E<0. || E2<mS2)
00953     {
00954         return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
00955     }
00956     G4double P=std::sqrt(E2-mS2);                   // Momentum in pseudo laboratory system
00957     // @@ Temporary NN t-dependence for all hadrons
00958     G4int PDG=2212;                                                // *TMP* instead of pPDG
00959     if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;               // *TMP* instead of pPDG
00960     if(!Z && N==1)                 // Change for Quasi-Elastic on neutron
00961     {
00962         Z=1;
00963         N=0;
00964         if     (PDG==2212) PDG=2112;
00965         else if(PDG==2112) PDG=2212;
00966     }
00967     G4double xSec=0.;                        // Prototype of Recalculated Cross Section *TMP*
00968     if(PDG==2212) xSec=PCSmanager->GetChipsCrossSection(P, Z, N, PDG); // P CrossSect *TMP*
00969     else          xSec=NCSmanager->GetChipsCrossSection(P, Z, N, PDG); // N CrossSect *TMP*
00970     // @@ check a possibility to separate p, n, or alpha (!)
00971     if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
00972     {
00973         return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
00974     }
00975     G4double mint=0.;                        // Prototype of functional rand -t (MeV^2) *TMP*
00976     if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
00977     else          mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
00978     G4double maxt=0.;                                    // Prototype of max possible -t
00979     if(PDG==2212) maxt=PCSmanager->GetHMaxT();           // max possible -t
00980     else          maxt=NCSmanager->GetHMaxT();           // max possible -t
00981     G4double cost=1.-mint/maxt;                          // cos(theta) in CMS
00982     if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
00983     {
00984         if     (cost>1.)  cost=1.;
00985         else if(cost<-1.) cost=-1.;
00986         else
00987         {
00988             G4cerr<<"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<",t="<<mint<<",tm="<<maxt<<G4endl;
00989             return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
00990         }
00991     }
00992     G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT);      // 4mom of the recoil nucleon
00993     pr4M=G4LorentzVector(0.,0.,0.,mS);                        // 4mom of the scattered hadron
00994     G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
00995     if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
00996     {
00997         G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl;
00998         //G4Exception("G4QFR::ChExer:","009",FatalException,"Decay of ElasticComp");
00999         return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
01000     }
01001     return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
01002 } // End of ChExer

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

Definition at line 588 of file G4QuasiElRatios.cc.

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

Referenced by GetChExFactor(), and GetElTot().

00589 {
00590     static const G4int    nlp=300;         // Number of steps in the S(lnp) logTable(5% step)
00591     static const G4int    mlp=nlp+1;       // Number of elements in the S(lnp) logTable
00592     static const G4double lpi=-5.;         // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
00593     static const G4double lpa=10.;         // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
00594     static const G4double mi=std::exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c)
00595     static const G4double ma=std::exp(lpa);// The max p of logTabEl(~ 22. TeV)
00596     static const G4double dl=(lpa-lpi)/nlp;// Step of the logarithmic Table
00597     static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table
00598     //static const G4double toler=.001;      // Relative Tolarence defining "theSameMomentum"
00599     static G4double lastP=0.;              // The last momentum for which XS was calculated
00600     static G4int    lastH=0;               // The last projPDG for which XS was calculated
00601     static G4bool   lastF=true;            // The last nucleon for which XS was calculated
00602     static std::pair<G4double,G4double> lastR(0.,0.); // The last result
00603     // Local Associative Data Base:
00604     static std::vector<G4int>     vI;      // Vector of index for which XS was calculated
00605     static std::vector<G4double>  vM;      // Vector of rel max ln(p) initialized in LogTable
00606     static std::vector<G4int>     vK;      // Vector of topBin number initialized in LogTable
00607     // Last values of the Associative Data Base:
00608     static G4int     lastI=0;              // The Last index for which XS was calculated
00609     static G4double  lastM=0.;             // The Last rel max ln(p) initialized in LogTable
00610     static G4int     lastK=0;             // The Last topBin number initialized in LogTable
00611     static std::pair<G4double,G4double>* lastX=0; // The Last ETPointers to LogTable in heap
00612     // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei
00613     G4int nDB=vI.size();                   // A number of hadrons already initialized in AMDB
00614     if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP) return lastR;// VI don't use toler.
00615     //  if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
00616     lastH=PDG;
00617     lastF=F;
00618     G4int ind=-1;                          // Prototipe of the index of the PDG/F combination
00619     // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p),
00620     // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp) 
00621     G4bool kfl=true;                             // Flag of K0/aK0 oscillation
00622     G4bool kf=false;
00623     if(PDG==130||PDG==310)
00624     {
00625         kf=true;
00626         if(G4UniformRand()>.5) kfl=false;
00627     }
00628     if      ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
00629     else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
00630     else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
00631     else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
00632     else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
00633     else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
00634     else if ( PDG >  3000 && PDG <  3335) ind=6; // @@ for all hyperons - take Lambda
00635     else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
00636     else {
00637         G4cout<<"*Error*G4QuasiElRatios::FetchElTot: PDG="<<PDG
00638         <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
00639         G4Exception("G4QuasiELRatio::FetchElTot:","22",FatalException,"QECrash");
00640     }
00641     if(nDB && lastI==ind && p>0. && p==lastP) return lastR;  // VI do not use toler
00642     //  if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
00643     if(p<=mi || p>=ma) return CalcElTot(p,ind);   // @@ Slow calculation ! (Warning?)
00644     G4bool found=false;
00645     G4int i=-1;
00646     if(nDB) for (i=0; i<nDB; i++) if(ind==vI[i])  // Sirch for this index in AMDB
00647     {
00648         found=true;                                 // The index is found
00649         break;
00650     }
00651     G4double lp=std::log(p);
00652     if(!nDB || !found)                            // Create new line in the AMDB
00653     {
00654         lastX = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot
00655         lastI = ind;                                // Remember the initialized inex
00656         lastK = static_cast<int>((lp-lpi)/dl)+1;    // MaxBin to be initialized in LogTaB
00657         if(lastK>nlp)
00658         {
00659             lastK=nlp;
00660             lastM=lpa-lpi;
00661         }
00662         else lastM = lastK*dl;               // Calculate max initialized ln(p)-lpi for LogTab
00663         G4double pv=mi;
00664         for(G4int j=0; j<=lastK; j++)        // Calculate LogTab values
00665         {
00666             lastX[j]=CalcElTot(pv,ind);
00667             if(j!=lastK) pv*=edl;
00668         }
00669         i++;                                 // Make a new record to AMDB and position on it
00670         vI.push_back(lastI);
00671         vM.push_back(lastM);
00672         vK.push_back(lastK);
00673         vX.push_back(lastX);
00674     }
00675     else                                   // The A value was found in AMDB
00676     {
00677         lastI=vI[i];
00678         lastM=vM[i];
00679         lastK=vK[i];
00680         lastX=vX[i];
00681         G4int nextK=lastK+1;
00682         G4double lpM=lastM+lpi;
00683         if(lp>lpM && lastK<nlp)              // LogTab must be updated
00684         {
00685             lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTab
00686             if(lastK>nlp)
00687             {
00688                 lastK=nlp;
00689                 lastM=lpa-lpi;
00690             }
00691             else lastM = lastK*dl;           // Calculate max initialized ln(p)-lpi for LogTab
00692             G4double pv=std::exp(lpM);       // momentum of the last calculated beam
00693             for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values
00694             {
00695                 pv*=edl;
00696                 lastX[j]=CalcElTot(pv,ind);
00697             }
00698         } // End of LogTab update
00699         if(lastK>=nextK)                   // The AMDB was apdated
00700         {
00701             vM[i]=lastM;
00702             vK[i]=lastK;
00703         }
00704     }
00705     // Now one can use tabeles to calculate the value
00706     G4double dlp=lp-lpi;                       // Shifted log(p) value
00707     G4int n=static_cast<int>(dlp/dl);          // Low edge number of the bin
00708     G4double d=dlp-n*dl;                       // Log shift
00709     G4double e=lastX[n].first;                 // E-Base
00710     lastR.first=e+d*(lastX[n+1].first-e)/dl;   // E-Result
00711     if(lastR.first<0.)  lastR.first = 0.;
00712     G4double t=lastX[n].second;                // T-Base
00713     lastR.second=t+d*(lastX[n+1].second-t)/dl; // T-Result
00714     if(lastR.second<0.) lastR.second= 0.;
00715     if(lastR.first>lastR.second) lastR.first = lastR.second;
00716     return lastR;
00717 } // End of FetchElTot

std::pair< G4double, G4double > G4QuasiElRatios::GetChExFactor	(	G4double	pIU,
		G4int	pPDG,
		G4int	Z,
		G4int	N
	)

Definition at line 736 of file G4QuasiElRatios.cc.

References FetchElTot(), G4cout, and G4endl.

00738 {
00739     G4double pGeV=pIU/gigaelectronvolt;
00740     G4double resP=0.;
00741     G4double resN=0.;
00742     if(Z<1 && N<1)
00743     {
00744         G4cout<<"-Warning-G4QuasiElRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
00745         return std::make_pair(resP,resN);
00746     }
00747     G4double A=Z+N;
00748     G4double pf=0.;                              // Possibility to interact with a proton
00749     G4double nf=0.;                              // Possibility to interact with a neutron
00750     if   (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+N);
00751     else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z);
00752     else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310)
00753     {
00754         G4double dA=A+A;
00755         pf=Z/(dA+N+N);
00756         nf=N/(dA+Z+Z);
00757     }
00758     G4double mult=1.;  // Factor of increasing multiplicity ( ? @@)
00759     if(pGeV>.5)
00760     {
00761         mult=1./(1.+std::log(pGeV+pGeV))/pGeV;
00762         if(mult>1.) mult=1.;
00763     }
00764     if(pf)
00765     {
00766         std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
00767         resP=pf*(hp.second/hp.first-1.)*mult;
00768     }
00769     if(nf)
00770     {
00771         std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
00772         resN=nf*(hn.second/hn.first-1.)*mult;
00773     }
00774     return std::make_pair(resP,resN);
00775 } // End of GetChExFactor

std::pair< G4double, G4double > G4QuasiElRatios::GetElTot	(	G4double	pIU,
		G4int	hPDG,
		G4int	Z,
		G4int	N
	)

Definition at line 720 of file G4QuasiElRatios.cc.

References FetchElTot(), G4cout, and G4endl.

Referenced by GetRatios().

00722 {
00723     G4double pGeV=pIU/gigaelectronvolt;
00724     if(Z<1 && N<1)
00725     {
00726         G4cout<<"-Warning-G4QuasiElRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
00727         return std::make_pair(0.,0.);
00728     }
00729     std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
00730     std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
00731     G4double A=(Z+N)/millibarn;                // To make the result in independent units(IU)
00732     return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
00733 } // End of GetElTot

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

Definition at line 561 of file G4QuasiElRatios.cc.

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

00562 {
00563     G4int ind=0;                                 // Prototype of the reaction index
00564     G4bool kfl=true;                             // Flag of K0/aK0 oscillation
00565     G4bool kf=false;
00566     if(PDG==130||PDG==310)
00567     {
00568         kf=true;
00569         if(G4UniformRand()>.5) kfl=false;
00570     }
00571     if      ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
00572     else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
00573     else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
00574     else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
00575     else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
00576     else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
00577     else if ( PDG >  3000 && PDG <  3335) ind=6; // @@ for all hyperons - take Lambda
00578     else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
00579     else {
00580         G4cout<<"*Error*G4QuasiElRatios::CalcElTotXS: PDG="<<PDG
00581         <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
00582         G4Exception("G4QuasiElRatio::CalcElTotXS:","22",FatalException,"QEcrash");
00583     }
00584     return CalcElTot(p,ind);
00585 }

G4QuasiElRatios * G4QuasiElRatios::GetPointer ( ) [static]

Definition at line 85 of file G4QuasiElRatios.cc.

00086 {
00087     static G4QuasiElRatios theRatios;   // *** Static body of the QEl Cross Section ***
00088     return &theRatios;
00089 }

std::pair< G4double, G4double > G4QuasiElRatios::GetRatios	(	G4double	pIU,
		G4int	prPDG,
		G4int	tgZ,
		G4int	tgN
	)

Definition at line 92 of file G4QuasiElRatios.cc.

References GetElTot().

Referenced by G4QuasiElasticChannel::GetFraction().

00094 {
00095     G4double R=0.;
00096     G4double QF2In=1.;                        // Prototype of QuasiFree/Inel ratio for hN_tot
00097     G4int tgA=tgZ+tgN;
00098     if(tgA<2) return std::make_pair(QF2In,R); // No quasi-elastic on the only nucleon
00099     std::pair<G4double,G4double> ElTot=GetElTot(pIU, pPDG, tgZ, tgN); // mean hN El&Tot(IU)
00100     //if( ( (pPDG>999 && pIU<227.) || pIU<27.) && tgA>1) R=1.; // @@ TMP to accelerate @lowE
00101     if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.;                // To accelerate @lowE
00102     else if(ElTot.second>0.)
00103     {
00104         R=ElTot.first/ElTot.second;             // El/Total ratio (does not depend on units
00105         QF2In=GetQF2IN_Ratio(ElTot.second/millibarn, tgZ+tgN);   // QuasiFree/Inelastic ratio
00106     }
00107     return std::make_pair(QF2In,R);
00108 }

G4bool G4QuasiElRatios::RelDecayIn2	(	G4LorentzVector &	theMomentum,
		G4LorentzVector &	f4Mom,
		G4LorentzVector &	s4Mom,
		G4LorentzVector &	dir,
		G4double	maxCost = `1.`,
		G4double	minCost = `-1.`
	)

Definition at line 1027 of file G4QuasiElRatios.cc.

References G4cerr, G4endl, and G4UniformRand.

Referenced by ChExer(), and Scatter().

01029 {
01030     G4double fM2 = f4Mom.m2();
01031     G4double fM  = std::sqrt(fM2);              // Mass of the 1st Hadron
01032     G4double sM2 = s4Mom.m2();
01033     G4double sM  = std::sqrt(sM2);              // Mass of the 2nd Hadron
01034     G4double iM2 = theMomentum.m2();
01035     G4double iM  = std::sqrt(iM2);              // Mass of the decaying hadron
01036     G4double vP  = theMomentum.rho();      // Momentum of the decaying hadron
01037     G4double dE  = theMomentum.e();        // Energy of the decaying hadron
01038     if(dE<vP)
01039     {
01040         G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
01041         G4double accuracy=.000001*vP;
01042         G4double emodif=std::fabs(dE-vP);
01043         //if(emodif<accuracy)
01044         //{
01045         G4cerr<<"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
01046         theMomentum.setE(vP+emodif+.01*accuracy);
01047         //}
01048     }
01049     G4ThreeVector ltb = theMomentum.boostVector();// Boost vector for backward Lorentz Trans.
01050     G4ThreeVector ltf = -ltb;              // Boost vector for forward Lorentz Trans.
01051     G4LorentzVector cdir = dir;            // A copy to make a transformation to CMS
01052     cdir.boost(ltf);                       // Direction transpormed to CMS of the Momentum
01053     G4ThreeVector vdir = cdir.vect();      // 3-Vector of the direction-particle
01054     G4ThreeVector vx(0.,0.,1.);            // Ort in the direction of the reference particle
01055     G4ThreeVector vy(0.,1.,0.);            // First ort orthogonal to the direction
01056     G4ThreeVector vz(1.,0.,0.);            // Second ort orthoganal to the direction
01057     if(vdir.mag2() > 0.)                   // the refference particle isn't at rest in CMS
01058     {
01059         vx = vdir.unit();                    // Ort in the direction of the reference particle
01060         G4ThreeVector vv= vx.orthogonal();   // Not normed orthogonal vector (!)
01061         vy = vv.unit();                      // First ort orthogonal to the direction
01062         vz = vx.cross(vy);                   // Second ort orthoganal to the direction
01063     }
01064     if(maxCost> 1.) maxCost= 1.;
01065     if(minCost<-1.) minCost=-1.;
01066     if(maxCost<-1.) maxCost=-1.;
01067     if(minCost> 1.) minCost= 1.;
01068     if(minCost> maxCost) minCost=maxCost;
01069     if(std::fabs(iM-fM-sM)<.00000001)
01070     {
01071         G4double fR=fM/iM;
01072         G4double sR=sM/iM;
01073         f4Mom=fR*theMomentum;
01074         s4Mom=sR*theMomentum;
01075         return true;
01076     }
01077     else if (iM+.001<fM+sM || iM==0.)
01078     {//@@ Later on make a quark content check for the decay
01079         G4cerr<<"***G4QH::RelDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl;
01080         return false;
01081     }
01082     G4double d2 = iM2-fM2-sM2;
01083     G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;    // Decay momentum(^2) in CMS of Quasmon
01084     if(p2<0.)
01085     {
01086         p2=0.;
01087     }
01088     G4double p  = std::sqrt(p2);
01089     G4double ct = maxCost;
01090     if(maxCost>minCost)
01091     {
01092         G4double dcost=maxCost-minCost;
01093         ct = minCost+dcost*G4UniformRand();
01094     }
01095     G4double phi= twopi*G4UniformRand();  // @@ Change 360.*deg to M_TWOPI (?)
01096     G4double ps=0.;
01097     if(std::fabs(ct)<1.) ps = p * std::sqrt(1.-ct*ct);
01098     else
01099     {
01100         if(ct>1.) ct=1.;
01101         if(ct<-1.) ct=-1.;
01102     }
01103     G4ThreeVector pVect=(ps*std::sin(phi))*vz+(ps*std::cos(phi))*vy+p*ct*vx;
01104     
01105     f4Mom.setVect(pVect);
01106     f4Mom.setE(std::sqrt(fM2+p2));
01107     s4Mom.setVect((-1)*pVect);
01108     s4Mom.setE(std::sqrt(sM2+p2));
01109     
01110     if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p="
01111         <<f4Mom.e()-f4Mom.rho()<<G4endl;
01112     f4Mom.boost(ltb);                        // Lor.Trans. of 1st hadron back to LS
01113     if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p="
01114         <<s4Mom.e()-s4Mom.rho()<<G4endl;
01115     s4Mom.boost(ltb);                        // Lor.Trans. of 2nd hadron back to LS
01116     return true;
01117 } // End of "RelDecayIn2"

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

Definition at line 779 of file G4QuasiElRatios.cc.

References G4Alpha::Alpha(), G4Deuteron::Deuteron(), FatalException, G4cerr, G4cout, G4endl, G4Exception(), G4ChipsNeutronElasticXS::GetChipsCrossSection(), G4ChipsProtonElasticXS::GetChipsCrossSection(), G4ChipsNeutronElasticXS::GetExchangeT(), G4ChipsProtonElasticXS::GetExchangeT(), G4ChipsNeutronElasticXS::GetHMaxT(), G4ChipsProtonElasticXS::GetHMaxT(), G4ParticleDefinition::GetPDGMass(), G4He3::He3(), G4Neutron::Neutron(), G4Proton::Proton(), RelDecayIn2(), and G4Triton::Triton().

Referenced by G4QuasiElasticChannel::Scatter().

00781 {    
00782     static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass();
00783     static const G4double mProt= G4Proton::Proton()->GetPDGMass();
00784     static const G4double mDeut= G4Deuteron::Deuteron()->GetPDGMass();
00785     static const G4double mTrit= G4Triton::Triton()->GetPDGMass();
00786     static const G4double mHel3= G4He3::He3()->GetPDGMass();
00787     static const G4double mAlph= G4Alpha::Alpha()->GetPDGMass();
00788     
00789     G4LorentzVector pr4M=p4M/megaelectronvolt;   // Convert 4-momenta in MeV (keep p4M)
00790     N4M/=megaelectronvolt;
00791     G4LorentzVector tot4M=N4M+p4M;
00792     G4double mT=mNeut;
00793     G4int Z=0;
00794     G4int N=1;
00795     if(NPDG==2212||NPDG==90001000)
00796     {
00797         mT=mProt;
00798         Z=1;
00799         N=0;
00800     }
00801     else if(NPDG==90001001)
00802     {
00803         mT=mDeut;
00804         Z=1;
00805         N=1;
00806     }
00807     else if(NPDG==90002001)
00808     {
00809         mT=mHel3;
00810         Z=2;
00811         N=1;
00812     }
00813     else if(NPDG==90001002)
00814     {
00815         mT=mTrit;
00816         Z=1;
00817         N=2;
00818     }
00819     else if(NPDG==90002002)
00820     {
00821         mT=mAlph;
00822         Z=2;
00823         N=2;
00824     }
00825     else if(NPDG!=2112&&NPDG!=90000001)
00826     {
00827         G4cout<<"Error:G4QuasiElRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
00828         G4Exception("G4QuasiElRatios::Scatter:","21",FatalException,"QEcomplain");
00829         //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
00830     }
00831     G4double mT2=mT*mT;
00832     G4double mP2=pr4M.m2();
00833     G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
00834     G4double E2=E*E;
00835     if(E<0. || E2<mP2)
00836     {
00837         return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
00838     }
00839     G4double P=std::sqrt(E2-mP2);                   // Momentum in pseudo laboratory system
00840     // @@ Temporary NN t-dependence for all hadrons
00841     if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QE::Scatter: pPDG="<<pPDG<<G4endl;
00842     G4int PDG=2212;                                                // *TMP* instead of pPDG
00843     if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112;               // *TMP* instead of pPDG
00844     if(!Z && N==1)                 // Change for Quasi-Elastic on neutron
00845     {
00846         Z=1;
00847         N=0;
00848         if     (PDG==2212) PDG=2112;
00849         else if(PDG==2112) PDG=2212;
00850     }
00851     G4double xSec=0.;                        // Prototype of Recalculated Cross Section *TMP*
00852     if(PDG==2212) xSec=PCSmanager->GetChipsCrossSection(P, Z, N, PDG); // P CrossSect *TMP*
00853     else          xSec=NCSmanager->GetChipsCrossSection(P, Z, N, PDG); // N CrossSect *TMP*
00854     // @@ check a possibility to separate p, n, or alpha (!)
00855     if(xSec <= 0.)                                    // The cross-section iz 0 -> Do Nothing
00856     {
00857         return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
00858     }
00859     G4double mint=0.;                        // Prototype of functional rand -t (MeV^2) *TMP*
00860     if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
00861     else          mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
00862     G4double maxt=0.;                                    // Prototype of max possible -t
00863     if(PDG==2212) maxt=PCSmanager->GetHMaxT();           // max possible -t
00864     else          maxt=NCSmanager->GetHMaxT();           // max possible -t
00865     G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS
00866     if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
00867     {
00868         if     (cost>1.)  cost=1.;
00869         else if(cost<-1.) cost=-1.;
00870         else
00871         {
00872             G4double tm=0.;
00873             if(PDG==2212) tm=PCSmanager->GetHMaxT();
00874             else          tm=NCSmanager->GetHMaxT();
00875             G4cerr<<"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<tm<<G4endl;
00876             return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
00877         }
00878     }
00879     G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT);      // 4mom of the recoil nucleon
00880     G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
00881     if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
00882     {
00883         G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl;
00884         //G4Exception("G4QFR::Scat:","009",FatalException,"Decay of ElasticComp");
00885         return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
00886     }
00887     return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
00888 } // End of Scatter


Public Member Functions
	~G4QuasiElRatios ()
std::pair< G4double, G4double >	GetRatios (G4double pIU, G4int prPDG, G4int tgZ, G4int tgN)
std::pair< G4double, G4double >	GetChExFactor (G4double pIU, G4int pPDG, G4int Z, G4int N)
std::pair< G4LorentzVector, G4LorentzVector >	Scatter (G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
std::pair< G4LorentzVector, G4LorentzVector >	ChExer (G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
std::pair< G4double, G4double >	GetElTot (G4double pIU, G4int hPDG, G4int Z, G4int N)
G4double	ChExElCoef (G4double p, G4int Z, G4int N, G4int pPDG)
std::pair< G4double, G4double >	GetElTotXS (G4double Mom, G4int PDG, G4bool F)
std::pair< G4double, G4double >	FetchElTot (G4double pGeV, G4int PDG, G4bool F)
G4bool	RelDecayIn2 (G4LorentzVector &theMomentum, G4LorentzVector &f4Mom, G4LorentzVector &s4Mom, G4LorentzVector &dir, G4double maxCost=1., G4double minCost=-1.)
Static Public Member Functions
static G4QuasiElRatios *	GetPointer ()
Protected Member Functions
	G4QuasiElRatios ()

G4QuasiElRatios Class Reference

Public Member Functions

Static Public Member Functions

Protected Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation