Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Protected Member Functions | List of all members
G4QuasiElRatios Class Reference

#include <G4QuasiElRatios.hh>

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 ()
 

Detailed Description

Definition at line 53 of file G4QuasiElRatios.hh.

Constructor & Destructor Documentation

◆ G4QuasiElRatios()

G4QuasiElRatios::G4QuasiElRatios ( )
protected

Definition at line 60 of file G4QuasiElRatios.cc.

61{
62
63 PCSmanager=(G4ChipsProtonElasticXS*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsProtonElasticXS::Default_Name());
64
65 NCSmanager=(G4ChipsNeutronElasticXS*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsNeutronElasticXS::Default_Name());
66}
G4ChipsNeutronElasticXS::Default_Name
static const char * Default_Name()
Definition: G4ChipsNeutronElasticXS.hh:56
G4ChipsProtonElasticXS::Default_Name
static const char * Default_Name()
Definition: G4ChipsProtonElasticXS.hh:54
G4CrossSectionDataSetRegistry::GetCrossSectionDataSet
G4VCrossSectionDataSet * GetCrossSectionDataSet(const G4String &name, G4bool warning=true)
Definition: G4CrossSectionDataSetRegistry.cc:137
G4CrossSectionDataSetRegistry::Instance
static G4CrossSectionDataSetRegistry * Instance()
Definition: G4CrossSectionDataSetRegistry.cc:74

◆ ~G4QuasiElRatios()

G4QuasiElRatios::~G4QuasiElRatios ( )

Definition at line 68 of file G4QuasiElRatios.cc.

69{
70 std::vector<G4double*>::iterator pos;
71 for(pos=vT.begin(); pos<vT.end(); pos++)
72 { delete [] *pos; }
73 vT.clear();
74 for(pos=vL.begin(); pos<vL.end(); pos++)
75 { delete [] *pos; }
76 vL.clear();
77
78 std::vector<std::pair<G4double,G4double>*>::iterator pos2;
79 for(pos2=vX.begin(); pos2<vX.end(); pos2++)
80 { delete [] *pos2; }
81 vX.clear();
82}

Member Function Documentation

◆ ChExElCoef()

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

Definition at line 1005 of file G4QuasiElRatios.cc.

1006{
1007 p/=MeV; // Converted from independent units
1008 G4double A=Z+N;
1009 if(A<1.5) return 0.;
1010 G4double C=0.;
1011 if (pPDG==2212) C=N/(A+Z);
1012 else if(pPDG==2112) C=Z/(A+N);
1013 else G4cout<<"*Warning*G4CohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl;
1014 C*=C; // Coherent processes squares the amplitude
1015 // @@ This is true only for nucleons: other projectiles must be treated differently
1016 G4double sp=std::sqrt(p);
1017 G4double p2=p*p;
1018 G4double p4=p2*p2;
1019 G4double dl1=std::log(p)-5.;
1020 G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013);
1021 G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p;
1022 G4double R=U/T;
1023 return C*R*R;
1024}
G4double
double G4double
Definition: G4Types.hh:64
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout

◆ ChExer()

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

Definition at line 893 of file G4QuasiElRatios.cc.

895{
896 static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass();
897 static const G4double mProt= G4Proton::Proton()->GetPDGMass();
898 G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV(keep p4M)
899 N4M/=megaelectronvolt;
900 G4LorentzVector tot4M=N4M+p4M;
901 G4int Z=0;
902 G4int N=1;
903 G4int sPDG=0; // PDG code of the scattered hadron
904 G4double mS=0.; // proto of mass of scattered hadron
905 G4double mT=mProt; // mass of the recoil nucleon
906 if(NPDG==2212)
907 {
908 mT=mNeut;
909 Z=1;
910 N=0;
911 if(pPDG==-211) sPDG=111; // pi+ -> pi0
912 else if(pPDG==-321)
913 {
914 sPDG=310; // K+ -> K0S
915 if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L
916 }
917 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321; // K0 -> K+ (?)
918 else if(pPDG==3112) sPDG=3212; // Sigma- -> Sigma0
919 else if(pPDG==3212) sPDG=3222; // Sigma0 -> Sigma+
920 else if(pPDG==3312) sPDG=3322; // Xi- -> Xi0
921 }
922 else if(NPDG==2112) // Default
923 {
924 if(pPDG==211) sPDG=111; // pi+ -> pi0
925 else if(pPDG==321)
926 {
927 sPDG=310; // K+ -> K0S
928 if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L
929 }
930 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321; // K0 -> K- (?)
931 else if(pPDG==3222) sPDG=3212; // Sigma+ -> Sigma0
932 else if(pPDG==3212) sPDG=3112; // Sigma0 -> Sigma-
933 else if(pPDG==3322) sPDG=3312; // Xi0 -> Xi-
934 }
935 else
936 {
937 G4cout<<"Error:G4QuasiElRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
938 G4Exception("G4QuasiElRatios::ChExer:","21",FatalException,"QE complain");
939 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
940 }
941 if(sPDG) mS=mNeut;
942 else
943 {
944 G4cout<<"Error:G4QuasiElRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl;
945 G4Exception("G4QuasiElRatios::ChExer:","21",FatalException,"QE complain");
946 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
947 }
948 G4double mT2=mT*mT;
949 G4double mS2=mS*mS;
950 G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);
951 G4double E2=E*E;
952 if(E<0. || E2<mS2)
953 {
954 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
955 }
956 G4double P=std::sqrt(E2-mS2); // Momentum in pseudo laboratory system
957 // @@ Temporary NN t-dependence for all hadrons
958 G4int PDG=2212; // *TMP* instead of pPDG
959 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG
960 if(!Z && N==1) // Change for Quasi-Elastic on neutron
961 {
962 Z=1;
963 N=0;
964 if (PDG==2212) PDG=2112;
965 else if(PDG==2112) PDG=2212;
966 }
967 G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP*
968 if(PDG==2212) xSec=PCSmanager->GetChipsCrossSection(P, Z, N, PDG); // P CrossSect *TMP*
969 else xSec=NCSmanager->GetChipsCrossSection(P, Z, N, PDG); // N CrossSect *TMP*
970 // @@ check a possibility to separate p, n, or alpha (!)
971 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
972 {
973 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
974 }
975 G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP*
976 if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
977 else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
978 G4double maxt=0.; // Prototype of max possible -t
979 if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t
980 else maxt=NCSmanager->GetHMaxT(); // max possible -t
981 G4double cost=1.-mint/maxt; // cos(theta) in CMS
982 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
983 {
984 if (cost>1.) cost=1.;
985 else if(cost<-1.) cost=-1.;
986 else
987 {
988 G4cerr<<"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<",t="<<mint<<",tm="<<maxt<<G4endl;
989 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
990 }
991 }
992 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon
993 pr4M=G4LorentzVector(0.,0.,0.,mS); // 4mom of the scattered hadron
994 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
995 if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
996 {
997 G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl;
998 //G4Exception("G4QFR::ChExer:","009",FatalException,"Decay of ElasticComp");
999 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
1000 }
1001 return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
1002} // End of ChExer
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4int
int G4int
Definition: G4Types.hh:66
G4cerr
G4DLLIMPORT std::ostream G4cerr
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4ChipsNeutronElasticXS::GetExchangeT
G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
Definition: G4ChipsNeutronElasticXS.cc:1853
G4ChipsNeutronElasticXS::GetChipsCrossSection
virtual G4double GetChipsCrossSection(G4double momentum, G4int Z, G4int N, G4int pdg)
Definition: G4ChipsNeutronElasticXS.cc:154
G4ChipsProtonElasticXS::GetExchangeT
G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
Definition: G4ChipsProtonElasticXS.cc:614
G4ChipsProtonElasticXS::GetChipsCrossSection
virtual G4double GetChipsCrossSection(G4double momentum, G4int Z, G4int N, G4int pdg)
Definition: G4ChipsProtonElasticXS.cc:160
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4QuasiElRatios::RelDecayIn2
G4bool RelDecayIn2(G4LorentzVector &theMomentum, G4LorentzVector &f4Mom, G4LorentzVector &s4Mom, G4LorentzVector &dir, G4double maxCost=1., G4double minCost=-1.)
Definition: G4QuasiElRatios.cc:1027
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

◆ FetchElTot()

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

Definition at line 588 of file G4QuasiElRatios.cc.

589{
590 static const G4int nlp=300; // Number of steps in the S(lnp) logTable(5% step)
591 static const G4int mlp=nlp+1; // Number of elements in the S(lnp) logTable
592 static const G4double lpi=-5.; // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
593 static const G4double lpa=10.; // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c)
594 static const G4double mi=std::exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c)
595 static const G4double ma=std::exp(lpa);// The max p of logTabEl(~ 22. TeV)
596 static const G4double dl=(lpa-lpi)/nlp;// Step of the logarithmic Table
597 static const G4double edl=std::exp(dl);// Multiplication step of the logarithmic Table
598 //static const G4double toler=.001; // Relative Tolarence defining "theSameMomentum"
599 static G4double lastP=0.; // The last momentum for which XS was calculated
600 static G4int lastH=0; // The last projPDG for which XS was calculated
601 static G4bool lastF=true; // The last nucleon for which XS was calculated
602 static std::pair<G4double,G4double> lastR(0.,0.); // The last result
603 // Local Associative Data Base:
604 static std::vector<G4int> vI; // Vector of index for which XS was calculated
605 static std::vector<G4double> vM; // Vector of rel max ln(p) initialized in LogTable
606 static std::vector<G4int> vK; // Vector of topBin number initialized in LogTable
607 // Last values of the Associative Data Base:
608 static G4int lastI=0; // The Last index for which XS was calculated
609 static G4double lastM=0.; // The Last rel max ln(p) initialized in LogTable
610 static G4int lastK=0; // The Last topBin number initialized in LogTable
611 static std::pair<G4double,G4double>* lastX=0; // The Last ETPointers to LogTable in heap
612 // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei
613 G4int nDB=vI.size(); // A number of hadrons already initialized in AMDB
614 if(nDB && lastH==PDG && lastF==F && p>0. && p==lastP) return lastR;// VI don't use toler.
615 // if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
616 lastH=PDG;
617 lastF=F;
618 G4int ind=-1; // Prototipe of the index of the PDG/F combination
619 // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p),
620 // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp)
621 G4bool kfl=true; // Flag of K0/aK0 oscillation
622 G4bool kf=false;
623 if(PDG==130||PDG==310)
624 {
625 kf=true;
626 if(G4UniformRand()>.5) kfl=false;
627 }
628 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
629 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
630 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
631 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
632 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
633 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
634 else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda
635 else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
636 else {
637 G4cout<<"*Error*G4QuasiElRatios::FetchElTot: PDG="<<PDG
638 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
639 G4Exception("G4QuasiELRatio::FetchElTot:","22",FatalException,"QECrash");
640 }
641 if(nDB && lastI==ind && p>0. && p==lastP) return lastR; // VI do not use toler
642 // if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR;
643 if(p<=mi || p>=ma) return CalcElTot(p,ind); // @@ Slow calculation ! (Warning?)
644 G4bool found=false;
645 G4int i=-1;
646 if(nDB) for (i=0; i<nDB; i++) if(ind==vI[i]) // Sirch for this index in AMDB
647 {
648 found=true; // The index is found
649 break;
650 }
651 G4double lp=std::log(p);
652 if(!nDB || !found) // Create new line in the AMDB
653 {
654 lastX = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot
655 lastI = ind; // Remember the initialized inex
656 lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTaB
657 if(lastK>nlp)
658 {
659 lastK=nlp;
660 lastM=lpa-lpi;
661 }
662 else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab
663 G4double pv=mi;
664 for(G4int j=0; j<=lastK; j++) // Calculate LogTab values
665 {
666 lastX[j]=CalcElTot(pv,ind);
667 if(j!=lastK) pv*=edl;
668 }
669 i++; // Make a new record to AMDB and position on it
670 vI.push_back(lastI);
671 vM.push_back(lastM);
672 vK.push_back(lastK);
673 vX.push_back(lastX);
674 }
675 else // The A value was found in AMDB
676 {
677 lastI=vI[i];
678 lastM=vM[i];
679 lastK=vK[i];
680 lastX=vX[i];
681 G4int nextK=lastK+1;
682 G4double lpM=lastM+lpi;
683 if(lp>lpM && lastK<nlp) // LogTab must be updated
684 {
685 lastK = static_cast<int>((lp-lpi)/dl)+1; // MaxBin to be initialized in LogTab
686 if(lastK>nlp)
687 {
688 lastK=nlp;
689 lastM=lpa-lpi;
690 }
691 else lastM = lastK*dl; // Calculate max initialized ln(p)-lpi for LogTab
692 G4double pv=std::exp(lpM); // momentum of the last calculated beam
693 for(G4int j=nextK; j<=lastK; j++)// Calculate LogTab values
694 {
695 pv*=edl;
696 lastX[j]=CalcElTot(pv,ind);
697 }
698 } // End of LogTab update
699 if(lastK>=nextK) // The AMDB was apdated
700 {
701 vM[i]=lastM;
702 vK[i]=lastK;
703 }
704 }
705 // Now one can use tabeles to calculate the value
706 G4double dlp=lp-lpi; // Shifted log(p) value
707 G4int n=static_cast<int>(dlp/dl); // Low edge number of the bin
708 G4double d=dlp-n*dl; // Log shift
709 G4double e=lastX[n].first; // E-Base
710 lastR.first=e+d*(lastX[n+1].first-e)/dl; // E-Result
711 if(lastR.first<0.) lastR.first = 0.;
712 G4double t=lastX[n].second; // T-Base
713 lastR.second=t+d*(lastX[n+1].second-t)/dl; // T-Result
714 if(lastR.second<0.) lastR.second= 0.;
715 if(lastR.first>lastR.second) lastR.first = lastR.second;
716 return lastR;
717} // End of FetchElTot
G4bool
bool G4bool
Definition: G4Types.hh:67

Referenced by GetChExFactor(), and GetElTot().

◆ GetChExFactor()

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

Definition at line 736 of file G4QuasiElRatios.cc.

738{
739 G4double pGeV=pIU/gigaelectronvolt;
740 G4double resP=0.;
741 G4double resN=0.;
742 if(Z<1 && N<1)
743 {
744 G4cout<<"-Warning-G4QuasiElRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
745 return std::make_pair(resP,resN);
746 }
747 G4double A=Z+N;
748 G4double pf=0.; // Possibility to interact with a proton
749 G4double nf=0.; // Possibility to interact with a neutron
750 if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+N);
751 else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z);
752 else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310)
753 {
754 G4double dA=A+A;
755 pf=Z/(dA+N+N);
756 nf=N/(dA+Z+Z);
757 }
758 G4double mult=1.; // Factor of increasing multiplicity ( ? @@)
759 if(pGeV>.5)
760 {
761 mult=1./(1.+std::log(pGeV+pGeV))/pGeV;
762 if(mult>1.) mult=1.;
763 }
764 if(pf)
765 {
766 std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
767 resP=pf*(hp.second/hp.first-1.)*mult;
768 }
769 if(nf)
770 {
771 std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
772 resN=nf*(hn.second/hn.first-1.)*mult;
773 }
774 return std::make_pair(resP,resN);
775} // End of GetChExFactor
G4QuasiElRatios::FetchElTot
std::pair< G4double, G4double > FetchElTot(G4double pGeV, G4int PDG, G4bool F)
Definition: G4QuasiElRatios.cc:588

◆ GetElTot()

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

Definition at line 720 of file G4QuasiElRatios.cc.

722{
723 G4double pGeV=pIU/gigaelectronvolt;
724 if(Z<1 && N<1)
725 {
726 G4cout<<"-Warning-G4QuasiElRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
727 return std::make_pair(0.,0.);
728 }
729 std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true);
730 std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false);
731 G4double A=(Z+N)/millibarn; // To make the result in independent units(IU)
732 return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
733} // End of GetElTot

Referenced by GetRatios().

◆ GetElTotXS()

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

Definition at line 561 of file G4QuasiElRatios.cc.

562{
563 G4int ind=0; // Prototype of the reaction index
564 G4bool kfl=true; // Flag of K0/aK0 oscillation
565 G4bool kf=false;
566 if(PDG==130||PDG==310)
567 {
568 kf=true;
569 if(G4UniformRand()>.5) kfl=false;
570 }
571 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
572 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
573 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
574 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
575 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
576 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
577 else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda
578 else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
579 else {
580 G4cout<<"*Error*G4QuasiElRatios::CalcElTotXS: PDG="<<PDG
581 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
582 G4Exception("G4QuasiElRatio::CalcElTotXS:","22",FatalException,"QEcrash");
583 }
584 return CalcElTot(p,ind);
585}

◆ GetPointer()

G4QuasiElRatios * G4QuasiElRatios::GetPointer ( )
static

Definition at line 85 of file G4QuasiElRatios.cc.

86{
87 static G4QuasiElRatios theRatios; // *** Static body of the QEl Cross Section ***
88 return &theRatios;
89}

◆ GetRatios()

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

Definition at line 92 of file G4QuasiElRatios.cc.

94{
95 G4double R=0.;
96 G4double QF2In=1.; // Prototype of QuasiFree/Inel ratio for hN_tot
97 G4int tgA=tgZ+tgN;
98 if(tgA<2) return std::make_pair(QF2In,R); // No quasi-elastic on the only nucleon
99 std::pair<G4double,G4double> ElTot=GetElTot(pIU, pPDG, tgZ, tgN); // mean hN El&Tot(IU)
100 //if( ( (pPDG>999 && pIU<227.) || pIU<27.) && tgA>1) R=1.; // @@ TMP to accelerate @lowE
101 if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.; // To accelerate @lowE
102 else if(ElTot.second>0.)
103 {
104 R=ElTot.first/ElTot.second; // El/Total ratio (does not depend on units
105 QF2In=GetQF2IN_Ratio(ElTot.second/millibarn, tgZ+tgN); // QuasiFree/Inelastic ratio
106 }
107 return std::make_pair(QF2In,R);
108}
G4QuasiElRatios::GetElTot
std::pair< G4double, G4double > GetElTot(G4double pIU, G4int hPDG, G4int Z, G4int N)
Definition: G4QuasiElRatios.cc:720

Referenced by G4QuasiElasticChannel::GetFraction().

◆ RelDecayIn2()

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.

1029{
1030 G4double fM2 = f4Mom.m2();
1031 G4double fM = std::sqrt(fM2); // Mass of the 1st Hadron
1032 G4double sM2 = s4Mom.m2();
1033 G4double sM = std::sqrt(sM2); // Mass of the 2nd Hadron
1034 G4double iM2 = theMomentum.m2();
1035 G4double iM = std::sqrt(iM2); // Mass of the decaying hadron
1036 G4double vP = theMomentum.rho(); // Momentum of the decaying hadron
1037 G4double dE = theMomentum.e(); // Energy of the decaying hadron
1038 if(dE<vP)
1039 {
1040 G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
1041 G4double accuracy=.000001*vP;
1042 G4double emodif=std::fabs(dE-vP);
1043 //if(emodif<accuracy)
1044 //{
1045 G4cerr<<"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
1046 theMomentum.setE(vP+emodif+.01*accuracy);
1047 //}
1048 }
1049 G4ThreeVector ltb = theMomentum.boostVector();// Boost vector for backward Lorentz Trans.
1050 G4ThreeVector ltf = -ltb; // Boost vector for forward Lorentz Trans.
1051 G4LorentzVector cdir = dir; // A copy to make a transformation to CMS
1052 cdir.boost(ltf); // Direction transpormed to CMS of the Momentum
1053 G4ThreeVector vdir = cdir.vect(); // 3-Vector of the direction-particle
1054 G4ThreeVector vx(0.,0.,1.); // Ort in the direction of the reference particle
1055 G4ThreeVector vy(0.,1.,0.); // First ort orthogonal to the direction
1056 G4ThreeVector vz(1.,0.,0.); // Second ort orthoganal to the direction
1057 if(vdir.mag2() > 0.) // the refference particle isn't at rest in CMS
1058 {
1059 vx = vdir.unit(); // Ort in the direction of the reference particle
1060 G4ThreeVector vv= vx.orthogonal(); // Not normed orthogonal vector (!)
1061 vy = vv.unit(); // First ort orthogonal to the direction
1062 vz = vx.cross(vy); // Second ort orthoganal to the direction
1063 }
1064 if(maxCost> 1.) maxCost= 1.;
1065 if(minCost<-1.) minCost=-1.;
1066 if(maxCost<-1.) maxCost=-1.;
1067 if(minCost> 1.) minCost= 1.;
1068 if(minCost> maxCost) minCost=maxCost;
1069 if(std::fabs(iM-fM-sM)<.00000001)
1070 {
1071 G4double fR=fM/iM;
1072 G4double sR=sM/iM;
1073 f4Mom=fR*theMomentum;
1074 s4Mom=sR*theMomentum;
1075 return true;
1076 }
1077 else if (iM+.001<fM+sM || iM==0.)
1078 {//@@ Later on make a quark content check for the decay
1079 G4cerr<<"***G4QH::RelDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl;
1080 return false;
1081 }
1082 G4double d2 = iM2-fM2-sM2;
1083 G4double p2 = (d2*d2/4.-fM2*sM2)/iM2; // Decay momentum(^2) in CMS of Quasmon
1084 if(p2<0.)
1085 {
1086 p2=0.;
1087 }
1088 G4double p = std::sqrt(p2);
1089 G4double ct = maxCost;
1090 if(maxCost>minCost)
1091 {
1092 G4double dcost=maxCost-minCost;
1093 ct = minCost+dcost*G4UniformRand();
1094 }
1095 G4double phi= twopi*G4UniformRand(); // @@ Change 360.*deg to M_TWOPI (?)
1096 G4double ps=0.;
1097 if(std::fabs(ct)<1.) ps = p * std::sqrt(1.-ct*ct);
1098 else
1099 {
1100 if(ct>1.) ct=1.;
1101 if(ct<-1.) ct=-1.;
1102 }
1103 G4ThreeVector pVect=(ps*std::sin(phi))*vz+(ps*std::cos(phi))*vy+p*ct*vx;
1104
1105 f4Mom.setVect(pVect);
1106 f4Mom.setE(std::sqrt(fM2+p2));
1107 s4Mom.setVect((-1)*pVect);
1108 s4Mom.setE(std::sqrt(sM2+p2));
1109
1110 if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p="
1111 <<f4Mom.e()-f4Mom.rho()<<G4endl;
1112 f4Mom.boost(ltb); // Lor.Trans. of 1st hadron back to LS
1113 if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p="
1114 <<s4Mom.e()-s4Mom.rho()<<G4endl;
1115 s4Mom.boost(ltb); // Lor.Trans. of 2nd hadron back to LS
1116 return true;
1117} // End of "RelDecayIn2"
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::orthogonal
Hep3Vector orthogonal() const
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:177
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:59
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
CLHEP::HepLorentzVector::setVect
void setVect(const Hep3Vector &)

Referenced by ChExer(), and Scatter().

◆ Scatter()

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

Definition at line 779 of file G4QuasiElRatios.cc.

781{
782 static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass();
783 static const G4double mProt= G4Proton::Proton()->GetPDGMass();
784 static const G4double mDeut= G4Deuteron::Deuteron()->GetPDGMass();
785 static const G4double mTrit= G4Triton::Triton()->GetPDGMass();
786 static const G4double mHel3= G4He3::He3()->GetPDGMass();
787 static const G4double mAlph= G4Alpha::Alpha()->GetPDGMass();
788
789 G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV (keep p4M)
790 N4M/=megaelectronvolt;
791 G4LorentzVector tot4M=N4M+p4M;
792 G4double mT=mNeut;
793 G4int Z=0;
794 G4int N=1;
795 if(NPDG==2212||NPDG==90001000)
796 {
797 mT=mProt;
798 Z=1;
799 N=0;
800 }
801 else if(NPDG==90001001)
802 {
803 mT=mDeut;
804 Z=1;
805 N=1;
806 }
807 else if(NPDG==90002001)
808 {
809 mT=mHel3;
810 Z=2;
811 N=1;
812 }
813 else if(NPDG==90001002)
814 {
815 mT=mTrit;
816 Z=1;
817 N=2;
818 }
819 else if(NPDG==90002002)
820 {
821 mT=mAlph;
822 Z=2;
823 N=2;
824 }
825 else if(NPDG!=2112&&NPDG!=90000001)
826 {
827 G4cout<<"Error:G4QuasiElRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
828 G4Exception("G4QuasiElRatios::Scatter:","21",FatalException,"QEcomplain");
829 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
830 }
831 G4double mT2=mT*mT;
832 G4double mP2=pr4M.m2();
833 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
834 G4double E2=E*E;
835 if(E<0. || E2<mP2)
836 {
837 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
838 }
839 G4double P=std::sqrt(E2-mP2); // Momentum in pseudo laboratory system
840 // @@ Temporary NN t-dependence for all hadrons
841 if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QE::Scatter: pPDG="<<pPDG<<G4endl;
842 G4int PDG=2212; // *TMP* instead of pPDG
843 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG
844 if(!Z && N==1) // Change for Quasi-Elastic on neutron
845 {
846 Z=1;
847 N=0;
848 if (PDG==2212) PDG=2112;
849 else if(PDG==2112) PDG=2212;
850 }
851 G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP*
852 if(PDG==2212) xSec=PCSmanager->GetChipsCrossSection(P, Z, N, PDG); // P CrossSect *TMP*
853 else xSec=NCSmanager->GetChipsCrossSection(P, Z, N, PDG); // N CrossSect *TMP*
854 // @@ check a possibility to separate p, n, or alpha (!)
855 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
856 {
857 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
858 }
859 G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP*
860 if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
861 else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
862 G4double maxt=0.; // Prototype of max possible -t
863 if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t
864 else maxt=NCSmanager->GetHMaxT(); // max possible -t
865 G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS
866 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
867 {
868 if (cost>1.) cost=1.;
869 else if(cost<-1.) cost=-1.;
870 else
871 {
872 G4double tm=0.;
873 if(PDG==2212) tm=PCSmanager->GetHMaxT();
874 else tm=NCSmanager->GetHMaxT();
875 G4cerr<<"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<tm<<G4endl;
876 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
877 }
878 }
879 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon
880 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
881 if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost))
882 {
883 G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl;
884 //G4Exception("G4QFR::Scat:","009",FatalException,"Decay of ElasticComp");
885 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
886 }
887 return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
888} // End of Scatter
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:94
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95

Referenced by G4QuasiElasticChannel::Scatter().

The documentation for this class was generated from the following files: