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

#include <G4QLowEnergy.hh>

+ Inheritance diagram for G4QLowEnergy:

Public Member Functions

 G4QLowEnergy (const G4String &processName="CHIPS_LowEnergyIonIonInelastic")
 
 ~G4QLowEnergy ()
 
G4bool IsApplicable (const G4ParticleDefinition &particle)
 
G4double GetMeanFreePath (const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
 
G4VParticleChange * PostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
 
G4LorentzVector GetEnegryMomentumConservation ()
 
G4int GetNumberOfNeutronsInTarget ()
 
void SwitchOnEvaporation ()
 
void SwitchOffEvaporation ()
 
- Public Member Functions inherited from G4VDiscreteProcess
 G4VDiscreteProcess (const G4String &, G4ProcessType aType=fNotDefined)
 
 G4VDiscreteProcess (G4VDiscreteProcess &)
 
virtual ~G4VDiscreteProcess ()
 
virtual G4double PostStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
 
virtual G4VParticleChange * PostStepDoIt (const G4Track &, const G4Step &)
 
virtual G4double AlongStepGetPhysicalInteractionLength (const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
 
virtual G4double AtRestGetPhysicalInteractionLength (const G4Track &, G4ForceCondition *)
 
virtual G4VParticleChange * AtRestDoIt (const G4Track &, const G4Step &)
 
virtual G4VParticleChange * AlongStepDoIt (const G4Track &, const G4Step &)
 
- Public Member Functions inherited from G4VProcess
 G4VProcess (const G4String &aName="NoName", G4ProcessType aType=fNotDefined)
 
 G4VProcess (const G4VProcess &right)
 
virtual ~G4VProcess ()
 
G4int operator== (const G4VProcess &right) const
 
G4int operator!= (const G4VProcess &right) const
 
virtual G4VParticleChange * PostStepDoIt (const G4Track &track, const G4Step &stepData)=0
 
virtual G4VParticleChange * AlongStepDoIt (const G4Track &track, const G4Step &stepData)=0
 
virtual G4VParticleChange * AtRestDoIt (const G4Track &track, const G4Step &stepData)=0
 
virtual G4double AlongStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &proposedSafety, G4GPILSelection *selection)=0
 
virtual G4double AtRestGetPhysicalInteractionLength (const G4Track &track, G4ForceCondition *condition)=0
 
virtual G4double PostStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)=0
 
G4double GetCurrentInteractionLength () const
 
void SetPILfactor (G4double value)
 
G4double GetPILfactor () const
 
G4double AlongStepGPIL (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &proposedSafety, G4GPILSelection *selection)
 
G4double AtRestGPIL (const G4Track &track, G4ForceCondition *condition)
 
G4double PostStepGPIL (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
 
virtual G4bool IsApplicable (const G4ParticleDefinition &)
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void PreparePhysicsTable (const G4ParticleDefinition &)
 
virtual G4bool StorePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)
 
virtual G4bool RetrievePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)
 
const G4String & GetPhysicsTableFileName (const G4ParticleDefinition *, const G4String &directory, const G4String &tableName, G4bool ascii=false)
 
const G4String & GetProcessName () const
 
G4ProcessType GetProcessType () const
 
void SetProcessType (G4ProcessType)
 
G4int GetProcessSubType () const
 
void SetProcessSubType (G4int)
 
virtual void StartTracking (G4Track *)
 
virtual void EndTracking ()
 
virtual void SetProcessManager (const G4ProcessManager *)
 
virtual const G4ProcessManager * GetProcessManager ()
 
virtual void ResetNumberOfInteractionLengthLeft ()
 
G4double GetNumberOfInteractionLengthLeft () const
 
G4double GetTotalNumberOfInteractionLengthTraversed () const
 
G4bool isAtRestDoItIsEnabled () const
 
G4bool isAlongStepDoItIsEnabled () const
 
G4bool isPostStepDoItIsEnabled () const
 
virtual void DumpInfo () const
 
void SetVerboseLevel (G4int value)
 
G4int GetVerboseLevel () const
 

Additional Inherited Members

- Static Public Member Functions inherited from G4VProcess
static const G4String & GetProcessTypeName (G4ProcessType)
 
virtual G4double GetMeanFreePath (const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)=0
 
- Protected Member Functions inherited from G4VProcess
void SubtractNumberOfInteractionLengthLeft (G4double previousStepSize)
 
void ClearNumberOfInteractionLengthLeft ()
 
- Protected Attributes inherited from G4VProcess
const G4ProcessManager * aProcessManager
 
G4VParticleChange * pParticleChange
 
G4ParticleChange aParticleChange
 
G4double theNumberOfInteractionLengthLeft
 
G4double currentInteractionLength
 
G4double theInitialNumberOfInteractionLength
 
G4String theProcessName
 
G4String thePhysicsTableFileName
 
G4ProcessType theProcessType
 
G4int theProcessSubType
 
G4double thePILfactor
 
G4bool enableAtRestDoIt
 
G4bool enableAlongStepDoIt
 
G4bool enablePostStepDoIt
 
G4int verboseLevel
 

Detailed Description

Definition at line 72 of file G4QLowEnergy.hh.

Constructor & Destructor Documentation

◆ G4QLowEnergy()

G4QLowEnergy::G4QLowEnergy ( const G4String &  processName = "CHIPS_LowEnergyIonIonInelastic")

Definition at line 60 of file G4QLowEnergy.cc.

60 :
61 G4VDiscreteProcess(processName, fHadronic), evaporate(true)
62{
63 G4HadronicDeprecate("G4QLowEnergy");
64
65#ifdef debug
66 G4cout<<"G4QLowEnergy::Constructor is called processName="<<processName<<G4endl;
67#endif
68 if (verboseLevel>0) G4cout<<GetProcessName()<<" process is created "<<G4endl;
69 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part. max)
70}
G4HadronicDeprecate
#define G4HadronicDeprecate(name)
Definition: G4HadronicDeprecate.hh:33
fHadronic
@ fHadronic
Definition: G4ProcessType.hh:49
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4QCHIPSWorld::Get
static G4QCHIPSWorld * Get()
Definition: G4QCHIPSWorld.cc:79
G4QCHIPSWorld::GetParticles
G4QParticleVector * GetParticles(G4int nOfParts=0)
Definition: G4QCHIPSWorld.cc:95
G4VProcess::verboseLevel
G4int verboseLevel
Definition: G4VProcess.hh:368
G4VProcess::GetProcessName
const G4String & GetProcessName() const
Definition: G4VProcess.hh:379

◆ ~G4QLowEnergy()

G4QLowEnergy::~G4QLowEnergy ( )

Definition at line 73 of file G4QLowEnergy.cc.

73{}

Member Function Documentation

◆ GetEnegryMomentumConservation()

G4LorentzVector G4QLowEnergy::GetEnegryMomentumConservation ( )

Definition at line 76 of file G4QLowEnergy.cc.

76{return EnMomConservation;}

◆ GetMeanFreePath()

G4double G4QLowEnergy::GetMeanFreePath ( const G4Track &  aTrack,
G4double  previousStepSize,
G4ForceCondition *  condition 
)
virtual

Implements G4VDiscreteProcess.

Definition at line 83 of file G4QLowEnergy.cc.

84{
85 *F = NotForced;
86 const G4DynamicParticle* incidentParticle = Track.GetDynamicParticle();
87 G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
88 if( !IsApplicable(*incidentParticleDefinition))
89 G4cout<<"-Warning-G4QLowEnergy::GetMeanFreePath for notImplemented Particle"<<G4endl;
90 // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
91 G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
92#ifdef debug
93 G4double KinEn = incidentParticle->GetKineticEnergy();
94 G4cout<<"G4QLowEnergy::GetMeanFreePath:Prpj, kinE="<<KinEn<<", Mom="<<Momentum<<G4endl;
95#endif
96 const G4Material* material = Track.GetMaterial(); // Get the current material
97 const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
98 const G4ElementVector* theElementVector = material->GetElementVector();
99 G4int nE=material->GetNumberOfElements();
100#ifdef debug
101 G4cout<<"G4QLowEnergy::GetMeanFreePath:"<<nE<<" Elems"<<G4endl;
102#endif
103 G4int pPDG=0;
104 G4int Z=static_cast<G4int>(incidentParticleDefinition->GetPDGCharge());
105 G4int A=incidentParticleDefinition->GetBaryonNumber();
106 if ( incidentParticleDefinition == G4Proton::Proton() ) pPDG = 2212;
107 else if ( incidentParticleDefinition == G4Deuteron::Deuteron() ) pPDG = 1000010020;
108 else if ( incidentParticleDefinition == G4Alpha::Alpha() ) pPDG = 1000020040;
109 else if ( incidentParticleDefinition == G4Triton::Triton() ) pPDG = 1000010030;
110 else if ( incidentParticleDefinition == G4He3::He3() ) pPDG = 1000020030;
111 else if ( incidentParticleDefinition == G4GenericIon::GenericIon() || (Z > 0 && A > 0))
112 {
113 pPDG=incidentParticleDefinition->GetPDGEncoding();
114#ifdef debug
115 G4int prPDG=1000000000+10000*A+10*Z;
116 G4cout<<"G4QIonIonElastic::GetMeanFreePath: PDG="<<prPDG<<"="<<pPDG<<G4endl;
117#endif
118 }
119 else G4cout<<"-Warning-G4QLowEnergy::GetMeanFreePath: only AA & pA implemented"<<G4endl;
120 G4VQCrossSection* CSmanager=G4QIonIonCrossSection::GetPointer();
121 if(pPDG == 2212) CSmanager=G4QProtonNuclearCrossSection::GetPointer(); // just to test
122 Momentum/=incidentParticleDefinition->GetBaryonNumber(); // Divide Mom by projectile A
123 G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
124 G4double sigma=0.; // Sums over elements for the material
125 G4int IPIE=IsoProbInEl.size(); // How many old elements?
126 if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
127 {
128 std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
129 SPI->clear();
130 delete SPI;
131 std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector
132 IsN->clear();
133 delete IsN;
134 }
135 ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE)
136 ElementZ.clear(); // Clear the body vector for Z of Elements
137 IsoProbInEl.clear(); // Clear the body vector for SPI
138 ElIsoN.clear(); // Clear the body vector for N of Isotopes
139 for(G4int i=0; i<nE; ++i)
140 {
141 G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
142 Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
143 ElementZ.push_back(Z); // Remember Z of the Element
144 G4int isoSize=0; // The default for the isoVectorLength is 0
145 G4int indEl=0; // Index of non-natural element or 0(default)
146 G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
147 if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
148#ifdef debug
149 G4cout<<"G4QLowEnergy::GetMeanFreePath: isovector Length="<<isoSize<<G4endl;
150#endif
151 if(isoSize) // The Element has non-trivial abundance set
152 {
153 indEl=pElement->GetIndex()+1; // Index of the non-trivial element is an order
154#ifdef debug
155 G4cout<<"G4QLowEn::GetMFP:iE="<<indEl<<",def="<<Isotopes->IsDefined(Z,indEl)<<G4endl;
156#endif
157 if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define
158 {
159 std::vector<std::pair<G4int,G4double>*>* newAbund =
160 new std::vector<std::pair<G4int,G4double>*>;
161 G4double* abuVector=pElement->GetRelativeAbundanceVector();
162 for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes
163 {
164 G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
165 if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QDiffract::GetMeanFreePath: Z="
166 <<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
167 G4double abund=abuVector[j];
168 std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
169#ifdef debug
170 G4cout<<"G4QLowEnergy::GetMeanFreePath:pair#"<<j<<",N="<<N<<",a="<<abund<<G4endl;
171#endif
172 newAbund->push_back(pr);
173 }
174#ifdef debug
175 G4cout<<"G4QLowEnergy::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl;
176#endif
177 indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
178 for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary
179 delete newAbund; // Was "new" in the beginning of the name space
180 }
181 }
182 std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
183 std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
184 IsoProbInEl.push_back(SPI);
185 std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
186 ElIsoN.push_back(IsN);
187 G4int nIs=cs->size(); // A#Of Isotopes in the Element
188#ifdef debug
189 G4cout<<"G4QLowEnergy::GetMFP:***=>,#isot="<<nIs<<", Z="<<Z<<", indEl="<<indEl<<G4endl;
190#endif
191 G4double susi=0.; // sum of CS over isotopes
192 if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El
193 {
194 std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
195 G4int N=curIs->first; // #of Neuterons in the isotope j of El i
196 IsN->push_back(N); // Remember Min N for the Element
197#ifdef debug
198 G4cout<<"G4QLoE::GMFP:true,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",pPDG="<<pPDG<<G4endl;
199#endif
200 G4bool ccsf=true; // Extract inelastic Ion-Ion cross-section
201#ifdef debug
202 G4cout<<"G4QLowEnergy::GMFP: GetCS #1 j="<<j<<G4endl;
203#endif
204 G4double CSI=CSmanager->GetCrossSection(ccsf,Momentum,Z,N,pPDG);//CS(j,i) for isotope
205#ifdef debug
206 G4cout<<"G4QLowEnergy::GetMeanFreePath: jI="<<j<<", Zt="<<Z<<", Nt="<<N<<", Mom="
207 <<Momentum<<", XSec="<<CSI/millibarn<<G4endl;
208#endif
209 curIs->second = CSI;
210 susi+=CSI; // Make a sum per isotopes
211 SPI->push_back(susi); // Remember summed cross-section
212 } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
213 sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV)
214#ifdef debug
215 G4cout<<"G4QLowEnergy::GetMeanFreePath:<XS>="<<Isotopes->GetMeanCrossSection(Z,indEl)
216 <<",AddSigm="<<Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i]<<G4endl;
217#endif
218 ElProbInMat.push_back(sigma);
219 } // End of LOOP over Elements
220 // Check that cross section is not zero and return the mean free path
221#ifdef debug
222 G4cout<<"G4QLowEnergy::GetMeanFreePath: MeanFreePath="<<1./sigma<<G4endl;
223#endif
224 if(sigma > 0.000000001) return 1./sigma; // Mean path [distance]
225 return DBL_MAX;
226}
G4ElementVector
std::vector< G4Element * > G4ElementVector
Definition: G4ElementVector.hh:44
NotForced
@ NotForced
Definition: G4ForceCondition.hh:56
G4IsotopeVector
std::vector< G4Isotope * > G4IsotopeVector
Definition: G4IsotopeVector.hh:44
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4cerr
G4DLLIMPORT std::ostream G4cerr
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
G4Element::GetRelativeAbundanceVector
G4double * GetRelativeAbundanceVector() const
Definition: G4Element.hh:166
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4Element::GetIsotope
const G4Isotope * GetIsotope(G4int iso) const
Definition: G4Element.hh:169
G4Element::GetIndex
size_t GetIndex() const
Definition: G4Element.hh:182
G4Element::GetIsotopeVector
G4IsotopeVector * GetIsotopeVector() const
Definition: G4Element.hh:162
G4GenericIon::GenericIon
static G4GenericIon * GenericIon()
Definition: G4GenericIon.cc:92
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:94
G4Isotope::GetZ
G4int GetZ() const
Definition: G4Isotope.hh:91
G4Isotope::GetN
G4int GetN() const
Definition: G4Isotope.hh:94
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:189
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:185
G4Material::GetVecNbOfAtomsPerVolume
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:205
G4Material::GetMaterial
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:576
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:119
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4QIonIonCrossSection::GetPointer
static G4VQCrossSection * GetPointer()
Definition: G4QIonIonCrossSection.cc:65
G4QIsotope::GetCSVector
std::vector< std::pair< G4int, G4double > * > * GetCSVector(G4int Z, G4int index=0)
Definition: G4QIsotope.cc:1737
G4QIsotope::IsDefined
G4bool IsDefined(G4int Z, G4int Ind)
Definition: G4QIsotope.cc:1660
G4QIsotope::GetMeanCrossSection
G4double GetMeanCrossSection(G4int Z, G4int index=0)
Definition: G4QIsotope.cc:1842
G4QIsotope::Get
static G4QIsotope * Get()
Definition: G4QIsotope.cc:720
G4QIsotope::InitElement
G4int InitElement(G4int Z, G4int index, std::vector< std::pair< G4int, G4double > * > *abund)
Definition: G4QIsotope.cc:1557
G4QLowEnergy::IsApplicable
G4bool IsApplicable(const G4ParticleDefinition &particle)
Definition: G4QLowEnergy.cc:228
G4QProtonNuclearCrossSection::GetPointer
static G4VQCrossSection * GetPointer()
Definition: G4QProtonNuclearCrossSection.cc:63
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4VQCrossSection::GetCrossSection
virtual G4double GetCrossSection(G4bool, G4double, G4int, G4int, G4int pPDG=0)
Definition: G4VQCrossSection.hh:99
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83

◆ GetNumberOfNeutronsInTarget()

G4int G4QLowEnergy::GetNumberOfNeutronsInTarget ( )

Definition at line 78 of file G4QLowEnergy.cc.

78{return nOfNeutrons;}

◆ IsApplicable()

G4bool G4QLowEnergy::IsApplicable ( const G4ParticleDefinition &  particle)
virtual

Reimplemented from G4VProcess.

Definition at line 228 of file G4QLowEnergy.cc.

229{
230 G4int Z=static_cast<G4int>(particle.GetPDGCharge());
231 G4int A=particle.GetBaryonNumber();
232 if (particle == *( G4Proton::Proton() )) return true;
233 else if (particle == *( G4Neutron::Neutron() )) return true;
234 else if (particle == *( G4Deuteron::Deuteron() )) return true;
235 else if (particle == *( G4Alpha::Alpha() )) return true;
236 else if (particle == *( G4Triton::Triton() )) return true;
237 else if (particle == *( G4He3::He3() )) return true;
238 else if (particle == *( G4GenericIon::GenericIon() )) return true;
239 else if (Z > 0 && A > 0) return true;
240#ifdef debug
241 G4cout<<"***>>G4QLowEnergy::IsApplicable: projPDG="<<particle.GetPDGEncoding()<<", A="
242 <<A<<", Z="<<Z<<G4endl;
243#endif
244 return false;
245}
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104

Referenced by GetMeanFreePath(), and PostStepDoIt().

◆ PostStepDoIt()

G4VParticleChange * G4QLowEnergy::PostStepDoIt ( const G4Track &  aTrack,
const G4Step &  aStep 
)
virtual

Reimplemented from G4VDiscreteProcess.

Definition at line 247 of file G4QLowEnergy.cc.

248{
249 static const G4double mProt= G4QPDGCode(2212).GetMass()/MeV; // CHIPS proton Mass in MeV
250 static const G4double mPro2= mProt*mProt; // CHIPS sq proton Mass
251 static const G4double mNeut= G4QPDGCode(2112).GetMass()/MeV; // CHIPS neutron Mass in MeV
252 static const G4double mNeu2= mNeut*mNeut; // CHIPS sq neutron Mass
253 static const G4double mLamb= G4QPDGCode(3122).GetMass()/MeV; // CHIPS Lambda Mass in MeV
254 static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0)/MeV;
255 static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0)/MeV;
256 static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0)/MeV;
257 static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0)/MeV;
258 static const G4double mFm= 250*MeV;
259 static const G4double third= 1./3.;
260 static const G4ThreeVector zeroMom(0.,0.,0.);
261 static G4ParticleDefinition* aGamma = G4Gamma::Gamma();
262 static G4ParticleDefinition* aPiZero = G4PionZero::PionZero();
263 static G4ParticleDefinition* aPiPlus = G4PionPlus::PionPlus();
264 static G4ParticleDefinition* aPiMinus = G4PionMinus::PionMinus();
265 static G4ParticleDefinition* aProton = G4Proton::Proton();
266 static G4ParticleDefinition* aNeutron = G4Neutron::Neutron();
267 static G4ParticleDefinition* aLambda = G4Lambda::Lambda();
268 static G4ParticleDefinition* aDeuteron = G4Deuteron::Deuteron();
269 static G4ParticleDefinition* aTriton = G4Triton::Triton();
270 static G4ParticleDefinition* aHe3 = G4He3::He3();
271 static G4ParticleDefinition* anAlpha = G4Alpha::Alpha();
272 static const G4int nCh=26; // #of combinations
273 static G4QNucleus Nuc; // A fake nucleus to call Evaporation
274 //
275 //-------------------------------------------------------------------------------------
276 static G4bool CWinit = true; // CHIPS Warld needs to be initted
277 if(CWinit)
278 {
279 CWinit=false;
280 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
281 }
282 //-------------------------------------------------------------------------------------
283 const G4DynamicParticle* projHadron = track.GetDynamicParticle();
284 const G4ParticleDefinition* particle=projHadron->GetDefinition();
285#ifdef pdebug
286 G4cout<<"G4QLowEnergy::PostStepDoIt: *** Called *** In4M="
287 <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type="
288 <<particle->GetParticleType()<<",SubType="<<particle->GetParticleSubType()<<G4endl;
289#endif
290 //G4ForceCondition cond=NotForced;
291 //GetMeanFreePath(track, -27., &cond); // @@ ?? jus to update parameters?
292#ifdef debug
293 G4cout<<"G4QLowEnergy::PostStepDoIt: After GetMeanFreePath is called"<<G4endl;
294#endif
295 std::vector<G4Track*> result;
296 G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV; // Convert to MeV!
297 G4double momentum = projHadron->GetTotalMomentum()/MeV; // 3-momentum of the Proj in MeV
298 G4double Momentum = proj4M.rho(); // @@ Just for the test purposes
299 if(std::fabs(Momentum-momentum)>.000001)
300 G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt:P_IU="<<Momentum<<"#"<<momentum<<G4endl;
301#ifdef debug
302 G4cout<<"G4QLowEnergy::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum<<",proj4M,m="
303 <<proj4M<<proj4M.m()<<G4endl;
304#endif
305 if (!IsApplicable(*particle)) // Check applicability
306 {
307 G4cerr<<"G4QLowEnergy::PostStepDoIt: Only NA is implemented."<<G4endl;
308 return 0;
309 }
310 const G4Material* material = track.GetMaterial(); // Get the current material
311 const G4ElementVector* theElementVector = material->GetElementVector();
312 G4int nE=material->GetNumberOfElements();
313#ifdef debug
314 G4cout<<"G4QLowEnergy::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
315#endif
316 G4int projPDG=0; // PDG Code prototype for the captured hadron
317 // Not all these particles are implemented yet (see Is Applicable)
318 G4int Z=static_cast<G4int>(particle->GetPDGCharge());
319 G4int A=particle->GetBaryonNumber();
320 if (particle == G4Proton::Proton() ) projPDG= 2212;
321 else if (particle == G4Neutron::Neutron() ) projPDG= 2112;
322 else if (particle == G4Deuteron::Deuteron() ) projPDG= 1000010020;
323 else if (particle == G4Alpha::Alpha() ) projPDG= 1000020040;
324 else if (particle == G4Triton::Triton() ) projPDG= 1000010030;
325 else if (particle == G4He3::He3() ) projPDG= 1000020030;
326 else if (particle == G4GenericIon::GenericIon() || (Z > 0 && A > 0))
327 {
328 projPDG=particle->GetPDGEncoding();
329#ifdef debug
330 G4int prPDG=1000000000+10000*Z+10*A; // just for the testing purposes
331 G4cout<<"G4QLowEnergy::PostStepDoIt: PDG="<<prPDG<<"="<<projPDG<<G4endl;
332#endif
333 }
334 else G4cout<<"-Warning-G4QLowEnergy::PostStepDoIt:Unknown projectile Ion"<<G4endl;
335#ifdef pdebug
336 G4int prPDG=particle->GetPDGEncoding();
337 G4cout<<"G4QLowEnergy::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
338#endif
339 if(!projPDG)
340 {
341 G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt:UndefProjHadron(PDG=0) ->ret 0"<<G4endl;
342 return 0;
343 }
344 // Element treatment
345 G4int EPIM=ElProbInMat.size();
346#ifdef debug
347 G4cout<<"G4QLowEn::PostStDoIt: m="<<EPIM<<", n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
348#endif
349 G4int i=0;
350 if(EPIM>1)
351 {
352 G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
353 for(i=0; i<nE; ++i)
354 {
355#ifdef debug
356 G4cout<<"G4QLowEn::PostStepDoIt: EPM["<<i<<"]="<<ElProbInMat[i]<<", r="<<rnd<<G4endl;
357#endif
358 if (rnd<ElProbInMat[i]) break;
359 }
360 if(i>=nE) i=nE-1; // Top limit for the Element
361 }
362 G4Element* pElement=(*theElementVector)[i];
363 G4int tZ=static_cast<G4int>(pElement->GetZ());
364#ifdef debug
365 G4cout<<"G4QLowEnergy::PostStepDoIt: i="<<i<<", Z(element)="<<tZ<<G4endl;
366#endif
367 if(tZ<=0)
368 {
369 G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Element with Z="<<tZ<<G4endl;
370 if(tZ<0) return 0;
371 }
372 std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
373 std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i]
374 G4int nofIsot=SPI->size(); // #of isotopes in the element i
375#ifdef debug
376 G4cout<<"G4QLowEnergy::PostStepDoIt: nI="<<nofIsot<<", T="<<(*SPI)[nofIsot-1]<<G4endl;
377#endif
378 G4int j=0;
379 if(nofIsot>1)
380 {
381 G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element
382 for(j=0; j<nofIsot; ++j)
383 {
384#ifdef debug
385 G4cout<<"G4QLowEnergy::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<",r="<<rndI<<G4endl;
386#endif
387 if(rndI < (*SPI)[j]) break;
388 }
389 if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope
390 }
391 G4int tN =(*IsN)[j]; ; // Randomized number of neutrons
392#ifdef debug
393 G4cout<<"G4QLowEnergy::PostStepDoIt: j="<<i<<", N(isotope)="<<tN<<", MeV="<<MeV<<G4endl;
394#endif
395 if(tN<0)
396 {
397 G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Isotope Z="<<tZ<<" has 0>N="<<tN<<G4endl;
398 return 0;
399 }
400 nOfNeutrons=tN; // Remember it for the energy-momentum check
401#ifdef debug
402 G4cout<<"G4QLowEnergy::PostStepDoIt: N="<<tN<<" for element with Z="<<tZ<<G4endl;
403#endif
404 if(tN<0)
405 {
406 G4cerr<<"***Warning***G4QLowEnergy::PostStepDoIt: Element with N="<<tN<< G4endl;
407 return 0;
408 }
409 aParticleChange.Initialize(track);
410#ifdef debug
411 G4cout<<"G4QLowEnergy::PostStepDoIt: track is initialized"<<G4endl;
412#endif
413 G4double weight = track.GetWeight();
414 G4double localtime = track.GetGlobalTime();
415 G4ThreeVector position = track.GetPosition();
416#ifdef debug
417 G4cout<<"G4QLowEnergy::PostStepDoIt: before Touchable extraction"<<G4endl;
418#endif
419 G4TouchableHandle trTouchable = track.GetTouchableHandle();
420#ifdef debug
421 G4cout<<"G4QLowEnergy::PostStepDoIt: Touchable is extracted"<<G4endl;
422#endif
423 G4int targPDG=90000000+tZ*1000+tN; // Target PDG code
424 G4QPDGCode targQPDG(targPDG); // To get CHIPS properties
425 G4double tM=targQPDG.GetMass(); // CHIPS target nucleus mass in MeV
426 G4int pL=particle->GetQuarkContent(3)-particle->GetAntiQuarkContent(3); // Strangeness
427 G4int pZ=static_cast<G4int>(particle->GetPDGCharge()); // Charge of the projectile
428 G4int pN=particle->GetBaryonNumber()-pZ-pL;// #of neutrons in projectile
429 G4double pM=targQPDG.GetNuclMass(pZ,pN,0); // CHIPS projectile nucleus mass in MeV
430 G4double cosp=-14*Momentum*(tM-pM)/tM/pM; // Asymmetry power for angular distribution
431#ifdef debug
432 G4cout<<"G4QLowEnergy::PoStDoIt: Proj("<<pZ<<","<<pN<<","<<pL<<")p="<<pM<<",Targ=("<<tZ
433 <<","<<tN<<"), cosp="<<cosp<<G4endl;
434#endif
435 G4double kinEnergy= projHadron->GetKineticEnergy()*MeV; // Kin energy in MeV (Is *MeV n?)
436 G4ParticleMomentum dir = projHadron->GetMomentumDirection();// It is a unit three-vector
437 G4LorentzVector targ4M=G4LorentzVector(0.,0.,0.,tM); // Target's 4-mom
438 G4LorentzVector tot4M =proj4M+targ4M; // Total 4-mom of the reaction
439#ifdef pdebug
440 G4cout<<"G4QLowEnergy::PostStepDoIt: tM="<<tM<<",p4M="<<proj4M<<",t4M="<<tot4M<<G4endl;
441#endif
442 EnMomConservation=tot4M; // Total 4mom of reaction for E/M conservation
443 // --- Start of Coulomb barrier check ---
444 G4double dER = tot4M.e() - tM - pM; // Energy of the reaction
445 G4double pA = pZ+pN; // Atomic weight of the projectile (chged blw)
446 G4double tA = tZ+tN; // Atomic weight of the target (changed below)
447 G4double CBE = Nuc.CoulombBarGen(tZ, tA, pZ, pA); // Coulomb Barrier of the reaction
448#ifdef debug
449 G4cout<<"G4QLowEnergy::PoStDoIt: Proj("<<pZ<<","<<pN<<","<<pL<<")p="<<pM<<",Targ=("<<tZ
450 <<","<<tN<<"), dE="<<dER<<", CB="<<CBE<<G4endl;
451#endif
452 if(dER<CBE) // The cross-section iz 0 -> Do Nothing
453 {
454#ifdef debug
455 G4cerr<<"-Warning-G4QLowEnergy::PSDoIt: *Below Coulomb barrier* PDG="<<projPDG
456 <<",Z="<<tZ<<",tN="<<tN<<",P="<<Momentum<<G4endl;
457#endif
458 //Do Nothing Action insead of the reaction
459 aParticleChange.ProposeEnergy(kinEnergy);
460 aParticleChange.ProposeLocalEnergyDeposit(0.);
461 aParticleChange.ProposeMomentumDirection(dir) ;
462 return G4VDiscreteProcess::PostStepDoIt(track,step);
463 }
464 // --- End of Coulomb barrier check ---
465 // @@ Probably this is not necessary any more
466#ifdef debug
467 G4cout<<"G4QLE::PSDI:false,P="<<Momentum<<",Z="<<pZ<<",N="<<pN<<",PDG="<<projPDG<<G4endl;
468#endif
469 G4double xSec=CalculateXS(Momentum, tZ, tN, projPDG); // Recalculate CrossSection
470#ifdef pdebug
471 G4cout<<"G4QLowEn::PSDI:PDG="<<projPDG<<",P="<<Momentum<<",tZ="<<tZ<<",N="<<tN<<",XS="
472 <<xSec/millibarn<<G4endl;
473#endif
474#ifdef nandebug
475 if(xSec>0. || xSec<0. || xSec==0);
476 else G4cout<<"-Warning-G4QLowEnergy::PostSDI: *NAN* xSec="<<xSec/millibarn<<G4endl;
477#endif
478 // @@ check a possibility to separate p, n, or alpha (!)
479 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
480 {
481#ifdef debug
482 G4cerr<<"-Warning-G4QLowEnergy::PSDoIt: *Zero cross-section* PDG="<<projPDG
483 <<",Z="<<tZ<<",tN="<<tN<<",P="<<Momentum<<G4endl;
484#endif
485 //Do Nothing Action insead of the reaction
486 aParticleChange.ProposeEnergy(kinEnergy);
487 aParticleChange.ProposeLocalEnergyDeposit(0.);
488 aParticleChange.ProposeMomentumDirection(dir) ;
489 return G4VDiscreteProcess::PostStepDoIt(track,step);
490 }
491 // Kill interacting hadron
492 aParticleChange.ProposeEnergy(0.);
493 aParticleChange.ProposeTrackStatus(fStopAndKill);
494 G4int tB=tZ+tN;
495 G4int pB=pZ+pN;
496#ifdef pdebug
497 G4cout<<"G4QLowEn::PSDI: Projectile track is killed"<<", tA="<<tB<<", pA="<<pB<<G4endl;
498#endif
499 // algorithm implementation STARTS HERE --- All calculations are in IU --------
500 tA=tB;
501 pA=pB;
502 G4double tR=1.1; // target nucleus R in fm
503 if(tB > 1) tR*=std::pow(tA,third); // in fm
504 G4double pR=1.1; // projectile nucleus R in fm
505 if(pB > 1) pR*=std::pow(pA,third); // in fm
506 G4double R=tR+pR; // total radius
507 G4double R2=R*R;
508 G4int tD=0;
509 G4int pD=0;
510 G4int nAt=0;
511 G4int nAtM=27;
512 G4int nSec = 0;
513 G4double tcM=0.;
514 G4double tnM=1.;
515#ifdef edebug
516 G4int totChg=0;
517 G4int totBaN=0;
518 G4LorentzVector tch4M =tot4M; // Total 4-mom of the reaction
519#endif
520 while((!tD && !pD && ++nAt<nAtM) || tcM<tnM)
521 {
522#ifdef edebug
523 totChg=tZ+pZ;
524 totBaN=tB+pB;
525 tch4M =tot4M;
526 G4cout<<">G4QLEn::PSDI:#"<<nAt<<",tC="<<totChg<<",tA="<<totBaN<<",t4M="<<tch4M<<G4endl;
527#endif
528 G4LorentzVector tt4M=tot4M;
529 G4int resultSize=result.size();
530 if(resultSize)
531 {
532 for(i=0; i<resultSize; ++i) delete result[i];
533 result.clear();
534 }
535 G4double D=std::sqrt(R2*G4UniformRand());
536#ifdef pdebug
537 G4cout<<"G4QLowEn::PSDI: D="<<D<<", tR="<<tR<<", pR="<<pR<<G4endl;
538#endif
539 if(D > std::fabs(tR-pR)) // leading parts can be separated
540 {
541 nSec = 0;
542 G4double DtR=D-tR;
543 G4double DpR=D-pR;
544 G4double DtR2=DtR*DtR;
545 G4double DpR2=DpR*DpR;
546 //G4double DtR3=DtR2*DtR;
547 G4double DpR3=DpR2*DpR;
548 //G4double DtR4=DtR3*DtR;
549 G4double DpR4=DpR3*DpR;
550 G4double tR2=tR*tR;
551 G4double pR2=pR*pR;
552 G4double tR3=tR2*tR;
553 //G4double pR3=pR2*pR;
554 G4double tR4=tR3*tR;
555 //G4double pR4=pR3*pR;
556 G4double HD=16.*D;
557 G4double tF=tA*(6*tR2*(pR2-DtR2)-4*D*(tR3-DpR3)+3*(tR4-DpR4))/HD/tR3;
558 G4double pF=tF;
559 tD=static_cast<G4int>(tF);
560 pD=static_cast<G4int>(pF);
561 //if(G4UniformRand() < tF-tD) ++tD; // Simple solution
562 //if(G4UniformRand() < pF-pD) ++pD;
563 // Enhance alphas solution
564 if(std::fabs(tF-4.) < 1.) tD=4; // @@ 1. is the enhancement parameter
565 else if(G4UniformRand() < tF-tD) ++tD;
566 if(std::fabs(pF-4.) < 1.) pD=4;
567 else if(G4UniformRand() < pF-pD) ++pD;
568 if(tD > tB) tD=tB;
569 if(pD > pB) pD=tB;
570 // @@ Quasi Free is not debugged @@ The following close it
571 if(tD < 1) tD=1;
572 if(pD < 1) pD=1;
573#ifdef pdebug
574 G4cout<<"G4QLE::PSDI:#"<<nAt<<",pF="<<pF<<",tF="<<tF<<",pD="<<pD<<",tD="<<tD<<G4endl;
575#endif
576 G4int pC=0; // charge of the projectile fraction
577 G4int tC=0; // charge of the target fraction
578 if((tD || pD) && tD<tB && pD<pB)// Periferal interaction
579 {
580 if(!tD || !pD) // Quasi-Elastic: B+A->(B-1)+N+A || ->B+N+(A-1)
581 {
582 G4VQCrossSection* PCSmanager=G4QProtonElasticCrossSection::GetPointer();
583 G4VQCrossSection* NCSmanager=G4QNeutronElasticCrossSection::GetPointer();
584 G4int pPDG=2112; // Proto of the nucleon PDG (proton)
585 G4double prM =mNeut; // Proto of the nucleon mass
586 G4double prM2=mNeu2; // Proto of the nucleon sq mass
587 if (!tD) // Quasi-elastic scattering of the proj QF nucleon
588 {
589 if(pD > 1) pD=1;
590 if(!pN || (pZ && pA*G4UniformRand() < pZ) ) // proj QF proton
591 {
592 pPDG=2212;
593 prM=mProt;
594 prM2=mPro2;
595 pC=1;
596 }
597 G4double Mom=0.;
598 G4LorentzVector com4M = targ4M; // Proto of 4mom of compound
599 G4double tgM=targ4M.e();
600 G4double rNM=0.;
601 G4LorentzVector rNuc4M(0.,0.,0.,0.);
602 if(pD==pB)
603 {
604 Mom=proj4M.rho();
605 com4M += proj4M; // Total 4-mom for scattering
606 rNM=prM;
607 }
608 else // It is necessary to split the nucleon (with fermiM)
609 {
610 G4ThreeVector fm=mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
611 rNM=G4QPDGCode(2112).GetNuclMass(pZ-pC, pN, 0);
612 G4double rNE=std::sqrt(fm*fm+rNM*rNM);
613 rNuc4M=G4LorentzVector(fm,rNE);
614 G4ThreeVector boostV=proj4M.vect()/proj4M.e();
615 rNuc4M.boost(boostV);
616 G4LorentzVector qfN4M=proj4M - rNuc4M;// 4Mom of the quasi-free nucleon in LS
617 com4M += qfN4M; // Calculate Total 4Mom for NA scattering
618 G4double pNE = qfN4M.e(); // Energy of the QF nucleon in LS
619 if(rNE >= prM) Mom = std::sqrt(pNE*pNE-prM2); // Mom(s) fake value
620 else break; // Break the while loop
621 }
622 xSec=0.;
623 if(pPDG==2212) xSec=PCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG);
624 else xSec=NCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG);
625 if( xSec <= 0. ) break; // Break the while loop
626 G4double mint=0.; // Prototype of functional randomized -t
627 if(pPDG==2212) mint=PCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t
628 else mint=NCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t
629 G4double maxt=0.; // Prototype of maximum -t
630 if(pPDG==2212) maxt=PCSmanager->GetHMaxT(); // maximum -t
631 else maxt=NCSmanager->GetHMaxT(); // maximum -t
632 G4double cost=1.-mint/maxt; // cos(theta) in CMS
633 if(cost>1. || cost<-1.) break; // Break the while loop
634 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,tgM); // 4mom of recoil target
635 G4LorentzVector scat4M=G4LorentzVector(0.,0.,0.,rNM); // 4mom of scattered N
636 G4LorentzVector dir4M=tt4M-G4LorentzVector(0.,0.,0.,(com4M.e()-rNM-prM)*.01);
637 if(!G4QHadron(com4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
638 {
639 G4cout<<"G4QLE::Pt="<<com4M.m()<<",p="<<prM<<",r="<<rNM<<",c="<<cost<<G4endl;
640 break; // Break the while loop
641 }
642 G4Track* projSpect = 0;
643 G4Track* aNucTrack = 0;
644 if(pB > pD) // Fill the proj residual nucleus
645 {
646 G4int rZ=pZ-pC;
647 G4int rA=pB-1;
648 G4ParticleDefinition* theDefinition;// Prototype of residualNucleusDefinition
649 if(rA==1)
650 {
651 if(!rZ) theDefinition = aNeutron;
652 else theDefinition = aProton;
653 }
654 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
655 G4DynamicParticle* resN = new G4DynamicParticle(theDefinition, rNuc4M);
656 projSpect = new G4Track(resN, localtime, position );
657 projSpect->SetWeight(weight); // weighted
658 projSpect->SetTouchableHandle(trTouchable);
659#ifdef pdebug
660 G4cout<<"G4QLowEn::PSDI:-->ProjQFSA="<<rA<<",rZ="<<rZ<<",4M="<<rNuc4M<<G4endl;
661#endif
662 ++nSec;
663 }
664 G4ParticleDefinition* theDefinition = G4Neutron::Neutron();
665 if(pPDG==2212) theDefinition = G4Proton::Proton();
666 G4DynamicParticle* scatN = new G4DynamicParticle(theDefinition, scat4M);
667 aNucTrack = new G4Track(scatN, localtime, position );
668 aNucTrack->SetWeight(weight); // weighted
669 aNucTrack->SetTouchableHandle(trTouchable);
670#ifdef pdebug
671 G4cout<<"G4QLowEn::PSDI:-->TgQFNPDG="<<pPDG<<", 4M="<<scat4M<<G4endl;
672#endif
673 ++nSec;
674 G4Track* aFraTrack=0;
675 if(tA==1)
676 {
677 if(!tZ) theDefinition = aNeutron;
678 else theDefinition = aProton;
679 }
680 else if(tA==8 && tC==4) // Be8 decay
681 {
682 theDefinition = anAlpha;
683 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha
684 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha
685 if(!G4QHadron(reco4M).DecayIn2(f4M, s4M))
686 {
687 G4cout<<"*G4QLE::TS->A+A,t="<<reco4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
688 }
689 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
690 aFraTrack = new G4Track(pAl, localtime, position );
691 aFraTrack->SetWeight(weight); // weighted
692 aFraTrack->SetTouchableHandle(trTouchable);
693#ifdef pdebug
694 G4cout<<"G4QLowEn::PSDI:-->TgRQFA4M="<<f4M<<G4endl;
695#endif
696 ++nSec;
697 reco4M=s4M;
698 }
699 else if(tA==5 && (tC==2 || tC==3)) // He5/Li5 decay
700 {
701 theDefinition = aProton;
702 G4double mNuc = mProt;
703 if(tC==2)
704 {
705 theDefinition = aNeutron;
706 mNuc = mNeut;
707 }
708 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon
709 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha
710 if(!G4QHadron(reco4M).DecayIn2(f4M, s4M))
711 {
712 G4cout<<"*G4QLE::TS->N+A,t="<<reco4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
713 }
714 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
715 aFraTrack = new G4Track(pAl, localtime, position );
716 aFraTrack->SetWeight(weight); // weighted
717 aFraTrack->SetTouchableHandle(trTouchable);
718#ifdef pdebug
719 G4cout<<"G4QLowEn::PSDI:-->TgQFRN4M="<<f4M<<G4endl;
720#endif
721 ++nSec;
722 theDefinition = anAlpha;
723 reco4M=s4M;
724 }
725 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(tZ,tB,0,tZ);
726 ++nSec;
727#ifdef pdebug
728 G4cout<<"G4QLowEn::PSDI:-->PQF_nSec="<<nSec<<G4endl;
729#endif
730 aParticleChange.SetNumberOfSecondaries(nSec);
731 if(projSpect) aParticleChange.AddSecondary(projSpect);
732 if(aNucTrack) aParticleChange.AddSecondary(aNucTrack);
733 if(aFraTrack) aParticleChange.AddSecondary(aFraTrack);
734 G4DynamicParticle* resA = new G4DynamicParticle(theDefinition, reco4M);
735 G4Track* aResTrack = new G4Track(resA, localtime, position );
736 aResTrack->SetWeight(weight); // weighted
737 aResTrack->SetTouchableHandle(trTouchable);
738#ifdef pdebug
739 G4cout<<"G4QLowEn::PSDI:-->TgR4M="<<reco4M<<", checkNSec="<<nSec<<G4endl;
740#endif
741 aParticleChange.AddSecondary(aResTrack);
742 }
743 else // !pD : QF target Nucleon on the whole Projectile
744 {
745 if(tD > 1) tD=1;
746 if(!tN || (tZ && tA*G4UniformRand() < tZ) ) // target QF proton
747 {
748 pPDG=2212;
749 prM=mProt;
750 prM2=mPro2;
751 tC=1;
752 }
753 G4double Mom=0.;
754 G4LorentzVector com4M=proj4M; // Proto of 4mom of compound
755 prM=proj4M.m();
756 G4double rNM=0.;
757 G4LorentzVector rNuc4M(0.,0.,0.,0.);
758 if(tD==tB)
759 {
760 Mom=prM*proj4M.rho()/proj4M.m();
761 com4M += targ4M; // Total 4-mom for scattering
762 rNM=prM;
763 }
764 else // It is necessary to split the nucleon (with fermiM)
765 {
766 G4ThreeVector fm=250.*std::pow(G4UniformRand(),third)*G4RandomDirection();
767 rNM=G4QPDGCode(2112).GetNuclMass(tZ-tC, tN, 0)/MeV;
768 G4double rNE=std::sqrt(fm*fm+rNM*rNM);
769 rNuc4M=G4LorentzVector(fm,rNE);
770 G4LorentzVector qfN4M=targ4M - rNuc4M;// 4Mom of the quasi-free nucleon in LS
771 com4M += qfN4M; // Calculate Total 4Mom for NA scattering
772 G4ThreeVector boostV=proj4M.vect()/proj4M.e();
773 qfN4M.boost(-boostV);
774 G4double tNE = qfN4M.e(); // Energy of the QF nucleon in LS
775 if(rNE >= prM) Mom = std::sqrt(tNE*tNE-prM2); // Mom(s) fake value
776 else break; // Break the while loop
777 }
778 xSec=0.;
779 if(pPDG==2212) xSec=PCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG);
780 else xSec=NCSmanager->GetCrossSection(false, Mom, tZ, tN, pPDG);
781 if( xSec <= 0. ) break; // Break the while loop
782 G4double mint=0.; // Prototype of functional randomized -t
783 if(pPDG==2212) mint=PCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t
784 else mint=NCSmanager->GetExchangeT(tZ,tN,pPDG); // randomized -t
785 G4double maxt=0.; // Prototype of maximum -t
786 if(pPDG==2212) maxt=PCSmanager->GetHMaxT(); // maximum -t
787 else maxt=NCSmanager->GetHMaxT(); // maximum -t
788 G4double cost=1.-mint/maxt; // cos(theta) in CMS
789 if(cost>1. || cost<-1.) break; // Break the while loop
790 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,prM); // 4mom of recoil target
791 G4LorentzVector scat4M=G4LorentzVector(0.,0.,0.,rNM); // 4mom of scattered N
792 G4LorentzVector dir4M=tt4M-G4LorentzVector(0.,0.,0.,(com4M.e()-rNM-prM)*.01);
793 if(!G4QHadron(com4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
794 {
795 G4cout<<"G4QLE::Tt="<<com4M.m()<<",p="<<prM<<",r="<<rNM<<",c="<<cost<<G4endl;
796 break; // Break the while loop
797 }
798 G4Track* targSpect = 0;
799 G4Track* aNucTrack = 0;
800 if(tB > tD) // Fill the residual nucleus
801 {
802 G4int rZ=tZ-tC;
803 G4int rA=tB-1;
804 G4ParticleDefinition* theDefinition;// Prototype of residualNucleusDefinition
805 if(rA==1)
806 {
807 if(!rZ) theDefinition = aNeutron;
808 else theDefinition = aProton;
809 }
810 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
811 G4DynamicParticle* resN = new G4DynamicParticle(theDefinition, rNuc4M);
812 targSpect = new G4Track(resN, localtime, position );
813 targSpect->SetWeight(weight); // weighted
814 targSpect->SetTouchableHandle(trTouchable);
815#ifdef pdebug
816 G4cout<<"G4QLowEn::PSDI:-->TargQFSA="<<rA<<",rZ="<<rZ<<",4M="<<rNuc4M<<G4endl;
817#endif
818 ++nSec;
819 }
820 G4ParticleDefinition* theDefinition = G4Neutron::Neutron();
821 if(pPDG==2212) theDefinition = G4Proton::Proton();
822 G4DynamicParticle* scatN = new G4DynamicParticle(theDefinition, scat4M);
823 aNucTrack = new G4Track(scatN, localtime, position );
824 aNucTrack->SetWeight(weight); // weighted
825 aNucTrack->SetTouchableHandle(trTouchable);
826#ifdef pdebug
827 G4cout<<"G4QLowEn::PSDI:-->PrQFNPDG="<<pPDG<<", 4M="<<scat4M<<G4endl;
828#endif
829 ++nSec;
830 G4Track* aFraTrack=0;
831 if(pA==1)
832 {
833 if(!pZ) theDefinition = aNeutron;
834 else theDefinition = aProton;
835 }
836 else if(pA==8 && pC==4) // Be8 decay
837 {
838 theDefinition = anAlpha;
839 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha
840 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha
841 if(!G4QHadron(reco4M).DecayIn2(f4M, s4M))
842 {
843 G4cout<<"*G4QLE::PS->A+A,t="<<reco4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
844 }
845 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
846 aFraTrack = new G4Track(pAl, localtime, position );
847 aFraTrack->SetWeight(weight); // weighted
848 aFraTrack->SetTouchableHandle(trTouchable);
849#ifdef pdebug
850 G4cout<<"G4QLowEn::PSDI:-->PrRQFA4M="<<f4M<<G4endl;
851#endif
852 ++nSec;
853 reco4M=s4M;
854 }
855 else if(pA==5 && (pC==2 || pC==3)) // He5/Li5 decay
856 {
857 theDefinition = aProton;
858 G4double mNuc = mProt;
859 if(pC==2)
860 {
861 theDefinition = aNeutron;
862 mNuc = mNeut;
863 }
864 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon
865 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha
866 if(!G4QHadron(reco4M).DecayIn2(f4M, s4M))
867 {
868 G4cout<<"*G4QLE::PS->N+A,t="<<reco4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
869 }
870 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
871 aFraTrack = new G4Track(pAl, localtime, position );
872 aFraTrack->SetWeight(weight); // weighted
873 aFraTrack->SetTouchableHandle(trTouchable);
874#ifdef pdebug
875 G4cout<<"G4QLowEn::PSDI:-->PrQFRN4M="<<f4M<<G4endl;
876#endif
877 ++nSec;
878 theDefinition = anAlpha;
879 reco4M=s4M;
880 }
881 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(pZ,pB,0,pZ);
882 ++nSec;
883#ifdef pdebug
884 G4cout<<"G4QLowEn::PSDI:-->TQF_nSec="<<nSec<<G4endl;
885#endif
886 aParticleChange.SetNumberOfSecondaries(nSec);
887 if(targSpect) aParticleChange.AddSecondary(targSpect);
888 if(aNucTrack) aParticleChange.AddSecondary(aNucTrack);
889 if(aFraTrack) aParticleChange.AddSecondary(aFraTrack);
890 G4DynamicParticle* resA = new G4DynamicParticle(theDefinition, reco4M);
891 G4Track* aResTrack = new G4Track(resA, localtime, position );
892 aResTrack->SetWeight(weight); // weighted
893 aResTrack->SetTouchableHandle(trTouchable);
894#ifdef pdebug
895 G4cout<<"G4QLowEn::PSDI:-->TgR4M="<<reco4M<<", checkNSec="<<nSec<<G4endl;
896#endif
897 aParticleChange.AddSecondary( aResTrack );
898 }
899#ifdef debug
900 G4cout<<"G4QLowEnergy::PostStepDoIt:***PostStepDoIt is done:Quasi-El***"<<G4endl;
901#endif
902 return G4VDiscreteProcess::PostStepDoIt(track, step); // ===> Reaction is DONE
903 }
904 else // The cental region compound can be created
905 {
906 // First calculate the isotopic state of the parts of the compound
907 if (!pZ) pC = 0;
908 else if(!pN) pC = pD;
909 else
910 {
911#ifdef pdebug
912 G4cout<<"G4QLowEn::PSDI: pD="<<pD<<", pZ="<<pZ<<", pA="<<pA<<G4endl;
913#endif
914 G4double C=pD*pZ/pA;
915 pC=static_cast<G4int>(C);
916 if(G4UniformRand() < C-pC) ++pC;
917 }
918 if(!tZ) tC=0;
919 else if(!tN) tC=tD;
920 else
921 {
922#ifdef pdebug
923 G4cout<<"G4QLowEn::PSDI: tD="<<tD<<", tZ="<<tZ<<", tA="<<tA<<G4endl;
924#endif
925 G4double C=tD*tZ/tA;
926 tC=static_cast<G4int>(C);
927 if(G4UniformRand() < C-tC) ++tC;
928 }
929 // calculate the transferred momentum
930 G4ThreeVector pFM(0.,0.,0.);
931 if(pD < pB) // The projectile nucleus must be splitted
932 {
933 G4int nc=pD;
934 if(pD+pD>pB) nc=pB-pD;
935 pFM = mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
936 for(i=1; i < nc; ++i)
937 pFM+= mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
938 }
939 G4ThreeVector tFM(0.,0.,0.);
940 if(tD<tB) // The projectile nucleus must be splitted
941 {
942 G4int nc=pD;
943 if(tD+tD>tB) nc=tB-tD;
944 tFM = mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
945 for(i=1; i < nc; ++i)
946 tFM+= mFm*std::pow(G4UniformRand(),third)*G4RandomDirection();
947 }
948#ifdef pdebug
949 G4cout<<"G4QLE::PSDI:pC="<<pC<<", tC="<<tC<<", pFM="<<pFM<<", tFM="<<tFM<<G4endl;
950#endif
951 // Split the projectile spectator
952 G4int rpZ=pZ-pC; // number of protons in the projectile spectator
953 G4int pF_value=pD-pC; // number of neutrons in the projectile part of comp
954 G4int rpN=pN-pF_value; // number of neutrons in the projectile spectator
955 G4double rpNM=G4QPDGCode(2112).GetNuclMass(rpZ, rpN, 0); // Mass of the spectator
956 G4ThreeVector boostV=proj4M.vect()/proj4M.e(); // Antilab Boost Vector
957 G4double rpE=std::sqrt(rpNM*rpNM+pFM.mag2());
958 G4LorentzVector rp4M(pFM,rpE);
959#ifdef pdebug
960 G4cout<<"G4QLE::PSDI: boostV="<<boostV<<",rp4M="<<rp4M<<",pr4M="<<proj4M<<G4endl;
961#endif
962 rp4M.boost(boostV);
963#ifdef pdebug
964 G4cout<<"G4QLE::PSDI: After boost, rp4M="<<rp4M<<G4endl;
965#endif
966 G4ParticleDefinition* theDefinition; // Prototype of projSpectatorNuclDefinition
967 G4int rpA=rpZ+rpN;
968 G4Track* aFraPTrack = 0;
969 theDefinition = 0;
970 if(rpA==1)
971 {
972 if(!rpZ) theDefinition = G4Neutron::Neutron();
973 else theDefinition = G4Proton::Proton();
974#ifdef pdebug
975 G4cout<<"G4QLE::PSDI: rpA=1, rpZ"<<rpZ<<G4endl;
976#endif
977 }
978 else if(rpA==2 && rpZ==0) // nn decay
979 {
980 theDefinition = aNeutron;
981 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 1st neutron
982 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 2nd neutron
983#ifdef pdebug
984 G4cout<<"G4QLE::CPS->n+n,nn="<<rp4M.m()<<" >? 2*MNeutron="<<2*mNeutron<<G4endl;
985#endif
986 if(!G4QHadron(rp4M).DecayIn2(f4M, s4M))
987 {
988 G4cout<<"*W*G4QLE::CPS->n+n,t="<<rp4M.m()<<" >? 2*Neutron="<<2*mAlph<<G4endl;
989 }
990 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
991 aFraPTrack = new G4Track(pNeu, localtime, position );
992 aFraPTrack->SetWeight(weight); // weighted
993 aFraPTrack->SetTouchableHandle(trTouchable);
994 tt4M-=f4M;
995#ifdef edebug
996 totBaN-=2;
997 tch4M -=f4M;
998 G4cout<<">>G4QLEn::PSDI:n,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
999#endif
1000#ifdef pdebug
1001 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
1002#endif
1003 rp4M=s4M;
1004 }
1005 else if(rpA>2 && rpZ==0) // Z=0 decay
1006 {
1007 theDefinition = aNeutron;
1008 G4LorentzVector f4M=rp4M/rpA; // 4mom of 1st neutron
1009#ifdef pdebug
1010 G4cout<<"G4QLE::CPS->Nn,M="<<rp4M.m()<<" >? N*MNeutron="<<rpA*mNeutron<<G4endl;
1011#endif
1012 for(G4int it=1; it<rpA; ++it) // Fill (N-1) neutrons to output
1013 {
1014 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
1015 G4Track* aNTrack = new G4Track(pNeu, localtime, position );
1016 aNTrack->SetWeight(weight); // weighted
1017 aNTrack->SetTouchableHandle(trTouchable);
1018 result.push_back(aNTrack);
1019 }
1020 G4int nesc = rpA-1;
1021 tt4M-=f4M*nesc;
1022#ifdef edebug
1023 totBaN-=nesc;
1024 tch4M -=f4M*nesc;
1025 G4cout<<">G4QLEn::PSDI:Nn,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1026#endif
1027#ifdef pdebug
1028 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
1029#endif
1030 rp4M=f4M;
1031 }
1032 else if(rpA==8 && rpZ==4) // Be8 decay
1033 {
1034 theDefinition = anAlpha;
1035 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha
1036 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha
1037#ifdef pdebug
1038 G4cout<<"G4QLE::CPS->A+A,mBe8="<<rp4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
1039#endif
1040 if(!G4QHadron(rp4M).DecayIn2(f4M, s4M))
1041 {
1042 G4cout<<"*W*G4QLE::CPS->A+A,t="<<rp4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
1043 }
1044 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
1045 aFraPTrack = new G4Track(pAl, localtime, position );
1046 aFraPTrack->SetWeight(weight); // weighted
1047 aFraPTrack->SetTouchableHandle(trTouchable);
1048 tt4M-=f4M;
1049#ifdef edebug
1050 totChg-=2;
1051 totBaN-=4;
1052 tch4M -=f4M;
1053 G4cout<<">>G4QLEn::PSDI:1,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1054#endif
1055#ifdef pdebug
1056 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
1057#endif
1058 rp4M=s4M;
1059 }
1060 else if(rpA==5 && (rpZ==2 || rpZ==3)) // He5/Li5 decay
1061 {
1062 theDefinition = aProton;
1063 G4double mNuc = mProt;
1064 if(rpZ==2)
1065 {
1066 theDefinition = aNeutron;
1067 mNuc = mNeut;
1068 }
1069 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon
1070 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha
1071#ifdef pdebug
1072 G4cout<<"G4QLowE::CPS->N+A, tM5="<<rp4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
1073#endif
1074 if(!G4QHadron(rp4M).DecayIn2(f4M, s4M))
1075 {
1076 G4cout<<"*W*G4QLE::CPS->N+A,t="<<rp4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
1077 }
1078 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
1079 aFraPTrack = new G4Track(pAl, localtime, position );
1080 aFraPTrack->SetWeight(weight); // weighted
1081 aFraPTrack->SetTouchableHandle(trTouchable);
1082 tt4M-=f4M;
1083#ifdef edebug
1084 if(theDefinition == aProton) totChg-=1;
1085 totBaN-=1;
1086 tch4M -=f4M;
1087 G4cout<<">>G4QLEn::PSDI:2,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1088#endif
1089#ifdef pdebug
1090 G4cout<<"G4QLowEn::PSDI:-->ProjSpectN4M="<<f4M<<G4endl;
1091#endif
1092 theDefinition = anAlpha;
1093 rp4M=s4M;
1094 }
1095 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rpZ,rpA,0,rpZ);
1096 if(!theDefinition)
1097 {
1098 // G4cout<<"-Warning-G4QLowEn::PSDI: pDef=0, Z="<<rpZ<<", A="<<rpA<<G4endl;
1099 // throw G4QException("G4QLowE::PoStDoIt: particle definition is a null pointer");
1100 G4ExceptionDescription ed;
1101 ed << "Particle definition is a null pointer: pDef=0, Z= " << rpZ
1102 << ", A=" << rpA << G4endl;
1103 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
1104 JustWarning, ed);
1105 }
1106#ifdef edebug
1107 if(theDefinition == anAlpha)
1108 {
1109 totChg-=2;
1110 totBaN-=4;
1111 }
1112 else
1113 {
1114 totChg-=rpZ;
1115 totBaN-=rpA;
1116 }
1117 tch4M -=rp4M;
1118 G4cout<<">>G4QLEn::PSDI:3, tZ="<<totChg<<",tB="<<totBaN<<", t4M="<<tch4M<<G4endl;
1119#endif
1120 G4DynamicParticle* rpD = new G4DynamicParticle(theDefinition, rp4M);
1121 G4Track* aNewPTrack = new G4Track(rpD, localtime, position);
1122 aNewPTrack->SetWeight(weight);// weighted
1123 aNewPTrack->SetTouchableHandle(trTouchable);
1124 tt4M-=rp4M;
1125#ifdef pdebug
1126 G4cout<<"G4QLowEn::PSDI:-->ProjSpectR4M="<<rp4M<<",Z="<<rpZ<<",A="<<rpA<<G4endl;
1127#endif
1128 //
1129 // Split the target spectator
1130 G4int rtZ=tZ-tC; // number of protons in the target spectator
1131 G4int tF_value=tD-tC; // number of neutrons in the target part of comp
1132 G4int rtN=tN-tF_value; // number of neutrons in the target spectator
1133 G4double rtNM=G4QPDGCode(2112).GetNuclMass(rtZ, rtN, 0); // Mass of the spectator
1134 G4double rtE=std::sqrt(rtNM*rtNM+tFM.mag2());
1135 G4LorentzVector rt4M(tFM,rtE);
1136 G4int rtA=rtZ+rtN;
1137 G4Track* aFraTTrack = 0;
1138 theDefinition = 0;
1139 if(rtA==1)
1140 {
1141 if(!rtZ) theDefinition = G4Neutron::Neutron();
1142 else theDefinition = G4Proton::Proton();
1143#ifdef pdebug
1144 G4cout<<"G4QLE::PSDI: rtA=1, rtZ"<<rtZ<<G4endl;
1145#endif
1146 }
1147 else if(rtA==2 && rtZ==0) // nn decay
1148 {
1149 theDefinition = aNeutron;
1150 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 1st neutron
1151 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mNeut); // 4mom of 2nd neutron
1152#ifdef pdebug
1153 G4cout<<"G4QLE::CPS->n+n,nn="<<rptM.m()<<" >? 2*MNeutron="<<2*mNeutron<<G4endl;
1154#endif
1155 if(!G4QHadron(rt4M).DecayIn2(f4M, s4M))
1156 G4cout<<"*W*G4QLE::CPS->n+n,t="<<rt4M.m()<<" >? 2*Neutron="<<2*mAlph<<G4endl;
1157 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
1158 aFraPTrack = new G4Track(pNeu, localtime, position );
1159 aFraPTrack->SetWeight(weight); // weighted
1160 aFraPTrack->SetTouchableHandle(trTouchable);
1161 tt4M-=f4M;
1162#ifdef edebug
1163 totBaN-=2;
1164 tch4M -=f4M;
1165 G4cout<<">>G4QLEn::PSDI:n,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1166#endif
1167#ifdef pdebug
1168 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
1169#endif
1170 rt4M=s4M;
1171 }
1172 else if(rtA>2 && rtZ==0) // Z=0 decay
1173 {
1174 theDefinition = aNeutron;
1175 G4LorentzVector f4M=rt4M/rtA; // 4mom of 1st neutron
1176#ifdef pdebug
1177 G4cout<<"G4QLE::CPS->Nn,M="<<rt4M.m()<<" >? N*MNeutron="<<rtA*mNeutron<<G4endl;
1178#endif
1179 for(G4int it=1; it<rtA; ++it) // Fill (N-1) neutrons to output
1180 {
1181 G4DynamicParticle* pNeu = new G4DynamicParticle(theDefinition, f4M);
1182 G4Track* aNTrack = new G4Track(pNeu, localtime, position );
1183 aNTrack->SetWeight(weight); // weighted
1184 aNTrack->SetTouchableHandle(trTouchable);
1185 result.push_back(aNTrack);
1186 }
1187 G4int nesc = rtA-1;
1188 tt4M-=f4M*nesc;
1189#ifdef edebug
1190 totBaN-=nesc;
1191 tch4M -=f4M*nesc;
1192 G4cout<<">G4QLEn::PSDI:Nn,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1193#endif
1194#ifdef pdebug
1195 G4cout<<"G4QLowEn::PSDI:-->ProjSpectA4M="<<f4M<<G4endl;
1196#endif
1197 rt4M=f4M;
1198 }
1199 else if(rtA==8 && rtZ==4) // Be8 decay
1200 {
1201 theDefinition = anAlpha;
1202 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 1st alpha
1203 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of 2nd alpha
1204#ifdef pdebug
1205 G4cout<<"G4QLE::CPS->A+A,mBe8="<<rt4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
1206#endif
1207 if(!G4QHadron(rt4M).DecayIn2(f4M, s4M))
1208 {
1209 G4cout<<"*W*G4QLE::CTS->A+A,t="<<rt4M.m()<<" >? 2*MAlpha="<<2*mAlph<<G4endl;
1210 }
1211 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
1212 aFraTTrack = new G4Track(pAl, localtime, position );
1213 aFraTTrack->SetWeight(weight); // weighted
1214 aFraTTrack->SetTouchableHandle(trTouchable);
1215 tt4M-=f4M;
1216#ifdef edebug
1217 totChg-=2;
1218 totBaN-=4;
1219 tch4M -=f4M;
1220 G4cout<<">>G4QLEn::PSDI:4,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1221#endif
1222#ifdef pdebug
1223 G4cout<<"G4QLowEn::PSDI:-->TargSpectA4M="<<f4M<<G4endl;
1224#endif
1225 rt4M=s4M;
1226 }
1227 else if(rtA==5 && (rtZ==2 || rtZ==3)) // He5/Li5 decay
1228 {
1229 theDefinition = aProton;
1230 G4double mNuc = mProt;
1231 if(rpZ==2)
1232 {
1233 theDefinition = aNeutron;
1234 mNuc = mNeut;
1235 }
1236 G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,mNuc); // 4mom of the nucleon
1237 G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,mAlph); // 4mom of the alpha
1238#ifdef pdebug
1239 G4cout<<"G4QLowE::CPS->N+A, tM5="<<rt4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
1240#endif
1241 if(!G4QHadron(rt4M).DecayIn2(f4M, s4M))
1242 {
1243 G4cout<<"*W*G4QLE::CTS->N+A,t="<<rt4M.m()<<" >? MA+MN="<<mAlph+mNuc<<G4endl;
1244 }
1245 G4DynamicParticle* pAl = new G4DynamicParticle(theDefinition, f4M);
1246 aFraTTrack = new G4Track(pAl, localtime, position );
1247 aFraTTrack->SetWeight(weight); // weighted
1248 aFraTTrack->SetTouchableHandle(trTouchable);
1249 tt4M-=f4M;
1250#ifdef edebug
1251 if(theDefinition == aProton) totChg-=1;
1252 totBaN-=1;
1253 tch4M -=f4M;
1254 G4cout<<">>G4QLEn::PSDI:5,tZ="<<totChg<<",tB="<<totBaN<<",t4M="<<tch4M<<G4endl;
1255#endif
1256#ifdef pdebug
1257 G4cout<<"G4QLowEn::PSDI:-->TargSpectN4M="<<f4M<<G4endl;
1258#endif
1259 theDefinition = anAlpha;
1260 rt4M=s4M;
1261 }
1262 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rtZ,rtA,0,rtZ);
1263 if(!theDefinition)
1264 {
1265 // G4cout<<"-Warning-G4QLowEn::PSDI: tDef=0, Z="<<rtZ<<", A="<<rtA<<G4endl;
1266 // throw G4QException("G4QLowE::PoStDoIt: particle definition is a null pointer");
1267 G4ExceptionDescription ed;
1268 ed << "Particle definition is a null pointer: tDef=0, Z= " << rtZ
1269 << ", A=" << rtA << G4endl;
1270 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
1271 JustWarning, ed);
1272 }
1273#ifdef edebug
1274 if(theDefinition == anAlpha)
1275 {
1276 totChg-=2;
1277 totBaN-=4;
1278 }
1279 else
1280 {
1281 totChg-=rtZ;
1282 totBaN-=rtA;
1283 }
1284 tch4M -=rt4M;
1285 G4cout<<">>G4QLEn::PSDI:6, tZ="<<totChg<<",tB="<<totBaN<<", t4M="<<tch4M<<G4endl;
1286#endif
1287 G4DynamicParticle* rtD = new G4DynamicParticle(theDefinition, rt4M);
1288 G4Track* aNewTTrack = new G4Track(rtD, localtime, position);
1289 aNewTTrack->SetWeight(weight); // weighted
1290 aNewTTrack->SetTouchableHandle(trTouchable);
1291 tt4M-=rt4M;
1292#ifdef pdebug
1293 G4cout<<"G4QLowEn::PSDI:-->TargSpectR4M="<<rt4M<<",Z="<<rtZ<<",A="<<rtA<<G4endl;
1294#endif
1295 if(aFraPTrack) result.push_back(aFraPTrack);
1296 if(aNewPTrack) result.push_back(aNewPTrack);
1297 if(aFraTTrack) result.push_back(aFraTTrack);
1298 if(aNewTTrack) result.push_back(aNewTTrack);
1299 tcM = tt4M.m(); // total CMS mass of the compound (after reduction!)
1300 G4int sN=tF_value+pF_value;
1301 G4int sZ=tC+pC;
1302 tnM = targQPDG.GetNuclMass(sZ,sN,0); // GSM
1303#ifdef pdebug
1304 G4cout<<"G4QLEn::PSDI:At#"<<nAt<<",totM="<<tcM<<",gsM="<<tnM<<",Z="<<sZ<<",N="
1305 <<sN<<G4endl;
1306#endif
1307 if(tcM > tnM) // !! The totalresidual reaction is modified !
1308 {
1309 pZ=pC;
1310 pN=pF_value;
1311 tZ=tC;
1312 tN=tF_value;
1313 tot4M=tt4M;
1314 }
1315 //else
1316 //{
1317 // G4cout<<"***G4QLowEnergy::PostStepDoIt: totM="<<tcM<<" <= GSM="<<tnM<<G4endl;
1318 // throw G4QException("G4QLowEnergy::PostStepDoIt: ResibualNucl below GSM shell");
1319 //}
1320 }
1321 } // At least one is splitted
1322 else if(tD==tB || pD==pB) // Total absorption
1323 {
1324 tD=tZ+tN;
1325 pD=pZ+pN;
1326 tcM=tnM+1.;
1327 }
1328 }
1329 else // Total compound (define tD to go out of the while)
1330 {
1331 tD=tZ+tN;
1332 pD=pZ+pN;
1333 tcM=tnM+1.;
1334 }
1335 } // End of the interaction WHILE
1336 G4double totM=tot4M.m(); // total CMS mass of the reaction
1337 G4int totN=tN+pN;
1338 G4int totZ=tZ+pZ;
1339#ifdef edebug
1340 G4cout<<">>>G4QLEn::PSDI: dZ="<<totChg-totZ<<", dB="<<totBaN-totN-totZ<<", d4M="
1341 <<tch4M-tot4M<<",N="<<totN<<",Z="<<totZ<<G4endl;
1342#endif
1343 // @@ Here mass[i] can be calculated if mass=0
1344 G4double mass[nCh]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
1345 0.,0.,0.,0.,0.,0.};
1346 mass[0] = tM = targQPDG.GetNuclMass(totZ, totN, 0); // gamma+gamma and GSM update
1347#ifdef pdebug
1348 G4cout<<"G4QLEn::PSDI:TotM="<<totM<<",NucM="<<tM<<",totZ="<<totZ<<",totN="<<totN<<G4endl;
1349#endif
1350 if (totN>0 && totZ>0)
1351 {
1352 mass[1] = targQPDG.GetNuclMass(totZ,totN-1,0); // gamma+neutron
1353 mass[2] = targQPDG.GetNuclMass(totZ-1,totN,0); // gamma+proton
1354 }
1355 if ( totZ > 1 && totN > 1 ) mass[3] = targQPDG.GetNuclMass(totZ-1,totN-1,0); //g+d
1356 if ( totZ > 1 && totN > 2 ) mass[4] = targQPDG.GetNuclMass(totZ-1,totN-2,0); //g+t
1357 if ( totZ > 2 && totN > 1 ) mass[5] = targQPDG.GetNuclMass(totZ-2,totN-1,0); //g+3
1358 if ( totZ > 2 && totN > 2 ) mass[6] = targQPDG.GetNuclMass(totZ-2,totN-2,0); //g+a
1359 if ( totZ > 0 && totN > 2 ) mass[7] = targQPDG.GetNuclMass(totZ ,totN-2,0); //n+n
1360 mass[ 8] = mass[3]; // neutron+proton (the same as a deuteron)
1361 if ( totZ > 2 ) mass[9] = targQPDG.GetNuclMass(totZ-2,totN ,0); //p+p
1362 mass[10] = mass[5]; // proton+deuteron (the same as He3)
1363 mass[11] = mass[4]; // neutron+deuteron (the same as t)
1364 mass[12] = mass[6]; // deuteron+deuteron (the same as alpha)
1365 mass[13] = mass[6]; // proton+tritium (the same as alpha)
1366 if ( totZ > 1 && totN > 3 ) mass[14] = targQPDG.GetNuclMass(totZ-1,totN-3,0);//n+t
1367 if ( totZ > 3 && totN > 1 ) mass[15] = targQPDG.GetNuclMass(totZ-3,totN-1,0);//He3+p
1368 mass[16] = mass[6]; // neutron+He3 (the same as alpha)
1369 if ( totZ > 3 && totN > 2 ) mass[17] = targQPDG.GetNuclMass(totZ-3,totN-2,0);//pa
1370 if ( totZ > 2 && totN > 3 ) mass[18] = targQPDG.GetNuclMass(totZ-2,totN-3,0);//na
1371 if(pL>0) // @@ Not debugged @@
1372 {
1373 G4int pL1=pL-1;
1374 if(totN>0||totZ>0) mass[19] = targQPDG.GetNuclMass(totZ ,totN ,pL1);// Lambda+gamma
1375 if( (totN > 0 && totZ > 0) || totZ > 1 )
1376 mass[20]=targQPDG.GetNuclMass(totZ-1,totN ,pL1);//Lp
1377 if( (totN > 0 && totZ > 0) || totN > 0 )
1378 mass[21]=targQPDG.GetNuclMass(totZ ,totN-1,pL1);//Ln
1379 if( (totN > 1 && totZ > 0) || (totN > 0 && totZ > 1) )
1380 mass[22]=targQPDG.GetNuclMass(totZ-1,totN-1,pL1);//Ld
1381 if( (totN > 2 && totZ > 0) || (totN > 1 && totZ > 1) )
1382 mass[23]=targQPDG.GetNuclMass(totZ-1,totN-2,pL1);//Lt
1383 if( (totN > 0 && totZ > 2) || (totN > 1 && totZ > 1) )
1384 mass[24]=targQPDG.GetNuclMass(totZ-2,totN-1,pL1);//L3
1385 if( (totN > 1 && totZ > 2) || (totN > 2 && totZ > 1) )
1386 mass[25]=targQPDG.GetNuclMass(totZ-2,totN-2,pL1);//La
1387 }
1388#ifdef debug
1389 G4cout<<"G4QLowEn::PSDI: Residual masses are calculated"<<G4endl;
1390#endif
1391 tA=tZ+tN;
1392 pA=pZ+pN;
1393 G4double prZ=pZ/pA+tZ/tA;
1394 G4double prN=pN/pA+tN/tA;
1395 G4double prD=prN*prZ;
1396 G4double prA=prD*prD;
1397 G4double prH=prD*prZ;
1398 G4double prT=prD*prN;
1399 G4double fhe3=6.*std::sqrt(tA);
1400 G4double prL=0.;
1401 if(pL>0) prL=pL/pA;
1402 G4double qval[nCh]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
1403 0.,0.,0.,0.,0.,0.};
1404 qval[ 0] = (totM - mass[ 0])/137./137.;
1405 qval[ 1] = (totM - mass[ 1] - mNeut)*prN/137.;
1406 qval[ 2] = (totM - mass[ 2] - mProt)*prZ/137.;
1407 qval[ 3] = (totM - mass[ 3] - mDeut)*prD/3./137.;
1408 qval[ 4] = (totM - mass[ 4] - mTrit)*prT/6./137.;
1409 qval[ 5] = (totM - mass[ 5] - mHel3)*prH/fhe3/137.;
1410 qval[ 6] = (totM - mass[ 6] - mAlph)*prA/9./137.;
1411 qval[ 7] = (totM - mass[ 7] - mNeut - mNeut)*prN*prN;
1412 qval[ 8] = (totM - mass[ 8] - mNeut - mProt)*prD;
1413 qval[ 9] = (totM - mass[ 9] - mProt - mProt)*prZ*prZ;
1414 qval[10] = (totM - mass[10] - mProt - mDeut)*prH/3.;
1415 qval[11] = (totM - mass[11] - mNeut - mDeut)*prT/3.;
1416 qval[12] = (totM - mass[12] - mDeut - mDeut)*prA/3./3.;
1417 qval[13] = (totM - mass[13] - mProt - mTrit)*prA/6.;
1418 qval[14] = (totM - mass[14] - mNeut - mTrit)*prT*prN/6.;
1419 qval[15] = (totM - mass[15] - mProt - mHel3)*prH*prZ/fhe3;
1420 qval[16] = (totM - mass[16] - mNeut - mHel3)*prA/fhe3;
1421 qval[17] = (totM - mass[17] - mProt - mAlph)*prZ*prA/9.;
1422 qval[18] = (totM - mass[18] - mNeut - mAlph)*prN*prA/9.;
1423 if(pZ>0)
1424 {
1425 qval[19] = (totM - mass[19] - mLamb)*prL;
1426 qval[20] = (totM - mass[20] - mProt - mLamb)*prL*prZ;
1427 qval[21] = (totM - mass[21] - mNeut - mLamb)*prL*prN;
1428 qval[22] = (totM - mass[22] - mDeut - mLamb)*prL*prD/2.;
1429 qval[23] = (totM - mass[23] - mTrit - mLamb)*prL*prT/3.;
1430 qval[24] = (totM - mass[24] - mHel3 - mLamb)*prL*prH/fhe3;
1431 qval[25] = (totM - mass[25] - mAlph - mLamb)*prL*prA/4;
1432 }
1433#ifdef debug
1434 G4cout<<"G4QLowEn::PSDI: Q-values are calculated, tgA="<<tA<<", prA="<<pA<<G4endl;
1435#endif
1436
1437 G4double qv = 0.0; // Total sum of probabilities (q-values)
1438 for(i=0; i<nCh; ++i )
1439 {
1440#ifdef sdebug
1441 G4cout<<"G4QLowEn::PSDI: i="<<i<<", q="<<qval[i]<<",m="<<mass[i]<<G4endl;
1442#endif
1443 if( mass[i] < 500.*MeV ) qval[i] = 0.0; // Close A/Z impossible channels (mesons)
1444 if( qval[i] < 0.0 ) qval[i] = 0.0; // Close the splitting impossible channels
1445 qv += qval[i];
1446 }
1447 // Select the channel
1448 G4double qv1 = 0.0;
1449 G4double ran = qv*G4UniformRand();
1450 G4int index = 0;
1451 for( index=0; index<nCh; ++index ) if( qval[index] > 0.0 )
1452 {
1453 qv1 += qval[index];
1454 if( ran <= qv1 ) break;
1455 }
1456#ifdef debug
1457 G4cout<<"G4QLowEn::PSDI: index="<<index<<" < "<<nCh<<G4endl;
1458#endif
1459 if(index == nCh)
1460 {
1461 G4cout<<"***G4QLowEnergy::PoStDI:Decay is impossible,totM="<<totM<<",GSM="<<tM<<G4endl;
1462 G4cout<<"G4QLowEn::PoStDI:Reaction "<<projPDG<<"+"<<targPDG<<", P="<<Momentum<<G4endl;
1463 G4int nRes=result.size();
1464 if(nRes)G4cout<<"G4QLowEnergy::PoStDI: result exists with N="<<nRes<<" tracks"<<G4endl;
1465// else throw G4QException("***G4QLowEnergy::PostStepDoIt: Can't decay ResidualCompound");
1466 else {
1467 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
1468 FatalException, "Can't decay ResidualCompound");
1469 }
1470 G4double minEx=1000000.; // Prototype of minimal excess
1471 G4bool found = false; // The solution is not found yet
1472 G4int cInd = 0; // Prototype of the best index
1473 G4int ctN = 0; // To
1474 G4int ctZ = 0; // avoid
1475 G4LorentzVector c4M=tot4M; // recalculation
1476 G4double ctM=0.; // of found
1477 for(G4int it=0; it<nRes; ++it)
1478 {
1479 G4Track* iTrack = result[it];
1480 const G4DynamicParticle* iHadron = iTrack->GetDynamicParticle();
1481 G4ParticleDefinition* iParticle=iHadron->GetDefinition();
1482 G4int iPDG = iParticle->GetPDGEncoding();
1483 G4LorentzVector ih4M = projHadron->Get4Momentum();
1484 G4cout<<" G4QLowEn::PoStDI: H["<<it<<"] = [ "<<iPDG<<", "<<ih4M<<" ]"<<G4endl;
1485 G4int iB = iParticle->GetBaryonNumber(); // A of the secondary
1486 G4int iC = G4int(iParticle->GetPDGCharge()); // Z of the secondary
1487 G4LorentzVector new4M = tot4M + ih4M; // Make temporary tot4M
1488 G4int ntN=totN + iB-iC; // Make temporary totN
1489 G4int ntZ=totZ + iC; // Make temporary totZ
1490 G4double ntM = targQPDG.GetNuclMass(ntZ, ntN, 0); // Make temporary GSM
1491 G4double ttM = new4M.m();
1492 if(ttM > ntM) // >>> This is at least one possible solution ! <<<
1493 {
1494 G4double nEx = ttM - ntM;
1495 if(found && nEx < minEx) // This one is better than the previous found
1496 {
1497 cInd = it;
1498 minEx= nEx;
1499 ctN = ntN;
1500 ctZ = ntZ;
1501 c4M = new4M;
1502 ctM = ttM;
1503 mass[0] = ntM;
1504 }
1505 found = true;
1506 }
1507 }
1508 if(found)
1509 {
1510 tot4M = c4M;
1511 totM = ctM;
1512 totN = ctN;
1513 totZ = ctZ;
1514 delete result[cInd];
1515 G4int nR1 = nRes -1;
1516 if(cInd < nR1) result[cInd] = result[nR1];
1517 result.pop_back();
1518 //nRes=nR1; // @@ can be used for the two pt correction
1519 index = 0; // @@ emergency gamma+gamma decay is selected
1520 }
1521 else
1522 {
1523 G4cout<<"***G4QLowEnergy::PoStDI: One hadron coddection did not succeed***"<<G4endl;
1524 if(nRes>1) G4cout<<"***G4QLowEnergy::PoStDI:nRes's big enough to use 2PtCor"<<G4endl;
1525 // throw G4QException("G4QLowEnergy::PostStepDoIt: Can't correct by one particle.");
1526 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
1527 FatalException, "Can't correct by one particle");
1528 }
1529 }
1530 // @@ Convert it to G4QHadrons
1531 G4DynamicParticle* ResSec = new G4DynamicParticle;
1532 G4DynamicParticle* FstSec = new G4DynamicParticle;
1533 G4DynamicParticle* SecSec = new G4DynamicParticle;
1534#ifdef debug
1535 G4cout<<"G4QLowEn::PSDI: Dynamic particles are created pL="<<pL<<G4endl;
1536#endif
1537
1538 G4LorentzVector res4Mom(zeroMom,mass[index]*MeV); // The recoil nucleus prototype
1539 G4double mF=0.;
1540 G4double mS=0.;
1541 G4int rA=totZ+totN;
1542 G4int rZ=totZ;
1543 G4int rL=pL;
1544 G4int complete=3;
1545 G4ParticleDefinition* theDefinition; // Prototype for qfNucleon
1546 switch( index )
1547 {
1548 case 0:
1549 if(!evaporate || rA<2)
1550 {
1551 if(!rZ) theDefinition=aNeutron;
1552 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1553 if(!theDefinition)
1554 {
1555 // G4cerr<<"-Warning-G4LE::PSDI: notDef(0), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1556 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1557 G4ExceptionDescription ed;
1558 ed << "Particle definition is a null pointer: notDef(0), Z= " << rZ
1559 << ", A=" << rA << ", L=" << rL << G4endl;
1560 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
1561 JustWarning, ed);
1562 }
1563 ResSec->SetDefinition( theDefinition );
1564 FstSec->SetDefinition( aGamma );
1565 SecSec->SetDefinition( aGamma );
1566 }
1567 else
1568 {
1569 delete ResSec;
1570 delete FstSec;
1571 delete SecSec;
1572 complete=0;
1573 }
1574 break;
1575 case 1:
1576 rA-=1; // gamma+n
1577 if(!evaporate || rA<2)
1578 {
1579 if(!rZ) theDefinition=aNeutron;
1580 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1581 if(!theDefinition)
1582 {
1583 // G4cerr<<"-Warning-G4LE::PSDI: notDef(1), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1584 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1585 G4ExceptionDescription ed;
1586 ed << "Particle definition is a null pointer: notDef(1), Z= " << rZ
1587 << ", A=" << rA << ", L=" << rL << G4endl;
1588 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0001",
1589 JustWarning, ed);
1590 }
1591 ResSec->SetDefinition( theDefinition );
1592 SecSec->SetDefinition( aGamma );
1593 }
1594 else
1595 {
1596 delete ResSec;
1597 delete SecSec;
1598 complete=1;
1599 }
1600 FstSec->SetDefinition( aNeutron );
1601 mF=mNeut; // First hadron 4-momentum
1602 break;
1603 case 2:
1604 rA-=1;
1605 rZ-=1; // gamma+p
1606 if(!evaporate || rA<2)
1607 {
1608 if(!rZ) theDefinition=aNeutron;
1609 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1610 if(!theDefinition)
1611 {
1612 // G4cerr<<"-Warning-G4LE::PSDI: notDef(2), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1613 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1614 G4ExceptionDescription ed;
1615 ed << "Particle definition is a null pointer: notDef(2), Z= " << rZ
1616 << ", A=" << rA << ", L="<< rL << G4endl;
1617 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0002",
1618 JustWarning, ed);
1619 }
1620 ResSec->SetDefinition( theDefinition );
1621 SecSec->SetDefinition( aGamma );
1622 }
1623 else
1624 {
1625 delete ResSec;
1626 delete SecSec;
1627 complete=1;
1628 }
1629 FstSec->SetDefinition( aProton );
1630 mF=mProt; // First hadron 4-momentum
1631 break;
1632 case 3:
1633 rA-=2;
1634 rZ-=1; // gamma+d
1635 if(!evaporate || rA<2)
1636 {
1637 if(!rZ) theDefinition=aNeutron;
1638 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1639 if(!theDefinition)
1640 {
1641 // G4cerr<<"-Warning-G4LE::PSDI: notDef(3), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1642 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1643 G4ExceptionDescription ed;
1644 ed << "Particle definition is a null pointer: notDef(3), Z= " << rZ
1645 << ", A=" << rA << ", L=" << rL << G4endl;
1646 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0003",
1647 JustWarning, ed);
1648 }
1649 ResSec->SetDefinition( theDefinition );
1650 SecSec->SetDefinition( aGamma );
1651 }
1652 else
1653 {
1654 delete ResSec;
1655 delete SecSec;
1656 complete=1;
1657 }
1658 FstSec->SetDefinition( aDeuteron );
1659 mF=mDeut; // First hadron 4-momentum
1660 break;
1661 case 4:
1662 rA-=3; // gamma+t
1663 rZ-=1;
1664 if(!evaporate || rA<2)
1665 {
1666 if(!rZ) theDefinition=aNeutron;
1667 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1668 if(!theDefinition)
1669 {
1670 // G4cerr<<"-Warning-G4LE::PSDI: notDef(4), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1671 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1672 G4ExceptionDescription ed;
1673 ed << "Particle definition is a null pointer: notDef(4), Z= " << rZ
1674 << ", A=" << rA << ", L=" << rL << G4endl;
1675 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0004",
1676 JustWarning, ed);
1677 }
1678 ResSec->SetDefinition( theDefinition );
1679 SecSec->SetDefinition( aGamma );
1680 }
1681 else
1682 {
1683 delete ResSec;
1684 delete SecSec;
1685 complete=1;
1686 }
1687 FstSec->SetDefinition( aTriton );
1688 mF=mTrit; // First hadron 4-momentum
1689 break;
1690 case 5: // gamma+He3
1691 rA-=3;
1692 rZ-=2;
1693 if(!evaporate || rA<2)
1694 {
1695 if(!rZ) theDefinition=aNeutron;
1696 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1697 if(!theDefinition)
1698 {
1699 // G4cerr<<"-Warning-G4LE::PSDI: notDef(5), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1700 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1701 G4ExceptionDescription ed;
1702 ed << "Particle definition is a null pointer: notDef(5), Z= " << rZ
1703 << ", A=" << rA << ", L=" << rL << G4endl;
1704 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0005",
1705 JustWarning, ed);
1706 }
1707 ResSec->SetDefinition( theDefinition );
1708 SecSec->SetDefinition( aGamma );
1709 }
1710 else
1711 {
1712 delete ResSec;
1713 delete SecSec;
1714 complete=1;
1715 }
1716 FstSec->SetDefinition( aHe3);
1717 mF=mHel3; // First hadron 4-momentum
1718 break;
1719 case 6:
1720 rA-=4;
1721 rZ-=2; // gamma+He4
1722 if(!evaporate || rA<2)
1723 {
1724 if(!rZ) theDefinition=aNeutron;
1725 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1726 if(!theDefinition)
1727 {
1728 // G4cerr<<"-Warning-G4LE::PSDI: notDef(6), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1729 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1730 G4ExceptionDescription ed;
1731 ed << "Particle definition is a null pointer: notDef(6), Z= " << rZ
1732 << ", A=" << rA << ", L=" << rL << G4endl;
1733 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0006",
1734 JustWarning, ed);
1735 }
1736 ResSec->SetDefinition( theDefinition );
1737 SecSec->SetDefinition( aGamma );
1738 }
1739 else
1740 {
1741 delete ResSec;
1742 delete SecSec;
1743 complete=1;
1744 }
1745 FstSec->SetDefinition( anAlpha );
1746 mF=mAlph; // First hadron 4-momentum
1747 break;
1748 case 7:
1749 rA-=2; // n+n
1750 if(rA==1 && !rZ) theDefinition=aNeutron;
1751 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1752 if(!theDefinition)
1753 {
1754 // G4cerr<<"-Warning-G4LE::PSDI: notDef(7), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1755 // throw G4QException("G4QLowE::PostStepDoIt: ParticleDefinition is a null pointer");
1756 G4ExceptionDescription ed;
1757 ed << "Particle definition is a null pointer: notDef(7), Z= " << rZ
1758 << ", A=" << rA << ", L=" << rL << G4endl;
1759 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0007",
1760 JustWarning, ed);
1761 }
1762 ResSec->SetDefinition( theDefinition );
1763 FstSec->SetDefinition( aNeutron );
1764 SecSec->SetDefinition( aNeutron );
1765 mF=mNeut; // First hadron 4-momentum
1766 mS=mNeut; // Second hadron 4-momentum
1767 break;
1768 case 8:
1769 rZ-=1;
1770 rA-=2; // n+p
1771 if(rA==1 && !rZ) theDefinition=aNeutron;
1772 else if(rA==2 && !rZ)
1773 {
1774 index=7;
1775 ResSec->SetDefinition( aDeuteron);
1776 FstSec->SetDefinition( aNeutron );
1777 SecSec->SetDefinition( aNeutron );
1778 mF=mNeut; // First hadron 4-momentum
1779 mS=mNeut; // Second hadron 4-momentum
1780 break;
1781 }
1782 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1783 if(!theDefinition)
1784 {
1785 // G4cerr<<"-Warning-G4LE::PSDI: notDef(8), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1786 // throw G4QException("G4QLowEn::PostStepDoIt: Darticle Definition is a null pointer");
1787 G4ExceptionDescription ed;
1788 ed << "Particle definition is a null pointer: notDef(8), Z= " << rZ
1789 << ", A=" << rA << ", L=" << rL << G4endl;
1790 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0008",
1791 JustWarning, ed);
1792 }
1793 ResSec->SetDefinition( theDefinition );
1794 FstSec->SetDefinition( aNeutron );
1795 SecSec->SetDefinition( aProton );
1796 mF=mNeut; // First hadron 4-momentum
1797 mS=mProt; // Second hadron 4-momentum
1798 break;
1799 case 9:
1800 rZ-=2;
1801 rA-=2; // p+p
1802 if(rA==1 && !rZ) theDefinition=aNeutron;
1803 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1804 if(!theDefinition)
1805 {
1806 // G4cerr<<"-Warning-G4LE::PSDI: notDef(9), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1807 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1808 G4ExceptionDescription ed;
1809 ed << "Particle definition is a null pointer: notDef(9), Z= " << rZ
1810 << ", A=" << rA << ", L=" << rL << G4endl;
1811 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0009",
1812 JustWarning, ed);
1813 }
1814 ResSec->SetDefinition( theDefinition );
1815 FstSec->SetDefinition( aProton );
1816 SecSec->SetDefinition( aProton );
1817 mF=mProt; // First hadron 4-momentum
1818 mS=mProt; // Second hadron 4-momentum
1819 break;
1820 case 10:
1821 rZ-=2;
1822 rA-=3; // p+d
1823 if(rA==1 && !rZ) theDefinition=aNeutron;
1824 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1825 if(!theDefinition)
1826 {
1827 // G4cerr<<"-Warning-G4LE::PSDI: notDef(10), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1828 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1829 G4ExceptionDescription ed;
1830 ed << "Particle definition is a null pointer: notDef(10), Z= " << rZ
1831 << ", A=" << rA << ", L=" << rL << G4endl;
1832 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0010",
1833 JustWarning, ed);
1834 }
1835 ResSec->SetDefinition( theDefinition );
1836 FstSec->SetDefinition( aProton );
1837 SecSec->SetDefinition( aDeuteron );
1838 mF=mProt; // First hadron 4-momentum
1839 mS=mDeut; // Second hadron 4-momentum
1840 break;
1841 case 11:
1842 rZ-=1;
1843 rA-=3; // n+d
1844 if(rA==1 && !rZ) theDefinition=aNeutron;
1845 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1846 if(!theDefinition)
1847 {
1848 // G4cerr<<"-Warning-G4LE::PSDI: notDef(11), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1849 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1850 G4ExceptionDescription ed;
1851 ed << "Particle definition is a null pointer: notDef(0), Z= " << rZ
1852 << ", A=" << rA << ", L=" << rL << G4endl;
1853 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0011",
1854 JustWarning, ed);
1855 }
1856 ResSec->SetDefinition( theDefinition );
1857 FstSec->SetDefinition( aNeutron );
1858 SecSec->SetDefinition( aDeuteron );
1859 mF=mNeut; // First hadron 4-momentum
1860 mS=mDeut; // Second hadron 4-momentum
1861 break;
1862 case 12:
1863 rZ-=2;
1864 rA-=4; // d+d
1865 if(rA==1 && !rZ) theDefinition=aNeutron;
1866 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1867 if(!theDefinition)
1868 {
1869 // G4cerr<<"-Warning-G4LE::PSDI: notDef(12), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1870 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1871 G4ExceptionDescription ed;
1872 ed << "Particle definition is a null pointer: notDef(12), Z= " << rZ
1873 << ", A=" << rA << ", L=" << rL << G4endl;
1874 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0012",
1875 JustWarning, ed);
1876 }
1877 ResSec->SetDefinition( theDefinition );
1878 FstSec->SetDefinition( aDeuteron );
1879 SecSec->SetDefinition( aDeuteron );
1880 mF=mDeut; // First hadron 4-momentum
1881 mS=mDeut; // Second hadron 4-momentum
1882 break;
1883 case 13:
1884 rZ-=2;
1885 rA-=4; // p+t
1886 if(rA==1 && !rZ) theDefinition=aNeutron;
1887 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1888 if(!theDefinition)
1889 {
1890 // G4cerr<<"-Warning-G4LE::PSDI: notDef(13), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1891 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1892 G4ExceptionDescription ed;
1893 ed << "Particle definition is a null pointer: notDef(13), Z= " << rZ
1894 << ", A=" << rA << ", L=" << rL << G4endl;
1895 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0013",
1896 JustWarning, ed);
1897 }
1898 ResSec->SetDefinition( theDefinition );
1899 FstSec->SetDefinition( aProton );
1900 SecSec->SetDefinition( aTriton );
1901 mF=mProt; // First hadron 4-momentum
1902 mS=mTrit; // Second hadron 4-momentum
1903 break;
1904 case 14:
1905 rZ-=1;
1906 rA-=4; // n+t
1907 if(rA==1 && !rZ) theDefinition=aNeutron;
1908 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1909 if(!theDefinition)
1910 {
1911 // G4cerr<<"-Warning-G4LE::PSDI: notDef(14), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1912 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1913 G4ExceptionDescription ed;
1914 ed << "Particle definition is a null pointer: notDef(14), Z= " << rZ
1915 << ", A=" << rA << ", L=" << rL << G4endl;
1916 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0014",
1917 JustWarning, ed);
1918 }
1919 ResSec->SetDefinition( theDefinition );
1920 FstSec->SetDefinition( aNeutron );
1921 SecSec->SetDefinition( aTriton );
1922 mF=mNeut; // First hadron 4-momentum
1923 mS=mTrit; // Second hadron 4-momentum
1924 break;
1925 case 15:
1926 rZ-=3;
1927 rA-=4; // p+He3
1928 if(rA==1 && !rZ) theDefinition=aNeutron;
1929 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1930 if(!theDefinition)
1931 {
1932 // G4cerr<<"-Warning-G4LE::PSDI: notDef(15), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1933 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1934 G4ExceptionDescription ed;
1935 ed << "Particle definition is a null pointer: notDef(15), Z= " << rZ
1936 << ", A=" << rA << ", L=" << rL << G4endl;
1937 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0015",
1938 JustWarning, ed);
1939 }
1940 ResSec->SetDefinition( theDefinition );
1941 FstSec->SetDefinition( aProton);
1942 SecSec->SetDefinition( aHe3 );
1943 mF=mProt; // First hadron 4-momentum
1944 mS=mHel3; // Second hadron 4-momentum
1945 break;
1946 case 16:
1947 rZ-=2;
1948 rA-=4; // n+He3
1949 if(rA==1 && !rZ) theDefinition=aNeutron;
1950 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1951 if(!theDefinition)
1952 {
1953 // G4cerr<<"-Warning-G4LE::PSDI: notDef(16), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1954 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1955 G4ExceptionDescription ed;
1956 ed << "Particle definition is a null pointer: notDef(16), Z= " << rZ
1957 << ", A=" << rA << ", L=" << rL << G4endl;
1958 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0016",
1959 JustWarning, ed);
1960 }
1961 ResSec->SetDefinition( theDefinition );
1962 FstSec->SetDefinition( aNeutron );
1963 SecSec->SetDefinition( aHe3 );
1964 mF=mNeut; // First hadron 4-momentum
1965 mS=mHel3; // Second hadron 4-momentum
1966 break;
1967 case 17:
1968 rZ-=3;
1969 rA-=5; // p+alph
1970 if(rA==1 && !rZ) theDefinition=aNeutron;
1971 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1972 if(!theDefinition)
1973 {
1974 // G4cerr<<"-Warning-G4LE::PSDI: notDef(17), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1975 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1976 G4ExceptionDescription ed;
1977 ed << "Particle definition is a null pointer: notDef(17), Z= " << rZ
1978 << ", A=" << rA << ", L=" << rL << G4endl;
1979 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0017",
1980 JustWarning, ed);
1981 }
1982 ResSec->SetDefinition( theDefinition );
1983 FstSec->SetDefinition( aProton );
1984 SecSec->SetDefinition( anAlpha );
1985 mF=mProt; // First hadron 4-momentum
1986 mS=mAlph; // Second hadron 4-momentum
1987 break;
1988 case 18:
1989 rZ-=2;
1990 rA-=5; // n+alph
1991 if(rA==1 && !rZ) theDefinition=aNeutron;
1992 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
1993 if(!theDefinition)
1994 {
1995 // G4cerr<<"-Warning-G4LE::PSDI: notDef(18), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
1996 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
1997 G4ExceptionDescription ed;
1998 ed << "Particle definition is a null pointer: notDef(18), Z= " << rZ
1999 << ", A=" << rA << ", L=" << rL << G4endl;
2000 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0018",
2001 JustWarning, ed);
2002 }
2003 ResSec->SetDefinition( theDefinition );
2004 FstSec->SetDefinition( aNeutron );
2005 SecSec->SetDefinition( anAlpha );
2006 mF=mNeut; // First hadron 4-momentum
2007 mS=mAlph; // Second hadron 4-momentum
2008 break;
2009 case 19:
2010 rL-=1; // L+gamma (@@ rA inludes rL?)
2011 rA-=1;
2012 if(rA==1 && !rZ) theDefinition=aNeutron;
2013 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2014 if(!theDefinition)
2015 {
2016 // G4cerr<<"-Warning-G4LE::PSDI: notDef(19), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2017 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
2018 G4ExceptionDescription ed;
2019 ed << "Particle definition is a null pointer: notDef(19), Z= " << rZ
2020 << ", A=" << rA << ", L=" << rL << G4endl;
2021 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0019",
2022 JustWarning, ed);
2023 }
2024 ResSec->SetDefinition( theDefinition );
2025 FstSec->SetDefinition( aLambda );
2026 SecSec->SetDefinition( aGamma );
2027 mF=mLamb; // First hadron 4-momentum
2028 break;
2029 case 20:
2030 rL-=1; // L+p (@@ rA inludes rL?)
2031 rZ-=1;
2032 rA-=2;
2033 if(rA==1 && !rZ) theDefinition=aNeutron;
2034 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2035 if(!theDefinition)
2036 {
2037 // G4cerr<<"-Warning-G4LE::PSDI: notDef(20), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2038 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
2039 G4ExceptionDescription ed;
2040 ed << "Particle definition is a null pointer: notDef(20), Z= " << rZ
2041 << ", A=" << rA << ", L=" << rL << G4endl;
2042 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0020",
2043 JustWarning, ed);
2044 }
2045 ResSec->SetDefinition( theDefinition );
2046 FstSec->SetDefinition( aProton );
2047 SecSec->SetDefinition( aLambda );
2048 mF=mProt; // First hadron 4-momentum
2049 mS=mLamb; // Second hadron 4-momentum
2050 break;
2051 case 21:
2052 rL-=1; // L+n (@@ rA inludes rL?)
2053 rA-=2;
2054 if(rA==1 && !rZ) theDefinition=aNeutron;
2055 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2056 if(!theDefinition)
2057 {
2058 // G4cerr<<"-Warning-G4LE::PSDI: notDef(21), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2059 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Pefinition is a null pointer");
2060 G4ExceptionDescription ed;
2061 ed << "Particle definition is a null pointer: notDef(21), Z= " << rZ
2062 << ", A=" << rA << ", L=" << rL << G4endl;
2063 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0021",
2064 JustWarning, ed);
2065 }
2066 ResSec->SetDefinition( theDefinition );
2067 FstSec->SetDefinition( aNeutron );
2068 SecSec->SetDefinition( aLambda );
2069 mF=mNeut; // First hadron 4-momentum
2070 mS=mLamb; // Second hadron 4-momentum
2071 break;
2072 case 22:
2073 rL-=1; // L+d (@@ rA inludes rL?)
2074 rZ-=1;
2075 rA-=3;
2076 if(rA==1 && !rZ) theDefinition=aNeutron;
2077 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2078 if(!theDefinition)
2079 {
2080 // G4cerr<<"-Warning-G4LE::PSDI: notDef(22), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2081 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
2082 G4ExceptionDescription ed;
2083 ed << "Particle definition is a null pointer: notDef(22), Z= " << rZ
2084 << ", A=" << rA << ", L=" << rL << G4endl;
2085 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0022",
2086 JustWarning, ed);
2087 }
2088 ResSec->SetDefinition( theDefinition );
2089 FstSec->SetDefinition( aDeuteron );
2090 SecSec->SetDefinition( aLambda );
2091 mF=mDeut; // First hadron 4-momentum
2092 mS=mLamb; // Second hadron 4-momentum
2093 break;
2094 case 23:
2095 rL-=1; // L+t (@@ rA inludes rL?)
2096 rZ-=1;
2097 rA-=4;
2098 if(rA==1 && !rZ) theDefinition=aNeutron;
2099 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2100 if(!theDefinition)
2101 {
2102 // G4cerr<<"-Warning-G4LE::PSDI: notDef(23), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2103 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
2104 G4ExceptionDescription ed;
2105 ed << "Particle definition is a null pointer: notDef(23), Z= " << rZ
2106 << ", A=" << rA << ", L=" << rL << G4endl;
2107 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0023",
2108 JustWarning, ed);
2109 }
2110 ResSec->SetDefinition( theDefinition );
2111 FstSec->SetDefinition( aTriton );
2112 SecSec->SetDefinition( aLambda );
2113 mF=mTrit; // First hadron 4-momentum
2114 mS=mLamb; // Second hadron 4-momentum
2115 break;
2116 case 24:
2117 rL-=1; // L+He3 (@@ rA inludes rL?)
2118 rZ-=2;
2119 rA-=4;
2120 if(rA==1 && !rZ) theDefinition=aNeutron;
2121 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2122 if(!theDefinition)
2123 {
2124 // G4cerr<<"-Warning-G4LE::PSDI: notDef(24), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2125 // throw G4QException("G4QLowEn::PostStepDoIt: particle definition is a null pointer");
2126 G4ExceptionDescription ed;
2127 ed << "Particle definition is a null pointer: notDef(24), Z= " << rZ
2128 << ", A=" << rA << ", L=" << rL << G4endl;
2129 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0024",
2130 JustWarning, ed);
2131 }
2132 ResSec->SetDefinition( theDefinition );
2133 FstSec->SetDefinition( aHe3 );
2134 SecSec->SetDefinition( aLambda );
2135 mF=mHel3; // First hadron 4-momentum
2136 mS=mLamb; // Second hadron 4-momentum
2137 break;
2138 case 25:
2139 rL-=1; // L+alph (@@ rA inludes rL?)
2140 rZ-=2;
2141 rA-=5;
2142 if(rA==1 && !rZ) theDefinition=aNeutron;
2143 else theDefinition=G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,rL,rZ);
2144 if(!theDefinition)
2145 {
2146 // G4cerr<<"-Warning-G4LE::PSDI: notDef(25), Z="<<rZ<<", A="<<rA<<", L="<<rL<<G4endl;
2147 // throw G4QException("G4QLowEn::PostStepDoIt: Particle Definition is a null pointer");
2148 G4ExceptionDescription ed;
2149 ed << "Particle definition is a null pointer: notDef(25), Z= " << rZ
2150 << ", A=" << rA << ", L=" << rL << G4endl;
2151 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0025",
2152 JustWarning, ed);
2153 }
2154 ResSec->SetDefinition( theDefinition );
2155 FstSec->SetDefinition( anAlpha );
2156 SecSec->SetDefinition( aLambda );
2157 mF=mAlph; // First hadron 4-momentum
2158 mS=mLamb; // Second hadron 4-momentum
2159 break;
2160 }
2161#ifdef debug
2162 G4cout<<"G4QLowEn::PSDI:F="<<mF<<",S="<<mS<<",com="<<complete<<",ev="<<evaporate<<G4endl;
2163#endif
2164 G4LorentzVector fst4Mom(zeroMom,mF); // Prototype of the first hadron 4-momentum
2165 G4LorentzVector snd4Mom(zeroMom,mS); // Prototype of the second hadron 4-momentum
2166 G4LorentzVector dir4Mom=tot4M; // Prototype of the resN decay direction 4-momentum
2167 dir4Mom.setE(tot4M.e()/2.); // Get half energy and total 3-momentum
2168 // @@ Can be repeated to take into account the Coulomb Barrier
2169 if(!G4QHadron(tot4M).CopDecayIn3(fst4Mom,snd4Mom,res4Mom,dir4Mom,cosp))
2170 {
2171 // G4cerr<<"**G4LowEnergy::PoStDoIt:i="<<index<<",tM="<<totM<<"->M1="<<res4Mom.m()<<"+M2="
2172 // <<fst4Mom.m()<<"+M3="<<snd4Mom.m()<<"=="<<res4Mom.m()+fst4Mom.m()+snd4Mom.m()<<G4endl;
2173 // throw G4QException("G4QLowEnergy::PostStepDoIt: Can't decay the Compound");
2174 G4ExceptionDescription ed;
2175 ed << "Can't decay the Compound: i=" << index << ",tM=" << totM << "->M1="
2176 << res4Mom.m() << "+M2=" << fst4Mom.m() << "+M3=" << snd4Mom.m() << "=="
2177 << res4Mom.m()+fst4Mom.m()+snd4Mom.m() << G4endl;
2178 G4Exception("G4QLowEnergy::PostStepDoIt()", "HAD_CHPS_0000",
2179 FatalException, ed);
2180 }
2181#ifdef debug
2182 G4cout<<"G4QLowEn::PSDI:r4M="<<res4Mom<<",f4M="<<fst4Mom<<",s4M="<<snd4Mom<<G4endl;
2183#endif
2184 G4Track* aNewTrack = 0;
2185 if(complete)
2186 {
2187 FstSec->Set4Momentum(fst4Mom);
2188 aNewTrack = new G4Track(FstSec, localtime, position );
2189 aNewTrack->SetWeight(weight); // weighted
2190 aNewTrack->SetTouchableHandle(trTouchable);
2191 result.push_back( aNewTrack );
2192 EnMomConservation-=fst4Mom;
2193#ifdef debug
2194 G4cout<<"G4QLowEn::PSDI: ***Filled*** 1stH4M="<<fst4Mom
2195 <<", PDG="<<FstSec->GetDefinition()->GetPDGEncoding()<<G4endl;
2196#endif
2197 if(complete>2) // Final solution
2198 {
2199 ResSec->Set4Momentum(res4Mom);
2200 aNewTrack = new G4Track(ResSec, localtime, position );
2201 aNewTrack->SetWeight(weight); // weighted
2202 aNewTrack->SetTouchableHandle(trTouchable);
2203 result.push_back( aNewTrack );
2204 EnMomConservation-=res4Mom;
2205#ifdef debug
2206 G4cout<<"G4QLowEn::PSDI: ***Filled*** rA4M="<<res4Mom<<",rZ="<<rZ<<",rA="<<rA<<",rL="
2207 <<rL<<G4endl;
2208#endif
2209 SecSec->Set4Momentum(snd4Mom);
2210 aNewTrack = new G4Track(SecSec, localtime, position );
2211 aNewTrack->SetWeight(weight); // weighted
2212 aNewTrack->SetTouchableHandle(trTouchable);
2213 result.push_back( aNewTrack );
2214 EnMomConservation-=snd4Mom;
2215#ifdef debug
2216 G4cout<<"G4QLowEn::PSDI: ***Filled*** 2ndH4M="<<snd4Mom
2217 <<", PDG="<<SecSec->GetDefinition()->GetPDGEncoding()<<G4endl;
2218#endif
2219 }
2220 else res4Mom+=snd4Mom;
2221 }
2222 else res4Mom=tot4M;
2223 if(complete<3) // Evaporation of the residual must be done
2224 {
2225 G4QHadron* rHadron = new G4QHadron(90000000+999*rZ+rA,res4Mom); // Input hadron-nucleus
2226 G4QHadronVector* evaHV = new G4QHadronVector; // Output vector of hadrons (delete!)
2227 Nuc.EvaporateNucleus(rHadron, evaHV); // here a pion can appear !
2228 G4int nOut=evaHV->size();
2229 for(i=0; i<nOut; i++)
2230 {
2231 G4QHadron* curH = (*evaHV)[i];
2232 G4int hPDG=curH->GetPDGCode();
2233 G4LorentzVector h4Mom=curH->Get4Momentum();
2234 EnMomConservation-=h4Mom;
2235#ifdef debug
2236 G4cout<<"G4QLowEn::PSDI: ***FillingCand#"<<i<<"*** evaH="<<hPDG<<h4Mom<<G4endl;
2237#endif
2238 if (hPDG==90000001 || hPDG==2112) theDefinition = aNeutron;
2239 else if(hPDG==90001000 || hPDG==2212) theDefinition = aProton;
2240 else if(hPDG==91000000 || hPDG==3122) theDefinition = aLambda;
2241 else if(hPDG== 22 ) theDefinition = aGamma;
2242 else if(hPDG== 111) theDefinition = aPiZero;
2243 else if(hPDG== 211) theDefinition = aPiPlus;
2244 else if(hPDG== -211) theDefinition = aPiMinus;
2245 else
2246 {
2247 G4int hZ=curH->GetCharge();
2248 G4int hA=curH->GetBaryonNumber();
2249 G4int hS=curH->GetStrangeness();
2250 theDefinition = G4ParticleTable::GetParticleTable()->FindIon(hZ,hA,hS,0); // ion
2251 }
2252 if(theDefinition)
2253 {
2254 G4DynamicParticle* theEQH = new G4DynamicParticle(theDefinition,h4Mom);
2255 G4Track* evaQH = new G4Track(theEQH, localtime, position );
2256 evaQH->SetWeight(weight); // weighted
2257 evaQH->SetTouchableHandle(trTouchable);
2258 result.push_back( evaQH );
2259 }
2260 else G4cerr<<"-Warning-G4QLowEnergy::PostStepDoIt: Bad secondary PDG="<<hPDG<<G4endl;
2261 }
2262 }
2263 // algorithm implementation --- STOPS HERE
2264 G4int nres=result.size();
2265 aParticleChange.SetNumberOfSecondaries(nres);
2266 for(i=0; i<nres; ++i) aParticleChange.AddSecondary(result[i]);
2267#ifdef debug
2268 G4cout<<"G4QLowEnergy::PostStepDoIt:*** PostStepDoIt is done ***"<<G4endl;
2269#endif
2270 return G4VDiscreteProcess::PostStepDoIt(track, step);
2271}
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:64
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4QHadronVector
std::vector< G4QHadron * > G4QHadronVector
Definition: G4QHadronVector.hh:56
G4RandomDirection
G4ThreeVector G4RandomDirection()
Definition: G4RandomDirection.hh:55
fStopAndKill
@ fStopAndKill
Definition: G4TrackStatus.hh:56
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:59
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::SetDefinition
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
Definition: G4DynamicParticle.cc:271
G4DynamicParticle::Get4Momentum
G4LorentzVector Get4Momentum() const
G4DynamicParticle::Set4Momentum
void Set4Momentum(const G4LorentzVector &momentum)
Definition: G4DynamicParticle.cc:351
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4Lambda::Lambda
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
G4ParticleChange::AddSecondary
void AddSecondary(G4Track *aSecondary)
Definition: G4ParticleChange.cc:218
G4ParticleChange::ProposeEnergy
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4ParticleChange::Initialize
virtual void Initialize(const G4Track &)
Definition: G4ParticleChange.cc:228
G4ParticleDefinition::GetParticleType
const G4String & GetParticleType() const
Definition: G4ParticleDefinition.hh:136
G4ParticleDefinition::GetQuarkContent
G4int GetQuarkContent(G4int flavor) const
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition: G4ParticleDefinition.hh:137
G4ParticleDefinition::GetAntiQuarkContent
G4int GetAntiQuarkContent(G4int flavor) const
G4ParticleTable::FindIon
G4ParticleDefinition * FindIon(G4int atomicNumber, G4int atomicMass, G4double excitationEnergy)
Definition: G4ParticleTable.cc:399
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:63
G4PionMinus::PionMinus
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
G4PionPlus::PionPlus
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4PionZero::PionZero
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
G4QHadron::Get4Momentum
G4LorentzVector Get4Momentum() const
Definition: G4QHadron.hh:79
G4QHadron::GetBaryonNumber
G4int GetBaryonNumber() const
Definition: G4QHadron.hh:181
G4QHadron::GetCharge
G4int GetCharge() const
Definition: G4QHadron.hh:179
G4QHadron::GetPDGCode
G4int GetPDGCode() const
Definition: G4QHadron.hh:170
G4QHadron::GetStrangeness
G4int GetStrangeness() const
Definition: G4QHadron.hh:180
G4QNucleus::CoulombBarGen
G4double CoulombBarGen(const G4double &rZ, const G4double &rA, const G4double &cZ, const G4double &cA)
Definition: G4QNucleus.cc:3358
G4QNucleus::EvaporateNucleus
void EvaporateNucleus(G4QHadron *hA, G4QHadronVector *oHV)
Definition: G4QNucleus.cc:4171
G4QPDGCode::GetMass
G4double GetMass()
Definition: G4QPDGCode.cc:693
G4QPDGCode::GetNuclMass
G4double GetNuclMass(G4int Z, G4int N, G4int S)
Definition: G4QPDGCode.cc:766
G4QProtonElasticCrossSection::GetPointer
static G4VQCrossSection * GetPointer()
Definition: G4QProtonElasticCrossSection.cc:143
G4Track::SetWeight
void SetWeight(G4double aValue)
G4Track::SetTouchableHandle
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4VDiscreteProcess::PostStepDoIt
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
Definition: G4VDiscreteProcess.cc:112
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4VParticleChange::SetNumberOfSecondaries
void SetNumberOfSecondaries(G4int totSecondaries)
G4VProcess::aParticleChange
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
G4VQCrossSection::GetExchangeT
virtual G4double GetExchangeT(G4int tZ, G4int tN, G4int pPDG)
Definition: G4VQCrossSection.cc:70
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
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76

◆ SwitchOffEvaporation()

void G4QLowEnergy::SwitchOffEvaporation ( )
inline

Definition at line 105 of file G4QLowEnergy.hh.

105{evaporate=false;} // Gammas instead of evaporation

◆ SwitchOnEvaporation()

void G4QLowEnergy::SwitchOnEvaporation ( )
inline

Definition at line 103 of file G4QLowEnergy.hh.

103{evaporate=true;} // Evaporation instead of gammas (default)
The documentation for this class was generated from the following files: