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
G4QuasiFreeRatios Class Reference

#include <G4QuasiFreeRatios.hh>

Public Member Functions

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

Static Public Member Functions

static G4QuasiFreeRatios * GetPointer ()
 

Protected Member Functions

 G4QuasiFreeRatios ()
 

Detailed Description

Definition at line 52 of file G4QuasiFreeRatios.hh.

Constructor & Destructor Documentation

◆ G4QuasiFreeRatios()

G4QuasiFreeRatios::G4QuasiFreeRatios ( )
protected

Definition at line 51 of file G4QuasiFreeRatios.cc.

52{
53#ifdef pdebug
54 G4cout<<"***^^^*** G4QuasiFreeRatios singletone is created ***^^^***"<<G4endl;
55#endif
56 QFreeScat=G4QFreeScattering::GetPointer();
57}
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4QFreeScattering::GetPointer
static G4QFreeScattering * GetPointer()
Definition: G4QFreeScattering.cc:64

◆ ~G4QuasiFreeRatios()

G4QuasiFreeRatios::~G4QuasiFreeRatios ( )

Definition at line 59 of file G4QuasiFreeRatios.cc.

60{
61 std::vector<G4double*>::iterator pos;
62 for(pos=vT.begin(); pos<vT.end(); ++pos) delete [] *pos;
63 vT.clear();
64 for(pos=vL.begin(); pos<vL.end(); ++pos) delete [] *pos;
65 vL.clear();
66}

Member Function Documentation

◆ ChExElCoef()

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

Definition at line 734 of file G4QuasiFreeRatios.cc.

735{
736 p/=MeV; // Converted from independent units
737 G4double A=Z+N;
738 if(A<1.5) return 0.;
739 G4double C=0.;
740 if (pPDG==2212) C=N/(A+Z);
741 else if(pPDG==2112) C=Z/(A+N);
742 else G4cout<<"*Warning*G4QCohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl;
743 C*=C; // Coherent processes squares the amplitude
744 // @@ This is true only for nucleons: other projectiles must be treated differently
745 G4double sp=std::sqrt(p);
746 G4double p2=p*p;
747 G4double p4=p2*p2;
748 G4double dl1=std::log(p)-5.;
749 G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013);
750 G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p;
751 G4double R=U/T;
752 return C*R*R;
753}
G4double
double G4double
Definition: G4Types.hh:64

Referenced by G4QInelastic::PostStepDoIt().

◆ ChExer()

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

Definition at line 590 of file G4QuasiFreeRatios.cc.

592{
593 static const G4double mNeut= G4QPDGCode(2112).GetMass();
594 static const G4double mProt= G4QPDGCode(2212).GetMass();
595 G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV(keep p4M)
596 N4M/=megaelectronvolt;
597 G4LorentzVector tot4M=N4M+p4M;
598 G4int Z=0;
599 G4int N=1;
600 G4int sPDG=0; // PDG code of the scattered hadron
601 G4double mS=0.; // proto of mass of scattered hadron
602 G4double mT=mProt; // mass of the recoil nucleon
603 G4cout<<"-Warning-G4QFR::ChEx: Impossible for Omega-"<<G4endl; // Omega-
604 if(NPDG==2212)
605 {
606 mT=mNeut;
607 Z=1;
608 N=0;
609 if(pPDG==-211) sPDG=111; // pi+ -> pi0
610 else if(pPDG==-321)
611 {
612 sPDG=310; // K+ -> K0S
613 if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L
614 }
615 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321; // K0 -> K+ (?)
616 else if(pPDG==3112) sPDG=3212; // Sigma- -> Sigma0
617 else if(pPDG==3212) sPDG=3222; // Sigma0 -> Sigma+
618 else if(pPDG==3312) sPDG=3322; // Xi- -> Xi0
619 }
620 else if(NPDG==2112) // Default
621 {
622 if(pPDG==211) sPDG=111; // pi+ -> pi0
623 else if(pPDG==321)
624 {
625 sPDG=310; // K+ -> K0S
626 if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L
627 }
628 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321; // K0 -> K- (?)
629 else if(pPDG==3222) sPDG=3212; // Sigma+ -> Sigma0
630 else if(pPDG==3212) sPDG=3112; // Sigma0 -> Sigma-
631 else if(pPDG==3322) sPDG=3312; // Xi0 -> Xi-
632 }
633 else
634 {
635 G4cout<<"Error:G4QuasiFreeRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
636 G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain");
637 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
638 }
639 if(sPDG) mS=mNeut;
640 else
641 {
642 G4cout<<"Error:G4QuasiFreeRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl;
643 G4Exception("G4QuasiFreeRatios::ChExer:","21",FatalException,"CHIPS complain");
644 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
645 }
646 G4double mT2=mT*mT;
647 G4double mS2=mS*mS;
648 G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT);
649 G4double E2=E*E;
650 if(E<0. || E2<mS2)
651 {
652#ifdef pdebug
653 G4cerr<<"-Warning-G4QFR::ChEx:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mS2<<G4endl;
654#endif
655 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
656 }
657 G4double P=std::sqrt(E2-mS2); // Momentum in pseudo laboratory system
658 G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
659 G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer();
660#ifdef debug
661 G4cout<<"G4QFR::ChExer: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl;
662#endif
663 // @@ Temporary NN t-dependence for all hadrons
664 G4int PDG=2212; // *TMP* instead of pPDG
665 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG
666 if(!Z && N==1) // Change for Quasi-Elastic on neutron
667 {
668 Z=1;
669 N=0;
670 if (PDG==2212) PDG=2112;
671 else if(PDG==2112) PDG=2212;
672 }
673 G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP*
674 if(PDG==2212) xSec=PCSmanager->GetCrossSection(false, P, Z, N, PDG); // P CrossSect *TMP*
675 else xSec=NCSmanager->GetCrossSection(false, P, Z, N, PDG); // N CrossSect *TMP*
676#ifdef debug
677 G4cout<<"G4QuasiFreeRat::ChExer: pPDG="<<pPDG<<",P="<<P<<", CS="<<xSec/millibarn<<G4endl;
678#endif
679#ifdef nandebug
680 if(xSec>0. || xSec<0. || xSec==0);
681 else G4cout<<"*Warning*G4QuasiFreeRatios::ChExer: xSec="<<xSec/millibarn<<G4endl;
682#endif
683 // @@ check a possibility to separate p, n, or alpha (!)
684 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
685 {
686#ifdef pdebug
687 G4cerr<<"-Warning-G4QFR::ChEx:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl;
688#endif
689 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
690 }
691 G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP*
692 if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
693 else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
694#ifdef pdebug
695 G4cout<<"G4QFR::ChEx:PDG="<<pPDG<<", P="<<P<<", CS="<<xSec<<", -t="<<mint<<G4endl;
696#endif
697#ifdef nandebug
698 if(mint>-.0000001);
699 else G4cout<<"*Warning*G4QFR::ChExer: -t="<<mint<<G4endl;
700#endif
701 G4double maxt=0.; // Prototype of max possible -t
702 if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t
703 else maxt=NCSmanager->GetHMaxT(); // max possible -t
704 G4double cost=1.-mint/maxt; // cos(theta) in CMS
705#ifdef pdebug
706 G4cout<<"G4QuasiFfreeRatio::ChExer: -t="<<mint<<", maxt="<<maxt<<", cost="<<cost<<G4endl;
707#endif
708 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
709 {
710 if (cost>1.) cost=1.;
711 else if(cost<-1.) cost=-1.;
712 else
713 {
714 G4cerr<<"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<",t="<<mint<<",tm="<<maxt<<G4endl;
715 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
716 }
717 }
718 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon
719 pr4M=G4LorentzVector(0.,0.,0.,mS); // 4mom of the scattered hadron
720 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
721 if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost))
722 {
723 G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl;
724 //G4Exception("G4QFR::ChExer:","009",FatalException,"Decay of ElasticComp");
725 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
726 }
727#ifdef debug
728 G4cout<<"G4QFR::ChEx:p4M="<<p4M<<"+r4M="<<reco4M<<"="<<p4M+reco4M<<"="<<tot4M<<G4endl;
729#endif
730 return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
731} // 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
G4QPDGCode::GetMass
G4double GetMass()
Definition: G4QPDGCode.cc:693
G4QProtonElasticCrossSection::GetPointer
static G4VQCrossSection * GetPointer()
Definition: G4QProtonElasticCrossSection.cc:143
G4VQCrossSection::GetExchangeT
virtual G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
Definition: G4VQCrossSection.cc:70
G4VQCrossSection::GetCrossSection
virtual G4double GetCrossSection(G4bool, G4double, G4int, G4int, G4int pPDG=0)
Definition: G4VQCrossSection.hh:99
G4VQCrossSection::GetHMaxT
virtual G4double GetHMaxT()
Definition: G4VQCrossSection.cc:72
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

◆ GetChExFactor()

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

Definition at line 389 of file G4QuasiFreeRatios.cc.

391{
392 G4double pGeV=pIU/gigaelectronvolt;
393 G4double resP=0.;
394 G4double resN=0.;
395 if(Z<1 && N<1)
396 {
397 G4cout<<"-Warning-G4QuasiFreeRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
398 return std::make_pair(resP,resN);
399 }
400 G4double A=Z+N;
401 G4double pf=0.; // Possibility to interact with a proton
402 G4double nf=0.; // Possibility to interact with a neutron
403 if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312||hPDG==3334) pf=Z/(A+N);
404 else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z);
405 else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310)
406 {
407 G4double dA=A+A;
408 pf=Z/(dA+N+N);
409 nf=N/(dA+Z+Z);
410 }
411 G4double mult=1.; // Factor of increasing multiplicity ( ? @@)
412 if(pGeV>.5)
413 {
414 mult=1./(1.+std::log(pGeV+pGeV))/pGeV;
415 if(mult>1.) mult=1.;
416 }
417 if(pf)
418 {
419 std::pair<G4double,G4double> hp=QFreeScat->FetchElTot(pGeV, hPDG, true);
420 resP=pf*(hp.second/hp.first-1.)*mult;
421 }
422 if(nf)
423 {
424 std::pair<G4double,G4double> hn=QFreeScat->FetchElTot(pGeV, hPDG, false);
425 resN=nf*(hn.second/hn.first-1.)*mult;
426 }
427 return std::make_pair(resP,resN);
428} // End of GetChExFactor
G4QFreeScattering::FetchElTot
std::pair< G4double, G4double > FetchElTot(G4double pGeV, G4int PDG, G4bool F)
Definition: G4QFreeScattering.cc:373

◆ GetElTot()

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

Definition at line 366 of file G4QuasiFreeRatios.cc.

368{
369 G4double pGeV=pIU/gigaelectronvolt;
370#ifdef pdebug
371 G4cout<<"-->G4QuasiFreeR::GetElTot: P="<<pIU<<",pPDG="<<hPDG<<",Z="<<Z<<",N="<<N<<G4endl;
372#endif
373 if(Z<1 && N<1)
374 {
375 G4cout<<"-Warning-G4QuasiFreeRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl;
376 return std::make_pair(0.,0.);
377 }
378 std::pair<G4double,G4double> hp=QFreeScat->FetchElTot(pGeV, hPDG, true);
379 std::pair<G4double,G4double> hn=QFreeScat->FetchElTot(pGeV, hPDG, false);
380#ifdef pdebug
381 G4cout<<"-OUT->G4QFRat::GetElTot: hp("<<hp.first<<","<<hp.second<<"), hn("<<hn.first<<","
382 <<hn.second<<")"<<G4endl;
383#endif
384 G4double A=(Z+N)/millibarn; // To make the result in independent units(IU)
385 return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A);
386} // End of GetElTot

Referenced by GetRatios().

◆ GetElTotXS()

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

Definition at line 339 of file G4QuasiFreeRatios.cc.

340{
341 G4int ind=0; // Prototype of the reaction index
342 G4bool kfl=true; // Flag of K0/aK0 oscillation
343 G4bool kf=false;
344 if(PDG==130||PDG==310)
345 {
346 kf=true;
347 if(G4UniformRand()>.5) kfl=false;
348 }
349 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn
350 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn
351 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn
352 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn
353 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) ) ind=4; // KmN/K0N
354 else if ( PDG == 321 || PDG == 311 || (kf && kfl) ) ind=5; // KpN/aK0N
355 else if ( PDG > 3000 && PDG < 3335) ind=6; // @@ for all hyperons - take Lambda
356 else if ( PDG > -3335 && PDG < -2000) ind=7; // @@ for all anti-baryons (anti-p/anti-n)
357 else {
358 G4cout<<"*Error*G4QuasiFreeRatios::CalcElTotXS: PDG="<<PDG
359 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
360 G4Exception("G4QuasiFreeRatio::CalcElTotXS:","22",FatalException,"CHIPScrash");
361 }
362 return QFreeScat->CalcElTot(p,ind);
363}
G4bool
bool G4bool
Definition: G4Types.hh:67
G4QFreeScattering::CalcElTot
std::pair< G4double, G4double > CalcElTot(G4double pGeV, G4int Index)
Definition: G4QFreeScattering.cc:71

◆ GetPointer()

G4QuasiFreeRatios * G4QuasiFreeRatios::GetPointer ( )
static

Definition at line 69 of file G4QuasiFreeRatios.cc.

70{
71 static G4QuasiFreeRatios theRatios; // *** Static body of the QEl Cross Section ***
72 return &theRatios;
73}

Referenced by G4QFragmentation::G4QFragmentation(), and G4QInelastic::PostStepDoIt().

◆ GetRatios()

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

Definition at line 76 of file G4QuasiFreeRatios.cc.

78{
79#ifdef pdebug
80 G4cout<<">>>IN>>>G4QFRat::GetQF:P="<<pIU<<",pPDG="<<pPDG<<",Z="<<tgZ<<",N="<<tgN<<G4endl;
81#endif
82 G4double R=0.;
83 G4double QF2In=1.; // Prototype of QuasiFree/Inel ratio for hN_tot
84 G4int tgA=tgZ+tgN;
85 if(tgA<2) return std::make_pair(QF2In,R); // No quasi-elastic on the only nucleon
86 std::pair<G4double,G4double> ElTot=GetElTot(pIU, pPDG, tgZ, tgN); // mean (IU)
87 //if( ( (pPDG>999 && pIU<227.) || pIU<27.) && tgA>1) R=1.; // @@ TMP to accelerate @lowE
88 if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.; // To accelerate @lowE
89 else if(ElTot.second>0.)
90 {
91 R=ElTot.first/ElTot.second; // El/Total ratio (does not depend on units
92 QF2In=GetQF2IN_Ratio(ElTot.second/millibarn, tgZ+tgN); // QuasiFree/Inelastic ratio
93 }
94#ifdef pdebug
95 G4cout<<">>>OUT>>>G4QuasiFreeRatio::GetQF2IN_Ratio: QF2In="<<QF2In<<", R="<<R<<G4endl;
96#endif
97 return std::make_pair(QF2In,R);
98}
G4QuasiFreeRatios::GetElTot
std::pair< G4double, G4double > GetElTot(G4double pIU, G4int hPDG, G4int Z, G4int N)
Definition: G4QuasiFreeRatios.cc:366

Referenced by G4QFragmentation::G4QFragmentation(), and G4QInelastic::PostStepDoIt().

◆ Scatter()

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

Definition at line 432 of file G4QuasiFreeRatios.cc.

434{
435 static const G4double mNeut= G4QPDGCode(2112).GetMass();
436 static const G4double mProt= G4QPDGCode(2212).GetMass();
437 static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron
438 static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium
439 static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3
440 static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);// Mass of alpha
441 static const G4double two3d= 2./3.;
442 static const G4double two3d2= std::pow(2.,two3d); // The slope reductions for fragments
443 static const G4double two3d3= std::pow(3.,two3d);
444 static const G4double two3d4= std::pow(4.,two3d);
445 G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV (keep p4M)
446 N4M/=megaelectronvolt;
447 G4LorentzVector tot4M=N4M+p4M;
448#ifdef ppdebug
449 G4cerr<<"->G4QFR::Scat:p4M="<<pr4M<<",N4M="<<N4M<<",t4M="<<tot4M<<",NPDG="<<NPDG<<G4endl;
450#endif
451 G4double mT=mNeut;
452 G4int Z=0;
453 G4int N=1;
454 if(NPDG==2212||NPDG==90001000)
455 {
456 mT=mProt;
457 Z=1;
458 N=0;
459 }
460 else if(NPDG==90001001)
461 {
462 mT=mDeut;
463 Z=1;
464 N=1;
465 }
466 else if(NPDG==90002001)
467 {
468 mT=mHel3;
469 Z=2;
470 N=1;
471 }
472 else if(NPDG==90001002)
473 {
474 mT=mTrit;
475 Z=1;
476 N=2;
477 }
478 else if(NPDG==90002002)
479 {
480 mT=mAlph;
481 Z=2;
482 N=2;
483 }
484 else if(NPDG!=2112&&NPDG!=90000001)
485 {
486 G4cout<<"Error:G4QuasiFreeRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl;
487 G4Exception("G4QuasiFreeRatios::Scatter:","21",FatalException,"CHIPScomplain");
488 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception
489 }
490 G4double mT2=mT*mT;
491 G4double mP2=pr4M.m2();
492 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT);
493#ifdef pdebug
494 G4cerr<<"G4QFR::Scat:qM="<<mT<<",qM2="<<mT2<<",pM2="<<mP2<<",totM2="<<tot4M.m2()<<G4endl;
495#endif
496 G4double E2=E*E;
497 if(E<0. || E2<mP2)
498 {
499#ifdef ppdebug
500 G4cerr<<"-Warning-G4QFR::Scat:*Negative Energy*E="<<E<<",E2="<<E2<<"<M2="<<mP2<<G4endl;
501#endif
502 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
503 }
504 G4double P=std::sqrt(E2-mP2); // Momentum in pseudo laboratory system
505 G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
506 G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer();
507#ifdef ppdebug
508 G4cout<<"G4QFR::Scatter: Before XS, P="<<P<<", Z="<<Z<<", N="<<N<<", PDG="<<pPDG<<G4endl;
509#endif
510 // @@ Temporary NN t-dependence for all hadrons
511 if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QuasiFree::Scatter: pPDG="<<pPDG<<G4endl;
512 G4int PDG=2212; // *TMP* instead of pPDG
513 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG
514 if(!Z && N==1) // Change for Quasi-Elastic on neutron
515 {
516 Z=1;
517 N=0;
518 if (PDG==2212) PDG=2112;
519 else if(PDG==2112) PDG=2212;
520 }
521 G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP*
522 if(PDG==2212) xSec=PCSmanager->GetCrossSection(false, P, Z, N, PDG); // P CrossSect *TMP*
523 else xSec=NCSmanager->GetCrossSection(false, P, Z, N, PDG); // N CrossSect *TMP*
524#ifdef ppdebug
525 G4cout<<"G4QuasiFreeRat::Scatter: pPDG="<<pPDG<<",P="<<P<<",CS="<<xSec/millibarn<<G4endl;
526#endif
527#ifdef nandebug
528 if(xSec>0. || xSec<0. || xSec==0);
529 else G4cout<<"*Warning*G4QuasiFreeRatios::Scatter: xSec="<<xSec/millibarn<<G4endl;
530#endif
531 // @@ check a possibility to separate p, n, or alpha (!)
532 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
533 {
534#ifdef ppdebug
535 G4cerr<<"-Warning-G4QFR::Scat:**Zero XS**PDG="<<pPDG<<",NPDG="<<NPDG<<",P="<<P<<G4endl;
536#endif
537 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action
538 }
539 G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP*
540 if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP*
541 else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP*
542 if (mT==mDeut) mint/=two3d2;
543 else if(mT==mTrit || mT==mHel3) mint/=two3d3;
544 else if(mT==mAlph) mint/=two3d4;
545 G4double maxt=0.; // Prototype of max possible -t
546 if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t
547 else maxt=NCSmanager->GetHMaxT(); // max possible -t
548#ifdef ppdebug
549 G4cout<<"G4QFR::Scat:PDG="<<PDG<<",P="<<P<<",X="<<xSec<<",-t="<<mint<<"<"<<maxt<<", Z="
550 <<Z<<",N="<<N<<G4endl;
551#endif
552#ifdef nandebug
553 if(mint>-.0000001);
554 else G4cout<<"*Warning*G4QFR::Scat: -t="<<mint<<G4endl;
555#endif
556 G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS
557#ifdef ppdebug
558 G4cout<<"G4QFR::Scat:-t="<<mint<<"<"<<maxt<<", cost="<<cost<<", Z="<<Z<<",N="<<N<<G4endl;
559#endif
560 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.))
561 {
562 if (cost>1.) cost=1.;
563 else if(cost<-1.) cost=-1.;
564 else
565 {
566 G4double tm=0.;
567 if(PDG==2212) tm=PCSmanager->GetHMaxT();
568 else tm=NCSmanager->GetHMaxT();
569 G4cerr<<"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<tm<<G4endl;
570 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
571 }
572 }
573 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon
574 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01);
575 if(!G4QHadron(tot4M).RelDecayIn2(pr4M, reco4M, dir4M, cost, cost))
576 {
577 G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl;
578 //G4Exception("G4QFR::Scat:","009",FatalException,"Decay of ElasticComp");
579 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action
580 }
581#ifdef ppdebug
582 G4cout<<"G4QFR::Scat:p4M="<<pr4M<<"+r4M="<<reco4M<<",dr="<<dir4M<<",t4M="<<tot4M<<G4endl;
583#endif
584 return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result
585} // End of Scatter
G4QPDGCode::GetNuclMass
G4double GetNuclMass(G4int Z, G4int N, G4int S)
Definition: G4QPDGCode.cc:766

Referenced by G4QFragmentation::G4QFragmentation(), and G4QInelastic::PostStepDoIt().

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