Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4QDiffraction Class Reference

#include <G4QDiffraction.hh>

+ Inheritance diagram for G4QDiffraction:

Public Member Functions

 G4QDiffraction (const G4String &processName="CHIPS_DiffractiveInteraction")
 
 ~G4QDiffraction ()
 
G4bool IsApplicable (const G4ParticleDefinition &particle)
 
G4double GetMeanFreePath (const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
 
G4VParticleChangePostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
 
G4LorentzVector GetEnegryMomentumConservation ()
 
G4int GetNumberOfNeutronsInTarget ()
 
- 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 G4VParticleChangePostStepDoIt (const G4Track &, const G4Step &)
 
virtual G4double AlongStepGetPhysicalInteractionLength (const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
 
virtual G4double AtRestGetPhysicalInteractionLength (const G4Track &, G4ForceCondition *)
 
virtual G4VParticleChangeAtRestDoIt (const G4Track &, const G4Step &)
 
virtual G4VParticleChangeAlongStepDoIt (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 G4VParticleChangePostStepDoIt (const G4Track &track, const G4Step &stepData)=0
 
virtual G4VParticleChangeAlongStepDoIt (const G4Track &track, const G4Step &stepData)=0
 
virtual G4VParticleChangeAtRestDoIt (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 G4StringGetPhysicsTableFileName (const G4ParticleDefinition *, const G4String &directory, const G4String &tableName, G4bool ascii=false)
 
const G4StringGetProcessName () 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 G4ProcessManagerGetProcessManager ()
 
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 G4StringGetProcessTypeName (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 G4ProcessManageraProcessManager
 
G4VParticleChangepParticleChange
 
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 73 of file G4QDiffraction.hh.

Constructor & Destructor Documentation

◆ G4QDiffraction()

G4QDiffraction::G4QDiffraction ( const G4String processName = "CHIPS_DiffractiveInteraction")

Definition at line 64 of file G4QDiffraction.cc.

64 :
65 G4VDiscreteProcess(processName, fHadronic)
66{
67 G4HadronicDeprecate("G4QDiffraction");
68
69#ifdef debug
70 G4cout<<"G4QDiffraction::Constructor is called processName="<<processName<<G4endl;
71#endif
72 if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
73 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part. max)
74}
#define G4HadronicDeprecate(name)
@ fHadronic
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
static G4QCHIPSWorld * Get()
G4QParticleVector * GetParticles(G4int nOfParts=0)
G4int verboseLevel
Definition: G4VProcess.hh:368
const G4String & GetProcessName() const
Definition: G4VProcess.hh:379

◆ ~G4QDiffraction()

G4QDiffraction::~G4QDiffraction ( )

Definition at line 77 of file G4QDiffraction.cc.

77{}

Member Function Documentation

◆ GetEnegryMomentumConservation()

G4LorentzVector G4QDiffraction::GetEnegryMomentumConservation ( )

Definition at line 80 of file G4QDiffraction.cc.

80{return EnMomConservation;}

◆ GetMeanFreePath()

G4double G4QDiffraction::GetMeanFreePath ( const G4Track aTrack,
G4double  previousStepSize,
G4ForceCondition condition 
)
virtual

Implements G4VDiscreteProcess.

Definition at line 87 of file G4QDiffraction.cc.

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

Referenced by PostStepDoIt().

◆ GetNumberOfNeutronsInTarget()

G4int G4QDiffraction::GetNumberOfNeutronsInTarget ( )

Definition at line 82 of file G4QDiffraction.cc.

82{return nOfNeutrons;}

◆ IsApplicable()

G4bool G4QDiffraction::IsApplicable ( const G4ParticleDefinition particle)
virtual

Reimplemented from G4VProcess.

Definition at line 215 of file G4QDiffraction.cc.

216{
217 if (particle == *( G4Proton::Proton() )) return true;
218 else if (particle == *( G4Neutron::Neutron() )) return true;
219 //else if (particle == *( G4MuonMinus::MuonMinus() )) return true;
220 //else if (particle == *( G4TauPlus::TauPlus() )) return true;
221 //else if (particle == *( G4TauMinus::TauMinus() )) return true;
222 //else if (particle == *( G4Electron::Electron() )) return true;
223 //else if (particle == *( G4Positron::Positron() )) return true;
224 //else if (particle == *( G4Gamma::Gamma() )) return true;
225 //else if (particle == *( G4MuonPlus::MuonPlus() )) return true;
226 //else if (particle == *(G4AntiNeutrinoMu::AntiNeutrinoMu())) return true;
227 //else if (particle == *( G4NeutrinoMu::NeutrinoMu() )) return true;
228 //else if (particle == *( G4PionMinus::PionMinus() )) return true;
229 //else if (particle == *( G4PionPlus::PionPlus() )) return true;
230 //else if (particle == *( G4KaonPlus::KaonPlus() )) return true;
231 //else if (particle == *( G4KaonMinus::KaonMinus() )) return true;
232 //else if (particle == *( G4KaonZeroLong::KaonZeroLong() )) return true;
233 //else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true;
234 //else if (particle == *( G4Lambda::Lambda() )) return true;
235 //else if (particle == *( G4SigmaPlus::SigmaPlus() )) return true;
236 //else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true;
237 //else if (particle == *( G4SigmaZero::SigmaZero() )) return true;
238 //else if (particle == *( G4XiMinus::XiMinus() )) return true;
239 //else if (particle == *( G4XiZero::XiZero() )) return true;
240 //else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true;
241 //else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true;
242 //else if (particle == *( G4AntiProton::AntiProton() )) return true;
243#ifdef debug
244 G4cout<<"***>>G4QDiffraction::IsApplicable: projPDG="<<particle.GetPDGEncoding()<<G4endl;
245#endif
246 return false;
247}

Referenced by GetMeanFreePath(), and PostStepDoIt().

◆ PostStepDoIt()

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

Reimplemented from G4VDiscreteProcess.

Definition at line 249 of file G4QDiffraction.cc.

250{
251 static const G4double mProt= G4QPDGCode(2212).GetMass(); // CHIPS proton Mass in MeV
252 static const G4double mNeut= G4QPDGCode(2112).GetMass(); // CHIPS neutron Mass in MeV
253 static const G4double mPion= G4QPDGCode(111).GetMass(); // CHIPS Pi0 Mass in MeV
254 static G4QDiffractionRatio* diffRatio;
255 //
256 //-------------------------------------------------------------------------------------
257 static G4bool CWinit = true; // CHIPS Warld needs to be initted
258 if(CWinit)
259 {
260 CWinit=false;
261 G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
263 }
264 //-------------------------------------------------------------------------------------
265 const G4DynamicParticle* projHadron = track.GetDynamicParticle();
266 const G4ParticleDefinition* particle=projHadron->GetDefinition();
267#ifdef debug
268 G4cout<<"G4QDiffraction::PostStepDoIt: Before the GetMeanFreePath is called In4M="
269 <<projHadron->Get4Momentum()<<" of PDG="<<particle->GetPDGEncoding()<<", Type="
270 <<particle->GetParticleType()<<",SubType="<<particle->GetParticleSubType()<<G4endl;
271#endif
273 GetMeanFreePath(track, -27., &cond); // @@ ?? jus to update parameters?
274#ifdef debug
275 G4cout<<"G4QDiffraction::PostStepDoIt: After GetMeanFreePath is called"<<G4endl;
276#endif
277 G4LorentzVector proj4M=(projHadron->Get4Momentum())/MeV; // Convert to MeV!
278 G4double momentum = projHadron->GetTotalMomentum()/MeV; // 3-momentum of the Proj in MeV
279 G4double Momentum = proj4M.rho(); // @@ Just for the test purposes
280 if(std::fabs(Momentum-momentum)>.000001)
281 G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt:P_IU="<<Momentum<<"#"<<momentum<<G4endl;
282#ifdef pdebug
283 G4cout<<"G4QDiffraction::PostStepDoIt: pP(IU)="<<Momentum<<"="<<momentum
284 <<",proj4M="<<proj4M<<", projM="<<proj4M.m()<<G4endl;
285#endif
286 if (!IsApplicable(*particle)) // Check applicability
287 {
288 G4cerr<<"G4QDiffraction::PostStepDoIt: Only NA is implemented."<<G4endl;
289 return 0;
290 }
291 const G4Material* material = track.GetMaterial(); // Get the current material
292 G4int Z=0;
293 const G4ElementVector* theElementVector = material->GetElementVector();
294 G4int nE=material->GetNumberOfElements();
295#ifdef debug
296 G4cout<<"G4QDiffraction::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
297#endif
298 G4int projPDG=0; // PDG Code prototype for the captured hadron
299 // Not all these particles are implemented yet (see Is Applicable)
300 if (particle == G4Proton::Proton() ) projPDG= 2212;
301 else if (particle == G4Neutron::Neutron() ) projPDG= 2112;
302 //else if (particle == G4PionMinus::PionMinus() ) projPDG= -211;
303 //else if (particle == G4PionPlus::PionPlus() ) projPDG= 211;
304 //else if (particle == G4KaonPlus::KaonPlus() ) projPDG= 321;
305 //else if (particle == G4KaonMinus::KaonMinus() ) projPDG= -321;
306 //else if (particle == G4KaonZeroLong::KaonZeroLong() ) projPDG= 130;
307 //else if (particle == G4KaonZeroShort::KaonZeroShort() ) projPDG= 310;
308 //else if (particle == G4MuonPlus::MuonPlus() ) projPDG= -13;
309 //else if (particle == G4MuonMinus::MuonMinus() ) projPDG= 13;
310 //else if (particle == G4NeutrinoMu::NeutrinoMu() ) projPDG= 14;
311 //else if (particle == G4AntiNeutrinoMu::AntiNeutrinoMu() ) projPDG= -14;
312 //else if (particle == G4Electron::Electron() ) projPDG= 11;
313 //else if (particle == G4Positron::Positron() ) projPDG= -11;
314 //else if (particle == G4NeutrinoE::NeutrinoE() ) projPDG= 12;
315 //else if (particle == G4AntiNeutrinoE::AntiNeutrinoE() ) projPDG= -12;
316 //else if (particle == G4Gamma::Gamma() ) projPDG= 22;
317 //else if (particle == G4TauPlus::TauPlus() ) projPDG= -15;
318 //else if (particle == G4TauMinus::TauMinus() ) projPDG= 15;
319 //else if (particle == G4NeutrinoTau::NeutrinoTau() ) projPDG= 16;
320 //else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) projPDG= -16;
321 //else if (particle == G4Lambda::Lambda() ) projPDG= 3122;
322 //else if (particle == G4SigmaPlus::SigmaPlus() ) projPDG= 3222;
323 //else if (particle == G4SigmaMinus::SigmaMinus() ) projPDG= 3112;
324 //else if (particle == G4SigmaZero::SigmaZero() ) projPDG= 3212;
325 //else if (particle == G4XiMinus::XiMinus() ) projPDG= 3312;
326 //else if (particle == G4XiZero::XiZero() ) projPDG= 3322;
327 //else if (particle == G4OmegaMinus::OmegaMinus() ) projPDG= 3334;
328 //else if (particle == G4AntiNeutron::AntiNeutron() ) projPDG=-2112;
329 //else if (particle == G4AntiProton::AntiProton() ) projPDG=-2212;
330#ifdef debug
331 G4int prPDG=particle->GetPDGEncoding();
332 G4cout<<"G4QDiffraction::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
333#endif
334 if(!projPDG)
335 {
336 G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt:UndefProjHadron(PDG=0) ->ret 0"<<G4endl;
337 return 0;
338 }
339 //G4double pM2=proj4M.m2(); // in MeV^2
340 //G4double pM=std::sqrt(pM2); // in MeV
341 G4double pM=mNeut;
342 if(projPDG==2112) pM=mProt;
343 // Element treatment
344 G4int EPIM=ElProbInMat.size();
345#ifdef debug
346 G4cout<<"G4QDiffra::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
347#endif
348 G4int i=0;
349 if(EPIM>1)
350 {
351 G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
352 for(i=0; i<nE; ++i)
353 {
354#ifdef debug
355 G4cout<<"G4QDiffra::PostStepDoIt: EPM["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
356#endif
357 if (rnd<ElProbInMat[i]) break;
358 }
359 if(i>=nE) i=nE-1; // Top limit for the Element
360 }
361 G4Element* pElement=(*theElementVector)[i];
362 Z=static_cast<G4int>(pElement->GetZ());
363#ifdef debug
364 G4cout<<"G4QDiffraction::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
365#endif
366 if(Z<=0)
367 {
368 G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt: Element with Z="<<Z<<G4endl;
369 if(Z<0) return 0;
370 }
371 std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
372 std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i]
373 G4int nofIsot=SPI->size(); // #of isotopes in the element i
374#ifdef debug
375 G4cout<<"G4QDiffraction::PostStepDoIt: nI="<<nofIsot<<", T="<<(*SPI)[nofIsot-1]<<G4endl;
376#endif
377 G4int j=0;
378 if(nofIsot>1)
379 {
380 G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element
381 for(j=0; j<nofIsot; ++j)
382 {
383#ifdef debug
384 G4cout<<"G4QDiffraction::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<",r="<<rndI<<G4endl;
385#endif
386 if(rndI < (*SPI)[j]) break;
387 }
388 if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope
389 }
390 G4int N =(*IsN)[j]; ; // Randomized number of neutrons
391#ifdef debug
392 G4cout<<"G4QDiffraction::PostStepDoIt: j="<<i<<", N(isotope)="<<N<<", MeV="<<MeV<<G4endl;
393#endif
394 if(N<0)
395 {
396 G4cerr<<"-Warning-G4QDiffraction::PostStepDoIt: Isotope Z="<<Z<<" has 0>N="<<N<<G4endl;
397 return 0;
398 }
399 nOfNeutrons=N; // Remember it for the energy-momentum check
400#ifdef debug
401 G4cout<<"G4QDiffraction::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
402#endif
403 if(N<0)
404 {
405 G4cerr<<"*Warning*G4QDiffraction::PostStepDoIt:Element with N="<<N<< G4endl;
406 return 0;
407 }
409#ifdef debug
410 G4cout<<"G4QDiffraction::PostStepDoIt: track is initialized"<<G4endl;
411#endif
412 G4double weight = track.GetWeight();
413 G4double localtime = track.GetGlobalTime();
414 G4ThreeVector position = track.GetPosition();
415#ifdef debug
416 G4cout<<"G4QDiffraction::PostStepDoIt: before Touchable extraction"<<G4endl;
417#endif
418 G4TouchableHandle trTouchable = track.GetTouchableHandle();
419#ifdef debug
420 G4cout<<"G4QDiffraction::PostStepDoIt: Touchable is extracted"<<G4endl;
421#endif
422 G4int targPDG=90000000+Z*1000+N; // CHIPS PDG Code of the target nucleus
423 G4QPDGCode targQPDG(targPDG); // @@ use G4Ion and get rid of CHIPS World
424 G4double tM=targQPDG.GetMass(); // CHIPS final nucleus mass in MeV
425 G4double kinEnergy= projHadron->GetKineticEnergy()*MeV; // Kin energy in MeV (Is *MeV n?)
426 G4ParticleMomentum dir = projHadron->GetMomentumDirection();// It is a unit three-vector
427 G4LorentzVector tot4M=proj4M+G4LorentzVector(0.,0.,0.,tM); // Total 4-mom of the reaction
428#ifdef debug
429 G4cout<<"G4QDiffraction::PostStepDoIt: tM="<<tM<<",p4M="<<proj4M<<",t4M="<<tot4M<<G4endl;
430#endif
431 EnMomConservation=tot4M; // Total 4-mom of reaction for E/M conservation
432 // @@ Probably this is not necessary any more
433#ifdef debug
434 G4cout<<"G4QDiff::PSDI:false,P="<<Momentum<<",Z="<<Z<<",N="<<N<<",PDG="<<projPDG<<G4endl;
435#endif
436 G4double xSec=CalculateXS(Momentum, Z, N, projPDG); // Recalculate CrossSection
437#ifdef debug
438 G4cout<<"G4QDiffra::PSDI:PDG="<<projPDG<<",P="<<Momentum<<",XS="<<xSec/millibarn<<G4endl;
439#endif
440#ifdef nandebug
441 if(xSec>0. || xSec<0. || xSec==0);
442 else G4cout<<"-Warning-G4QDiffraction::PostSDI: *NAN* xSec="<<xSec/millibarn<<G4endl;
443#endif
444 // @@ check a possibility to separate p, n, or alpha (!)
445 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing
446 {
447#ifdef pdebug
448 G4cerr<<"*Warning*G4QDiffraction::PSDoIt:*Zero cross-section* PDG="<<projPDG
449 <<",tPDG="<<targPDG<<",P="<<Momentum<<G4endl;
450#endif
451 //Do Nothing Action insead of the reaction
455 return G4VDiscreteProcess::PostStepDoIt(track,step);
456 }
457 G4double totCMMass=tot4M.m(); // Total CM mass, pM=projectileMass, tM=targetMass
458 if(totCMMass < mPion+pM+tM) // The diffraction reaction is impossible -> Do Nothing
459 {
460#ifdef pdebug
461 G4cerr<<"*Warning*G4QDiffraction::PSDoIt:*Below Diffraction Threshold* cmM="<<totCMMass
462 <<">pM="<<pM<<"+tM="<<tM<<"+pi0="<<mPion<<"=="<<pM+tM+mPion<<G4endl;
463#endif
464 //Do Nothing Action insead of the reaction
468 return G4VDiscreteProcess::PostStepDoIt(track,step);
469 }
470 // Kill interacting hadron
472 G4QHadronVector* out=diffRatio->TargFragment(projPDG, proj4M, Z, N);
473 G4int nSec=out->size(); // #of secondaries in the diffraction reaction
474 G4DynamicParticle* theSec=0; // A prototype for secondary for the secondary
475 G4LorentzVector dif4M(0.,0.,0.,0.); // Prototype for the secondary 4-momentum
476 G4int difPDG=0; // PDG code of the secondary
477 G4QHadron* difQH=0; // Prototype for a Q-secondary
478#ifdef pdebug
479 G4cout<<"G4QDiffraction::PostStepDoIt: =---=found=---= nSecondaries="<<nSec<<G4endl;
480#endif
481 for(i=0; i<nSec; i++)
482 {
483 difQH = (*out)[i];
484 difPDG= difQH->GetPDGCode();
485 G4ParticleDefinition* theDefinition=0;
486 if (difPDG==2212 || difPDG==90001000) theDefinition=G4Proton::Proton();
487 else if(difPDG==2112 || difPDG==90000001) theDefinition=G4Neutron::Neutron();
488 else if(difPDG== 22) theDefinition=G4Gamma::Gamma();
489 else if(difPDG== 111) theDefinition=G4PionZero::PionZero();
490 else if(difPDG==-211 || difPDG==89999001) theDefinition=G4PionMinus::PionMinus();
491 else if(difPDG== 211 || difPDG==90000999) theDefinition=G4PionPlus::PionPlus();
492 else if(difPDG== 321 || difPDG==89001000) theDefinition=G4KaonPlus::KaonPlus();
493 else if(difPDG==-321 || difPDG==90999000) theDefinition=G4KaonMinus::KaonMinus();
494 else if(difPDG== 130 || difPDG==-311 || difPDG==89000001)
495 theDefinition=G4KaonZeroLong::KaonZeroLong();
496 else if(difPDG== 310 || difPDG== 311 || difPDG==90999999)
497 theDefinition=G4KaonZeroShort::KaonZeroShort();
498 else if(difPDG==3122 || difPDG==91000000) theDefinition=G4Lambda::Lambda();
499 else if(difPDG== 3222) theDefinition=G4SigmaPlus::SigmaPlus();
500 else if(difPDG== 3112) theDefinition=G4SigmaMinus::SigmaMinus();
501 else if(difPDG== 3212) theDefinition=G4SigmaZero::SigmaZero();
502 else if(difPDG== 3312) theDefinition=G4XiMinus::XiMinus();
503 else if(difPDG== 3322) theDefinition=G4XiZero::XiZero();
504 else if(difPDG== 3334) theDefinition=G4OmegaMinus::OmegaMinus();
505 else if(difPDG==-2112) theDefinition=G4AntiNeutron::AntiNeutron();
506 else if(difPDG==-2212) theDefinition=G4AntiProton::AntiProton();
507 else if(difPDG==-3122) theDefinition=G4AntiLambda::AntiLambda();
508 else if(difPDG==-3222) theDefinition=G4AntiSigmaPlus::AntiSigmaPlus();
509 else if(difPDG==-3112) theDefinition=G4AntiSigmaMinus::AntiSigmaMinus();
510 else if(difPDG==-3212) theDefinition=G4AntiSigmaZero::AntiSigmaZero();
511 else if(difPDG==-3312) theDefinition=G4AntiXiMinus::AntiXiMinus();
512 else if(difPDG==-3322) theDefinition=G4AntiXiZero::AntiXiZero();
513 else if(difPDG==-3334) theDefinition=G4AntiOmegaMinus::AntiOmegaMinus();
514 else if(difPDG== -11) theDefinition=G4Electron::Electron();
515 else if(difPDG== -13) theDefinition=G4MuonMinus::MuonMinus();
516 else if(difPDG== 11) theDefinition=G4Positron::Positron();
517 else if(difPDG== 13) theDefinition=G4MuonPlus::MuonPlus();
518 else
519 {
520 Z = difQH->GetCharge();
521 G4int B = difQH->GetBaryonNumber();
522 G4int S = difQH->GetStrangeness();
523 if(S||Z>B||Z<0)G4cout<<"-Warning-G4QDif::PoStDoIt:Z="<<Z<<",A="<<B<<",S="<<S<<G4endl;
524 theDefinition = G4ParticleTable::GetParticleTable()->FindIon(Z,B,0,0);
525#ifdef pdebug
526 G4cout<<"G4QDiffraction::PoStDoIt:Ion,Z="<<Z<<",A="<<B<<",D="<<theDefinition<<G4endl;
527#endif
528 }
529 if(theDefinition)
530 {
531 theSec = new G4DynamicParticle; // A secondary for the recoil hadron
532 theSec->SetDefinition(theDefinition);
533 dif4M = difQH->Get4Momentum();
534 EnMomConservation-=dif4M;
535 theSec->Set4Momentum(dif4M);
536 G4Track* aNewTrack = new G4Track(theSec, localtime, position );
537 aNewTrack->SetWeight(weight); // weighted
538 aNewTrack->SetTouchableHandle(trTouchable);
539 aParticleChange.AddSecondary( aNewTrack );
540#ifdef pdebug
541 G4cout<<"G4QDiffraction::PostStepDoIt: Filled 4M="<<dif4M<<", PDG="<<difPDG<<G4endl;
542#endif
543 }
544 else G4cout<<"-Warning-G4QDif::PSDI: Lost PDG="<<difPDG<<", Z="<<difQH->GetCharge()
545 <<", A="<<difQH->GetBaryonNumber()<<",S ="<<difQH->GetStrangeness()<<G4endl;
546 delete difQH; // Clean up the output QHadrons
547 }
548 delete out; // Delete the output QHadron-vector
549#ifdef debug
550 G4cout<<"G4QDiffraction::PostStepDoIt:*** PostStepDoIt is done ***"<<G4endl;
551#endif
552 return G4VDiscreteProcess::PostStepDoIt(track, step);
553}
G4ForceCondition
CLHEP::HepLorentzVector G4LorentzVector
std::vector< G4QHadron * > G4QHadronVector
@ fStopAndKill
bool G4bool
Definition: G4Types.hh:67
#define G4UniformRand()
Definition: Randomize.hh:53
static G4AntiLambda * AntiLambda()
static G4AntiNeutron * AntiNeutron()
static G4AntiOmegaMinus * AntiOmegaMinus()
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
static G4AntiSigmaMinus * AntiSigmaMinus()
static G4AntiSigmaPlus * AntiSigmaPlus()
static G4AntiSigmaZero * AntiSigmaZero()
static G4AntiXiMinus * AntiXiMinus()
static G4AntiXiZero * AntiXiZero()
const G4ThreeVector & GetMomentumDirection() const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
static G4Electron * Electron()
Definition: G4Electron.cc:94
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
static G4KaonZeroLong * KaonZeroLong()
static G4KaonZeroShort * KaonZeroShort()
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
static G4MuonMinus * MuonMinus()
Definition: G4MuonMinus.cc:100
static G4MuonPlus * MuonPlus()
Definition: G4MuonPlus.cc:99
static G4OmegaMinus * OmegaMinus()
void AddSecondary(G4Track *aSecondary)
void ProposeEnergy(G4double finalEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
virtual void Initialize(const G4Track &)
const G4String & GetParticleType() const
const G4String & GetParticleSubType() const
G4ParticleDefinition * FindIon(G4int atomicNumber, G4int atomicMass, G4double excitationEnergy)
static G4ParticleTable * GetParticleTable()
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
static G4Positron * Positron()
Definition: G4Positron.cc:94
static G4QDiffractionRatio * GetPointer()
G4QHadronVector * TargFragment(G4int pPDG, G4LorentzVector p4M, G4int tgZ, G4int tgN)
G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
G4LorentzVector Get4Momentum() const
Definition: G4QHadron.hh:79
G4int GetBaryonNumber() const
Definition: G4QHadron.hh:181
G4int GetCharge() const
Definition: G4QHadron.hh:179
G4int GetPDGCode() const
Definition: G4QHadron.hh:170
G4int GetStrangeness() const
Definition: G4QHadron.hh:180
G4double GetMass()
Definition: G4QPDGCode.cc:693
static G4SigmaMinus * SigmaMinus()
static G4SigmaPlus * SigmaPlus()
Definition: G4SigmaPlus.cc:108
static G4SigmaZero * SigmaZero()
Definition: G4SigmaZero.cc:99
void SetWeight(G4double aValue)
void SetTouchableHandle(const G4TouchableHandle &apValue)
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
static G4XiMinus * XiMinus()
Definition: G4XiMinus.cc:106
static G4XiZero * XiZero()
Definition: G4XiZero.cc:106

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