ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QuasiElasticChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4QuasiElasticChannel.cc
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 //
28 
29 // Author : Gunter Folger March 2007
30 // Modified by Mikhail Kossov. Apr2009, E/M conservation: ResidualNucleus is added (ResNuc)
31 // Class Description
32 // Final state production model for theoretical models of hadron inelastic
33 // quasi elastic scattering in geant4;
34 // Class Description - End
35 //
36 // Modified:
37 // 20110805 M. Kelsey -- Follow change to G4V3DNucleus::GetNucleons()
38 // 20110808 M. Kelsey -- Move #includes from .hh, add many missing
39 
40 #include "G4QuasiElasticChannel.hh"
41 
42 #include "G4Fancy3DNucleus.hh"
43 #include "G4DynamicParticle.hh"
44 #include "G4HadTmpUtil.hh" /* lrint */
45 #include "G4KineticTrack.hh"
46 #include "G4KineticTrackVector.hh"
47 #include "G4LorentzVector.hh"
48 #include "G4Neutron.hh"
49 #include "G4Nucleon.hh"
50 #include "G4Nucleus.hh"
51 #include "G4ParticleDefinition.hh"
52 #include "G4ParticleTable.hh"
53 #include "G4IonTable.hh"
54 #include "G4QuasiElRatios.hh"
55 #include "globals.hh"
56 #include <vector>
57 
58 //#define debug_scatter
59 
60 
62  : theQuasiElastic(new G4QuasiElRatios),
63  the3DNucleus(new G4Fancy3DNucleus) {
64 }
65 
67 {
68  delete the3DNucleus;
69  delete theQuasiElastic;
70 }
71 
73  const G4DynamicParticle & thePrimary)
74 {
75  #ifdef debug_scatter
76  G4cout << "G4QuasiElasticChannel:: P=" << thePrimary.GetTotalMomentum()
77  << ", pPDG=" << thePrimary.GetDefinition()->GetPDGEncoding()
78  << ", Z = " << theNucleus.GetZ_asInt())
79  << ", N = " << theNucleus.GetN_asInt())
80  << ", A = " << theNucleus.GetA_asInt() << G4endl;
81  #endif
82 
83  std::pair<G4double,G4double> ratios;
84  ratios=theQuasiElastic->GetRatios(thePrimary.GetTotalMomentum(),
85  thePrimary.GetDefinition()->GetPDGEncoding(),
86  theNucleus.GetZ_asInt(),
87  theNucleus.GetN_asInt());
88  #ifdef debug_scatter
89  G4cout << "G4QuasiElasticChannel::ratios " << ratios.first << " x " <<ratios.second
90  << " = " << ratios.first*ratios.second << G4endl;
91  #endif
92 
93  return ratios.first*ratios.second;
94 }
95 
97  const G4DynamicParticle & thePrimary)
98 {
99  G4int A=theNucleus.GetA_asInt();
100  G4int Z=theNucleus.GetZ_asInt();
101  // build Nucleus and choose random nucleon to scatter with
102  the3DNucleus->Init(theNucleus.GetA_asInt(),theNucleus.GetZ_asInt());
103  const std::vector<G4Nucleon>& nucleons=the3DNucleus->GetNucleons();
104  G4double targetNucleusMass=the3DNucleus->GetMass();
105  G4LorentzVector targetNucleus4Mom(0.,0.,0.,targetNucleusMass);
106  G4int index;
107  do {
108  index=G4lrint((A-1)*G4UniformRand());
109  } while (index < 0 || index >= static_cast<G4int>(nucleons.size())); /* Loop checking, 07.08.2015, A.Ribon */
110 
111  const G4ParticleDefinition * pDef= nucleons[index].GetDefinition();
112 
113  G4int resA=A - 1;
114  G4int resZ=Z - static_cast<int>(pDef->GetPDGCharge());
115  const G4ParticleDefinition* resDef;
116  G4double residualNucleusMass;
117  if(resZ)
118  {
119  resDef=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(resZ,resA,0);
120  residualNucleusMass=resDef->GetPDGMass();
121  }
122  else {
123  resDef=G4Neutron::Neutron();
124  residualNucleusMass=resA * G4Neutron::Neutron()->GetPDGMass();
125  }
126  #ifdef debug_scatter
127  G4cout<<"G4QElChan::Scatter: neutron - proton? A ="<<A<<", Z="<<Z<<", projName="
128  <<pDef->GetParticleName()<<G4endl;
129  #endif
130 
131  G4LorentzVector pNucleon=nucleons[index].Get4Momentum();
132  G4double residualNucleusEnergy=std::sqrt(sqr(residualNucleusMass) +
133  pNucleon.vect().mag2());
134  pNucleon.setE(targetNucleusMass-residualNucleusEnergy);
135  G4LorentzVector residualNucleus4Mom=targetNucleus4Mom-pNucleon;
136 
137  std::pair<G4LorentzVector,G4LorentzVector> result;
138 
139  result=theQuasiElastic->Scatter(pDef->GetPDGEncoding(),pNucleon,
140  thePrimary.GetDefinition()->GetPDGEncoding(),
141  thePrimary.Get4Momentum());
142  G4LorentzVector scatteredHadron4Mom;
143  if (result.first.e() > 0.)
144  scatteredHadron4Mom=result.second;
145  else { //scatter failed
146  //G4cout << "Warning - G4QuasiElasticChannel::Scatter no scattering" << G4endl;
147  //return 0; //no scatter
148  scatteredHadron4Mom=thePrimary.Get4Momentum();
149  residualNucleus4Mom=G4LorentzVector(0.,0.,0.,targetNucleusMass);
151  }
152 
153 #ifdef debug_scatter
154  G4LorentzVector EpConservation=pNucleon+thePrimary.Get4Momentum()
155  - result.first - result.second;
156  if ( (EpConservation.vect().mag2() > .01*MeV*MeV )
157  || (std::abs(EpConservation.e()) > 0.1 * MeV ) )
158  {
159  G4cout << "Warning - G4QuasiElasticChannel::Scatter E-p non conservation : "
160  << EpConservation << G4endl;
161  }
162 #endif
163 
165  G4KineticTrack * sPrim=new G4KineticTrack(thePrimary.GetDefinition(),
166  0.,G4ThreeVector(0), scatteredHadron4Mom);
167  ktv->push_back(sPrim);
168  if (result.first.e() > 0.)
169  {
170  G4KineticTrack * sNuc=new G4KineticTrack(pDef, 0.,G4ThreeVector(0), result.first);
171  ktv->push_back(sNuc);
172  }
173  if(resZ || resA==1) // For the only neutron or for tnuclei with Z>0
174  {
175  G4KineticTrack * rNuc=new G4KineticTrack(resDef,
176  0.,G4ThreeVector(0), residualNucleus4Mom);
177  ktv->push_back(rNuc);
178  }
179  else // The residual nucleus consists of only neutrons
180  {
181  residualNucleus4Mom/=resA; // Split 4-mom of A*n system equally
182  for(G4int in=0; in<resA; in++) // Loop over neutrons in A*n system.
183  {
184  G4KineticTrack* rNuc=new G4KineticTrack(resDef,
185  0.,G4ThreeVector(0), residualNucleus4Mom);
186  ktv->push_back(rNuc);
187  }
188  }
189 #ifdef debug_scatter
190  G4cout<<"G4QElC::Scat: Nucleon: "<<result.first <<" mass "<<result.first.mag() << G4endl;
191  G4cout<<"G4QElC::Scat: Project: "<<result.second<<" mass "<<result.second.mag()<< G4endl;
192 #endif
193  return ktv;
194 }