49G4int G4QElastic::nPartCWorld=85;
50std::vector<G4int> G4QElastic::ElementZ;
51std::vector<G4double> G4QElastic::ElProbInMat;
52std::vector<std::vector<G4int>*> G4QElastic::ElIsoN;
53std::vector<std::vector<G4double>*>G4QElastic::IsoProbInEl;
62 G4cout<<
"G4QElastic::Constructor is called processName="<<processName<<
G4endl;
85 G4cout<<
"*W*G4QElastic::GetMeanFreePath: is called for notImplementedParticle"<<
G4endl;
90 G4cout<<
"G4QElastic::GetMeanFreePath: kinE="<<KinEn<<
",Mom="<<Momentum<<
G4endl;
122 <<
" isn't implemented (only SU(3) particles)"<<
G4endl;
124 G4cout<<
"G4QElastic::GetMeanFreePath:"<<nE<<
" Elem's in theMaterial,proj="<<pPDG<<
G4endl;
134 else if(pPDG== 130 || pPDG== 310)
142 else G4cout<<
"*Warning*G4QElastic::GetMeanFreePath: wrong PDG="<<pPDG<<
G4endl;
145 G4int IPIE=IsoProbInEl.size();
146 if(IPIE)
for(
G4int ip=0; ip<IPIE; ++ip)
148 std::vector<G4double>* SPI=IsoProbInEl[ip];
151 std::vector<G4int>* IsN=ElIsoN[ip];
159 for(
G4int i=0; i<nE; ++i)
161 G4Element* pElement=(*theElementVector)[i];
163 ElementZ.push_back(Z);
167 if(isoVector) isoSize=isoVector->size();
169 G4cout<<
"G4QElastic::GetMeanFreePath: isovectorLength="<<isoSize<<
G4endl;
179 std::vector<std::pair<G4int,G4double>*>* newAbund =
180 new std::vector<std::pair<G4int,G4double>*>;
182 for(
G4int j=0; j<isoSize; j++)
188 std::pair<G4int,G4double>* pr=
new std::pair<G4int,G4double>(N,abund);
190 G4cout<<
"G4QElastic::GetMeanFreePath:pair#="<<j<<
",N="<<N<<
",ab="<<abund<<
G4endl;
192 newAbund->push_back(pr);
195 G4cout<<
"G4QElastic::GetMeanFreePath: pairVectorLength="<<newAbund->size()<<
G4endl;
198 for(
G4int k=0; k<isoSize; k++)
delete (*newAbund)[k];
202 std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->
GetCSVector(Z,indEl);
203 std::vector<G4double>* SPI =
new std::vector<G4double>;
204 IsoProbInEl.push_back(SPI);
205 std::vector<G4int>* IsN =
new std::vector<G4int>;
206 ElIsoN.push_back(IsN);
207 G4int nIs=cs->size();
209 G4cout<<
"G4QEl::GMFP:=***=>,#isot="<<nIs<<
", Z="<<Z<<
", indEl="<<indEl<<
G4endl;
212 if(nIs)
for(
G4int j=0; j<nIs; j++)
214 std::pair<G4int,G4double>* curIs=(*cs)[j];
215 G4int N=curIs->first;
218 G4cout<<
"G4QEl::GMFP:*true*,P="<<Momentum<<
",Z="<<Z<<
",N="<<N<<
",PDG="<<pPDG<<
G4endl;
221 if(Q==-27.) ccsf=
false;
226 if(!CSmanager2) CSI=CSmanager->
GetCrossSection(ccsf,Momentum,Z,N,pPDG);
230 G4cout<<
"G4QEl::GMFP: jI="<<j<<
", Zt="<<Z<<
", Nt="<<N<<
", Mom="<<Momentum<<
", XSec="
235 SPI->push_back(susi);
242 ElProbInMat.push_back(sigma);
246 G4cout<<
"G4QElastic::GetMeanFreePath: MeanFreePath="<<1./sigma<<
G4endl;
248 if(sigma > 0.)
return 1./sigma;
294 static G4bool CWinit =
true;
304 G4cout<<
"G4QElastic::PostStepDoIt: Before the GetMeanFreePath is called In4M="
311 G4cout<<
"G4QElastic::PostStepDoIt: After the GetMeanFreePath is called"<<
G4endl;
317 if(std::fabs(Momentum-momentum)>.000001)
318 G4cerr<<
"*War*G4QElastic::PostStepDoIt:P(IU)="<<Momentum<<
"#"<<momentum<<
G4endl;
322 G4cout<<
"G4QElastic::PostStepDoIt: pP(IU)="<<Momentum<<
"="<<momentum<<
",pM="<<pM
323 <<
",scat4M="<<scat4M<<scat4M.
m()<<
G4endl;
327 G4cerr<<
"G4QElastic::PostStepDoIt: Only NA elastic is implemented."<<
G4endl;
334 G4cout<<
"G4QElastic::PostStepDoIt: "<<nE<<
" elements in the material."<<
G4endl;
377 G4cout<<
"G4QElastic::PostStepDoIt: projPDG="<<projPDG<<
", stPDG="<<prPDG<<
G4endl;
381 G4cerr<<
"*Warning*G4QElastic::PostStepDoIt: Undefined interacting hadron"<<
G4endl;
384 G4int EPIM=ElProbInMat.size();
386 G4cout<<
"G4QElastic::PostStDoIt:m="<<EPIM<<
",n="<<nE<<
",T="<<ElProbInMat[EPIM-1]<<
G4endl;
395 G4cout<<
"G4QElastic::PostStepDoIt:EPM["<<i<<
"]="<<ElProbInMat[i]<<
",r="<<rnd<<
G4endl;
397 if (rnd<ElProbInMat[i])
break;
401 G4Element* pElement=(*theElementVector)[i];
404 G4cout<<
"G4QElastic::PostStepDoIt: i="<<i<<
", Z(element)="<<Z<<
G4endl;
408 G4cerr<<
"---Warning---G4QElastic::PostStepDoIt: Element with Z="<<Z<<
G4endl;
411 std::vector<G4double>* SPI = IsoProbInEl[i];
412 std::vector<G4int>* IsN = ElIsoN[i];
413 G4int nofIsot=SPI->size();
415 G4cout<<
"G4QElastic::PosStDoIt: nI="<<nofIsot<<
",T="<<(*SPI)[nofIsot-1]<<
G4endl;
421 for(j=0; j<nofIsot; ++j)
424 G4cout<<
"G4QElastic::PostStepDoIt: SP["<<j<<
"]="<<(*SPI)[j]<<
", r="<<rndI<<
G4endl;
426 if(rndI < (*SPI)[j])
break;
428 if(j>=nofIsot) j=nofIsot-1;
430 G4int N =(*IsN)[j]; ;
432 G4cout<<
"G4QElastic::PostStepDoIt: j="<<i<<
", N(isotope)="<<N<<
", MeV="<<MeV<<
G4endl;
436 G4cerr<<
"*Warning*G4QElastic::PostStepDoIt: Isotope with Z="<<Z<<
", 0>N="<<N<<
G4endl;
441 G4cout<<
"G4QElastic::PostStepDoIt: N="<<N<<
" for element with Z="<<Z<<
G4endl;
445 G4cout<<
"G4QElastic::PostStepDoIt: track is initialized"<<
G4endl;
451 G4cout<<
"G4QElastic::PostStepDoIt: before Touchable extraction"<<
G4endl;
455 G4cout<<
"G4QElastic::PostStepDoIt: Touchable is extracted"<<
G4endl;
458 G4int targPDG=90000000+Z*1000+N;
465 G4cout<<
"G4QElastic::PostStepDoIt: tM="<<tM<<
", p4M="<<proj4M<<
", t4M="<<tot4M<<
G4endl;
467 EnMomConservation=tot4M;
476 else if(projPDG== 130 || projPDG== 310)
483 else if(projPDG>-3335 && projPDG<-1110)
485 else G4cout<<
"*Warning*G4QElastic::PostStepDoIt: wrong PDG="<<projPDG<<
G4endl;
487 G4cout<<
"G4QElas::PSDI:false,P="<<Momentum<<
",Z="<<Z<<
",N="<<N<<
",PDG="<<projPDG<<
G4endl;
490 if(!CSmanager2) xSec=CSmanager->
GetCrossSection(
false,Momentum,Z,N,projPDG);
494 G4cout<<
"G4QElast::PSDI:pPDG="<<projPDG<<
",P="<<Momentum<<
",CS="<<xSec/millibarn<<
G4endl;
497 if(xSec>0. || xSec<0. || xSec==0);
498 else G4cout<<
"*Warning*G4QElastic::PostStDoIt: xSec="<<xSec/millibarn<<
G4endl;
504 G4cerr<<
"*Warning*G4QElastic::PostStDoIt: *Zero cross-section* PDG="<<projPDG<<
",tPDG="
505 <<targPDG<<
",P="<<Momentum<<
G4endl;
515 G4cout<<
"G4QElast::PSDI:pPDG="<<projPDG<<
",tPDG="<<targPDG<<
",P="<<Momentum<<
",CS="
516 <<xSec<<
",-t="<<mint<<
G4endl;
520 else G4cout<<
"*Warning*G4QElastic::PostStDoIt:-t="<<mint<<
G4endl;
526 G4cout<<
"G4QElastic::PoStDoI:t="<<mint<<
", maxt="<<maxt<<
",Ek="<<kinEnergy<<
",tM="<<tM
527 <<
",pM="<<pM<<
",cost="<<cost<<
G4endl;
529 if(cost>1. || cost<-1. || !(cost>-1. || cost<1.))
531 if(cost>1.000001 || cost<-1.000001 || !(cost>-1. || cost<1.))
536 G4double twop2cm=(tM2+tM2)*(pEn*pEn-pM2)/sM;
537 G4cout<<
"*Warning*G4QElastic::PostStepDoIt:cos="<<cost<<
",t="<<mint<<
",T="<<kinEnergy
538 <<
",tM="<<tM<<
",tmax="<<2*kinEnergy*tM<<
",p="<<projPDG<<
",t="<<targPDG<<
G4endl;
539 G4cout<<
"..G4QElastic::PoStDoI: dpcm2="<<twop2cm<<
"="<<maxt<<
G4endl;
541 if (cost>1.) cost=1.;
542 else if(cost<-1.) cost=-1.;
546 if(!
G4QHadron(tot4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
548 G4cerr<<
"G4QElastic::PSD:t4M="<<tot4M<<
",pM="<<pM<<
",tM="<<tM<<
",cost="<<cost<<
G4endl;
552 G4cout<<
"G4QElastic::PoStDoIt:s4M="<<scat4M<<
"+r4M="<<reco4M<<
"="<<scat4M+reco4M<<
G4endl;
553 G4cout<<
"G4QElastic::PoStDoIt: scatE="<<scat4M.
e()-pM<<
", recoE="<<reco4M.
e()-tM<<
",d4M="
554 <<tot4M-scat4M-reco4M<<
G4endl;
561 if(finE<-1.e-8 || !(finE>-1.||finE<1.))
562 G4cerr<<
"*Warning*G4QElastic::PostStDoIt: Zero or negative scattered E="<<finE
563 <<
", s4M="<<scat4M<<
", r4M="<<reco4M<<
", d4M="<<tot4M-scat4M-reco4M<<
G4endl;
570 EnMomConservation-=scat4M;
575 G4cout<<
"G4QElastic::PostStepDoIt: Ion Z="<<Z<<
", A="<<aA<<
G4endl;
579 if(!theDefinition)
G4cout<<
"*Warning*G4QElastic::PostStepDoIt:drop PDG="<<targPDG<<
G4endl;
585 EnMomConservation-=reco4M;
587 G4cout<<
"G4QElastic::PostSDoIt:"<<targPDG<<reco4M<<reco4M.
m()<<EnMomConservation<<
G4endl;
594 G4cout<<
"G4QElastic::PSDoIt:p="<<curD<<curD.
mag()<<
",e="<<curE<<
",m="<<curM<<
G4endl;
602 G4cout<<
"G4QElastic::PostStepDoIt: **** PostStepDoIt is done ****"<<
G4endl;
639 else if(pPDG==2112) res=NCSmanager->
GetCrossSection(
true, p, Z, N, pPDG);
640 else if(pPDG== 211) res=PiPCSmanager->
GetCrossSection(
true, p, Z, N, pPDG);
641 else if(pPDG==-211) res=PiMCSmanager->
GetCrossSection(
true, p, Z, N, pPDG);
642 else if(pPDG== 321) res=KaPCSmanager->
GetCrossSection(
true, p, Z, N, pPDG);
643 else if(pPDG==-321) res=KaMCSmanager->
GetCrossSection(
true, p, Z, N, pPDG);
644 else if(pPDG==310||pPDG==130) res=(KaMCSmanager->
GetCrossSection(
true, p, Z, N, pPDG)+
646 else if(pPDG==3222) res=HyPCSmanager->
GetCrossSection(
true, p, Z, N, pPDG);
647 else if(pPDG>3110 && pPDG<3335) res=HypCSmanager->
GetCrossSection(
true, p, Z, N, pPDG);
648 else if(pPDG>-3335 && pPDG<-1110) res=aBaCSmanager->
GetCrossSection(
true,p,Z,N,pPDG);
649 else G4cout<<
"*Warning*G4QElastic::CalculateXSt: (o) wrong PDG="<<pPDG<<
G4endl;
654 else if(pPDG==2112) res=NCSmanager->
GetCrossSection(
false, p, Z, N, pPDG);
655 else if(pPDG== 211) res=PiPCSmanager->
GetCrossSection(
false, p, Z, N, pPDG);
656 else if(pPDG==-211) res=PiMCSmanager->
GetCrossSection(
false, p, Z, N, pPDG);
657 else if(pPDG== 321) res=KaPCSmanager->
GetCrossSection(
false, p, Z, N, pPDG);
658 else if(pPDG==-321) res=KaMCSmanager->
GetCrossSection(
false, p, Z, N, pPDG);
659 else if(pPDG==310||pPDG==130) res=(KaMCSmanager->
GetCrossSection(
false, p, Z, N, pPDG)+
661 else if(pPDG==3222) res=HyPCSmanager->
GetCrossSection(
false, p, Z, N, pPDG);
662 else if(pPDG>3110 && pPDG<3335) res=HypCSmanager->
GetCrossSection(
false, p, Z,N, pPDG);
663 else if(pPDG>-3335 && pPDG<-1110) res=aBaCSmanager->
GetCrossSection(
false,p,Z,N,pPDG);
664 else G4cout<<
"*Warning*G4QElastic::CalculateXSt: (x) wrong PDG="<<pPDG<<
G4endl;
682 if(oxs) res=PiPCSmanager->
GetHMaxT();
687 if(oxs) res=PiMCSmanager->
GetHMaxT();
690 else if(pPDG==321 || pPDG==310 || pPDG==130)
692 if(oxs) res=KaPCSmanager->
GetHMaxT();
697 if(oxs) res=KaMCSmanager->
GetHMaxT();
702 if(oxs) res=HyPCSmanager->
GetHMaxT();
705 else if(pPDG<3335 && pPDG>1110)
707 if(oxs) res=HypCSmanager->
GetHMaxT();
710 else if(pPDG>-3335 && pPDG<-1110)
712 if(oxs) res=aBaCSmanager->
GetHMaxT();
715 else G4cout<<
"*Warning*G4QElastic::CalculateXSt: (i) wrong PDG="<<pPDG<<
G4endl;
717 else G4cout<<
"*Warning*G4QElastic::CalculateXSt:*NotInitiatedScattering*"<<
G4endl;
std::vector< G4Element * > G4ElementVector
#define G4HadronicDeprecate(name)
std::vector< G4Isotope * > G4IsotopeVector
CLHEP::HepLorentzVector G4LorentzVector
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
static G4AntiLambda * AntiLambda()
static G4AntiNeutrinoMu * AntiNeutrinoMu()
static G4AntiNeutron * AntiNeutron()
static G4AntiOmegaMinus * AntiOmegaMinus()
static G4AntiProton * AntiProton()
static G4AntiSigmaMinus * AntiSigmaMinus()
static G4AntiSigmaPlus * AntiSigmaPlus()
static G4AntiSigmaZero * AntiSigmaZero()
static G4AntiXiMinus * AntiXiMinus()
static G4AntiXiZero * AntiXiZero()
const G4ThreeVector & GetMomentumDirection() const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4ParticleDefinition * GetDefinition() const
G4LorentzVector Get4Momentum() const
G4double GetKineticEnergy() const
void Set4Momentum(const G4LorentzVector &momentum)
G4double GetTotalMomentum() const
static G4Electron * Electron()
G4double * GetRelativeAbundanceVector() const
const G4Isotope * GetIsotope(G4int iso) const
G4IsotopeVector * GetIsotopeVector() const
static G4KaonMinus * KaonMinus()
static G4KaonPlus * KaonPlus()
static G4KaonZeroLong * KaonZeroLong()
static G4KaonZeroShort * KaonZeroShort()
static G4Lambda * Lambda()
const G4ElementVector * GetElementVector() const
size_t GetNumberOfElements() const
const G4double * GetVecNbOfAtomsPerVolume() const
static G4MuonMinus * MuonMinus()
static G4MuonPlus * MuonPlus()
static G4NeutrinoMu * NeutrinoMu()
static G4Neutron * Neutron()
static G4OmegaMinus * OmegaMinus()
void AddSecondary(G4Track *aSecondary)
void ProposeEnergy(G4double finalEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
virtual void Initialize(const G4Track &)
const G4String & GetParticleType() const
G4int GetPDGEncoding() const
const G4String & GetParticleName() const
const G4String & GetParticleSubType() const
G4ParticleDefinition * FindIon(G4int atomicNumber, G4int atomicMass, G4double excitationEnergy)
static G4ParticleTable * GetParticleTable()
static G4PionMinus * PionMinus()
static G4PionPlus * PionPlus()
static G4Positron * Positron()
static G4Proton * Proton()
static G4VQCrossSection * GetPointer()
static G4QCHIPSWorld * Get()
G4QParticleVector * GetParticles(G4int nOfParts=0)
G4LorentzVector GetEnegryMomentumConservation()
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
G4QElastic(const G4String &processName="CHIPSElasticScattering")
G4int GetNumberOfNeutronsInTarget()
G4bool IsApplicable(const G4ParticleDefinition &particle)
static G4VQCrossSection * GetPointer()
static G4VQCrossSection * GetPointer()
std::vector< std::pair< G4int, G4double > * > * GetCSVector(G4int Z, G4int index=0)
G4bool IsDefined(G4int Z, G4int Ind)
G4double GetMeanCrossSection(G4int Z, G4int index=0)
static G4QIsotope * Get()
G4int InitElement(G4int Z, G4int index, std::vector< std::pair< G4int, G4double > * > *abund)
static G4VQCrossSection * GetPointer()
static G4VQCrossSection * GetPointer()
static G4VQCrossSection * GetPointer()
static G4VQCrossSection * GetPointer()
static G4VQCrossSection * GetPointer()
static G4VQCrossSection * GetPointer()
static G4SigmaMinus * SigmaMinus()
static G4SigmaPlus * SigmaPlus()
static G4SigmaZero * SigmaZero()
static G4TauMinus * TauMinus()
static G4TauPlus * TauPlus()
G4double GetWeight() const
void SetWeight(G4double aValue)
const G4ThreeVector & GetPosition() const
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4double GetGlobalTime() const
G4Material * GetMaterial() const
const G4DynamicParticle * GetDynamicParticle() const
const G4TouchableHandle & GetTouchableHandle() const
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4ParticleChange aParticleChange
const G4String & GetProcessName() const
virtual G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
virtual G4double GetCrossSection(G4bool, G4double, G4int, G4int, G4int pPDG=0)
virtual G4double GetHMaxT()
static G4XiMinus * XiMinus()
static G4XiZero * XiZero()