Geant4  10.00.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 70688 2013-06-04 09:01:01Z 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(G4QuasiElRatios::GetPointer()),
64  the3DNucleus(new G4Fancy3DNucleus) {}
65 
67 {
68  delete the3DNucleus;
69 }
70 
72  const G4DynamicParticle & thePrimary)
73 {
74  #ifdef debug_scatter
75  G4cout << "G4QuasiElasticChannel:: P=" << thePrimary.GetTotalMomentum()
76  << ", pPDG=" << thePrimary.GetDefinition()->GetPDGEncoding()
77  << ", Z = " << theNucleus.GetZ_asInt())
78  << ", N = " << theNucleus.GetN_asInt())
79  << ", A = " << theNucleus.GetA_asInt() << G4endl;
80  #endif
81 
82  std::pair<G4double,G4double> ratios;
83  ratios=theQuasiElastic->GetRatios(thePrimary.GetTotalMomentum(),
84  thePrimary.GetDefinition()->GetPDGEncoding(),
85  theNucleus.GetZ_asInt(),
86  theNucleus.GetN_asInt());
87  #ifdef debug_scatter
88  G4cout << "G4QuasiElasticChannel::ratios " << ratios.first << " x " <<ratios.second
89  << " = " << ratios.first*ratios.second << G4endl;
90  #endif
91 
92  return ratios.first*ratios.second;
93 }
94 
96  const G4DynamicParticle & thePrimary)
97 {
98  G4int A=theNucleus.GetA_asInt();
99  G4int Z=theNucleus.GetZ_asInt();
100  // build Nucleus and choose random nucleon to scatter with
101  the3DNucleus->Init(theNucleus.GetA_asInt(),theNucleus.GetZ_asInt());
102  const std::vector<G4Nucleon>& nucleons=the3DNucleus->GetNucleons();
103  G4double targetNucleusMass=the3DNucleus->GetMass();
104  G4LorentzVector targetNucleus4Mom(0.,0.,0.,targetNucleusMass);
105  G4int index;
106  do {
107  index=G4lrint((A-1)*G4UniformRand());
108  } while (index < 0 || index >= static_cast<G4int>(nucleons.size()));
109 
110  G4ParticleDefinition * pDef= nucleons[index].GetDefinition();
111 
112  G4int resA=A - 1;
113  G4int resZ=Z - static_cast<int>(pDef->GetPDGCharge());
114  G4ParticleDefinition* resDef;
115  G4double residualNucleusMass;
116  if(resZ)
117  {
118  resDef=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(resZ,resA,0);
119  residualNucleusMass=resDef->GetPDGMass();
120  }
121  else {
122  resDef=G4Neutron::Neutron();
123  residualNucleusMass=resA * G4Neutron::Neutron()->GetPDGMass();
124  }
125  #ifdef debug_scatter
126  G4cout<<"G4QElChan::Scatter: neutron - proton? A ="<<A<<", Z="<<Z<<", projName="
127  <<pDef->GetParticleName()<<G4endl;
128  #endif
129 
130  G4LorentzVector pNucleon=nucleons[index].Get4Momentum();
131  G4double residualNucleusEnergy=std::sqrt(sqr(residualNucleusMass) +
132  pNucleon.vect().mag2());
133  pNucleon.setE(targetNucleusMass-residualNucleusEnergy);
134  G4LorentzVector residualNucleus4Mom=targetNucleus4Mom-pNucleon;
135 
136  std::pair<G4LorentzVector,G4LorentzVector> result;
137 
138  result=theQuasiElastic->Scatter(pDef->GetPDGEncoding(),pNucleon,
139  thePrimary.GetDefinition()->GetPDGEncoding(),
140  thePrimary.Get4Momentum());
141  G4LorentzVector scatteredHadron4Mom;
142  if (result.first.e() > 0.)
143  scatteredHadron4Mom=result.second;
144  else { //scatter failed
145  //G4cout << "Warning - G4QuasiElasticChannel::Scatter no scattering" << G4endl;
146  //return 0; //no scatter
147  scatteredHadron4Mom=thePrimary.Get4Momentum();
148  residualNucleus4Mom=G4LorentzVector(0.,0.,0.,targetNucleusMass);
150  }
151 
152 #ifdef debug_scatter
153  G4LorentzVector EpConservation=pNucleon+thePrimary.Get4Momentum()
154  - result.first - result.second;
155  if ( (EpConservation.vect().mag2() > .01*MeV*MeV )
156  || (std::abs(EpConservation.e()) > 0.1 * MeV ) )
157  {
158  G4cout << "Warning - G4QuasiElasticChannel::Scatter E-p non conservation : "
159  << EpConservation << G4endl;
160  }
161 #endif
162 
164  G4KineticTrack * sPrim=new G4KineticTrack(thePrimary.GetDefinition(),
165  0.,G4ThreeVector(0), scatteredHadron4Mom);
166  ktv->push_back(sPrim);
167  if (result.first.e() > 0.)
168  {
169  G4KineticTrack * sNuc=new G4KineticTrack(pDef, 0.,G4ThreeVector(0), result.first);
170  ktv->push_back(sNuc);
171  }
172  if(resZ || resA==1) // For the only neutron or for tnuclei with Z>0
173  {
174  G4KineticTrack * rNuc=new G4KineticTrack(resDef,
175  0.,G4ThreeVector(0), residualNucleus4Mom);
176  ktv->push_back(rNuc);
177  }
178  else // The residual nucleus consists of only neutrons
179  {
180  residualNucleus4Mom/=resA; // Split 4-mom of A*n system equally
181  for(G4int in=0; in<resA; in++) // Loop over neutrons in A*n system.
182  {
183  G4KineticTrack* rNuc=new G4KineticTrack(resDef,
184  0.,G4ThreeVector(0), residualNucleus4Mom);
185  ktv->push_back(rNuc);
186  }
187  }
188 #ifdef debug_scatter
189  G4cout<<"G4QElC::Scat: Nucleon: "<<result.first <<" mass "<<result.first.mag() << G4endl;
190  G4cout<<"G4QElC::Scat: Project: "<<result.second<<" mass "<<result.second.mag()<< G4endl;
191 #endif
192  return ktv;
193 }
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:193
CLHEP::Hep3Vector G4ThreeVector
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:449
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:87
G4GLOB_DLL std::ostream G4cout
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)
static const G4double A[nN]
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