#include <G4QuasiElRatios.hh>
Definition at line 53 of file G4QuasiElRatios.hh.
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 }
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
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
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