Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4QuasiElasticChannel Class Reference

#include <G4QuasiElasticChannel.hh>

Public Member Functions

 G4QuasiElasticChannel ()
 
 ~G4QuasiElasticChannel ()
 
G4double GetFraction (G4Nucleus &theNucleus, const G4DynamicParticle &thePrimary)
 
G4KineticTrackVectorScatter (G4Nucleus &theNucleus, const G4DynamicParticle &thePrimary)
 

Detailed Description

Definition at line 52 of file G4QuasiElasticChannel.hh.

Constructor & Destructor Documentation

G4QuasiElasticChannel::G4QuasiElasticChannel ( )

Definition at line 62 of file G4QuasiElasticChannel.cc.

63  : theQuasiElastic(new G4QuasiElRatios),
64  the3DNucleus(new G4Fancy3DNucleus) {
65 }
G4QuasiElasticChannel::~G4QuasiElasticChannel ( )

Definition at line 67 of file G4QuasiElasticChannel.cc.

68 {
69  delete the3DNucleus;
70  delete theQuasiElastic;
71 }

Member Function Documentation

G4double G4QuasiElasticChannel::GetFraction ( G4Nucleus theNucleus,
const G4DynamicParticle thePrimary 
)

Definition at line 73 of file G4QuasiElasticChannel.cc.

75 {
76  #ifdef debug_scatter
77  G4cout << "G4QuasiElasticChannel:: P=" << thePrimary.GetTotalMomentum()
78  << ", pPDG=" << thePrimary.GetDefinition()->GetPDGEncoding()
79  << ", Z = " << theNucleus.GetZ_asInt())
80  << ", N = " << theNucleus.GetN_asInt())
81  << ", A = " << theNucleus.GetA_asInt() << G4endl;
82  #endif
83 
84  std::pair<G4double,G4double> ratios;
85  ratios=theQuasiElastic->GetRatios(thePrimary.GetTotalMomentum(),
86  thePrimary.GetDefinition()->GetPDGEncoding(),
87  theNucleus.GetZ_asInt(),
88  theNucleus.GetN_asInt());
89  #ifdef debug_scatter
90  G4cout << "G4QuasiElasticChannel::ratios " << ratios.first << " x " <<ratios.second
91  << " = " << ratios.first*ratios.second << G4endl;
92  #endif
93 
94  return ratios.first*ratios.second;
95 }
std::pair< G4double, G4double > GetRatios(G4double pIU, G4int prPDG, G4int tgZ, G4int tgN)
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4ParticleDefinition * GetDefinition() const
G4double GetTotalMomentum() const
G4GLOB_DLL std::ostream G4cout
G4int GetN_asInt() const
Definition: G4Nucleus.hh:112
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

G4KineticTrackVector * G4QuasiElasticChannel::Scatter ( G4Nucleus theNucleus,
const G4DynamicParticle thePrimary 
)

Definition at line 97 of file G4QuasiElasticChannel.cc.

99 {
100  G4int A=theNucleus.GetA_asInt();
101  G4int Z=theNucleus.GetZ_asInt();
102  // build Nucleus and choose random nucleon to scatter with
103  the3DNucleus->Init(theNucleus.GetA_asInt(),theNucleus.GetZ_asInt());
104  const std::vector<G4Nucleon>& nucleons=the3DNucleus->GetNucleons();
105  G4double targetNucleusMass=the3DNucleus->GetMass();
106  G4LorentzVector targetNucleus4Mom(0.,0.,0.,targetNucleusMass);
107  G4int index;
108  do {
109  index=G4lrint((A-1)*G4UniformRand());
110  } while (index < 0 || index >= static_cast<G4int>(nucleons.size())); /* Loop checking, 07.08.2015, A.Ribon */
111 
112  const G4ParticleDefinition * pDef= nucleons[index].GetDefinition();
113 
114  G4int resA=A - 1;
115  G4int resZ=Z - static_cast<int>(pDef->GetPDGCharge());
116  const G4ParticleDefinition* resDef;
117  G4double residualNucleusMass;
118  if(resZ)
119  {
120  resDef=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(resZ,resA,0);
121  residualNucleusMass=resDef->GetPDGMass();
122  }
123  else {
124  resDef=G4Neutron::Neutron();
125  residualNucleusMass=resA * G4Neutron::Neutron()->GetPDGMass();
126  }
127  #ifdef debug_scatter
128  G4cout<<"G4QElChan::Scatter: neutron - proton? A ="<<A<<", Z="<<Z<<", projName="
129  <<pDef->GetParticleName()<<G4endl;
130  #endif
131 
132  G4LorentzVector pNucleon=nucleons[index].Get4Momentum();
133  G4double residualNucleusEnergy=std::sqrt(sqr(residualNucleusMass) +
134  pNucleon.vect().mag2());
135  pNucleon.setE(targetNucleusMass-residualNucleusEnergy);
136  G4LorentzVector residualNucleus4Mom=targetNucleus4Mom-pNucleon;
137 
138  std::pair<G4LorentzVector,G4LorentzVector> result;
139 
140  result=theQuasiElastic->Scatter(pDef->GetPDGEncoding(),pNucleon,
141  thePrimary.GetDefinition()->GetPDGEncoding(),
142  thePrimary.Get4Momentum());
143  G4LorentzVector scatteredHadron4Mom;
144  if (result.first.e() > 0.)
145  scatteredHadron4Mom=result.second;
146  else { //scatter failed
147  //G4cout << "Warning - G4QuasiElasticChannel::Scatter no scattering" << G4endl;
148  //return 0; //no scatter
149  scatteredHadron4Mom=thePrimary.Get4Momentum();
150  residualNucleus4Mom=G4LorentzVector(0.,0.,0.,targetNucleusMass);
152  }
153 
154 #ifdef debug_scatter
155  G4LorentzVector EpConservation=pNucleon+thePrimary.Get4Momentum()
156  - result.first - result.second;
157  if ( (EpConservation.vect().mag2() > .01*MeV*MeV )
158  || (std::abs(EpConservation.e()) > 0.1 * MeV ) )
159  {
160  G4cout << "Warning - G4QuasiElasticChannel::Scatter E-p non conservation : "
161  << EpConservation << G4endl;
162  }
163 #endif
164 
166  G4KineticTrack * sPrim=new G4KineticTrack(thePrimary.GetDefinition(),
167  0.,G4ThreeVector(0), scatteredHadron4Mom);
168  ktv->push_back(sPrim);
169  if (result.first.e() > 0.)
170  {
171  G4KineticTrack * sNuc=new G4KineticTrack(pDef, 0.,G4ThreeVector(0), result.first);
172  ktv->push_back(sNuc);
173  }
174  if(resZ || resA==1) // For the only neutron or for tnuclei with Z>0
175  {
176  G4KineticTrack * rNuc=new G4KineticTrack(resDef,
177  0.,G4ThreeVector(0), residualNucleus4Mom);
178  ktv->push_back(rNuc);
179  }
180  else // The residual nucleus consists of only neutrons
181  {
182  residualNucleus4Mom/=resA; // Split 4-mom of A*n system equally
183  for(G4int in=0; in<resA; in++) // Loop over neutrons in A*n system.
184  {
185  G4KineticTrack* rNuc=new G4KineticTrack(resDef,
186  0.,G4ThreeVector(0), residualNucleus4Mom);
187  ktv->push_back(rNuc);
188  }
189  }
190 #ifdef debug_scatter
191  G4cout<<"G4QElC::Scat: Nucleon: "<<result.first <<" mass "<<result.first.mag() << G4endl;
192  G4cout<<"G4QElC::Scat: Project: "<<result.second<<" mass "<<result.second.mag()<< G4endl;
193 #endif
194  return ktv;
195 }
G4double G4ParticleHPJENDLHEData::G4double result
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
CLHEP::Hep3Vector G4ThreeVector
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
std::pair< G4LorentzVector, G4LorentzVector > Scatter(G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
const G4String & GetParticleName() const
G4IonTable * GetIonTable() const
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
virtual void Init(G4int theA, G4int theZ)=0
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4LorentzVector Get4Momentum() const
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
int G4lrint(double ad)
Definition: templates.hh:163
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
double mag2() const
virtual G4double GetMass()=0
virtual const std::vector< G4Nucleon > & GetNucleons()=0
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
CLHEP::HepLorentzVector G4LorentzVector

Here is the call graph for this function:

Here is the caller graph for this function:


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