Geant4  10.02.p01
G4QuasiElasticChannel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4QuasiElasticChannel.cc 91859 2015-08-07 13:58:30Z gcosmo $
28 //
29 
30 // Author : Gunter Folger March 2007
31 // Modified by Mikhail Kossov. Apr2009, E/M conservation: ResidualNucleus is added (ResNuc)
32 // Class Description
33 // Final state production model for theoretical models of hadron inelastic
34 // quasi elastic scattering in geant4;
35 // Class Description - End
36 //
37 // Modified:
38 // 20110805 M. Kelsey -- Follow change to G4V3DNucleus::GetNucleons()
39 // 20110808 M. Kelsey -- Move #includes from .hh, add many missing
40 
41 #include "G4QuasiElasticChannel.hh"
42 
43 #include "G4Fancy3DNucleus.hh"
44 #include "G4DynamicParticle.hh"
45 #include "G4HadTmpUtil.hh" /* lrint */
46 #include "G4KineticTrack.hh"
47 #include "G4KineticTrackVector.hh"
48 #include "G4LorentzVector.hh"
49 #include "G4Neutron.hh"
50 #include "G4Nucleon.hh"
51 #include "G4Nucleus.hh"
52 #include "G4ParticleDefinition.hh"
53 #include "G4ParticleTable.hh"
54 #include "G4IonTable.hh"
55 #include "G4QuasiElRatios.hh"
56 #include "globals.hh"
57 #include <vector>
58 
59 //#define debug_scatter
60 
61 
63  : theQuasiElastic(new G4QuasiElRatios),
64  the3DNucleus(new G4Fancy3DNucleus) {
65 }
66 
68 {
69  delete the3DNucleus;
70  delete theQuasiElastic;
71 }
72 
74  const G4DynamicParticle & thePrimary)
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 }
96 
98  const G4DynamicParticle & thePrimary)
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 }
std::pair< G4double, G4double > GetRatios(G4double pIU, G4int prPDG, G4int tgZ, G4int tgN)
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
static const double MeV
Definition: G4SIunits.hh:211
CLHEP::Hep3Vector G4ThreeVector
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:491
G4QuasiElRatios * theQuasiElastic
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
G4KineticTrackVector * Scatter(G4Nucleus &theNucleus, const G4DynamicParticle &thePrimary)
std::pair< G4LorentzVector, G4LorentzVector > Scatter(G4int NPDG, G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M)
const G4String & GetParticleName() const
G4double GetTotalMomentum() const
G4IonTable * GetIonTable() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4int GetN_asInt() const
Definition: G4Nucleus.hh:112
virtual void Init(G4int theA, G4int theZ)=0
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double GetFraction(G4Nucleus &theNucleus, const G4DynamicParticle &thePrimary)
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
virtual G4double GetMass()=0
virtual const std::vector< G4Nucleon > & GetNucleons()=0
#define G4endl
Definition: G4ios.hh:61
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
CLHEP::HepLorentzVector G4LorentzVector