Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4ElectroVDNuclearModel.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 // $Id: $
27 //
28 // Author: D.H. Wright
29 // Date: 1 May 2012
30 //
31 // Description: model for electron and positron interaction with nuclei
32 // using the equivalent photon spectrum. A real gamma is
33 // produced from the virtual photon spectrum and is then
34 // interacted hadronically by the Bertini cascade at low
35 // energies. At high energies the gamma is treated as a
36 // pi0 and interacted with the nucleus using the FTFP model.
37 // The electro- and photo-nuclear cross sections of
38 // M. Kossov are used to generate the virtual photon
39 // spectrum.
40 //
41 
43 
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 
50 
51 #include "G4CascadeInterface.hh"
52 #include "G4TheoFSGenerator.hh"
54 #include "G4ExcitationHandler.hh"
55 #include "G4PreCompoundModel.hh"
57 #include "G4ExcitedStringDecay.hh"
58 #include "G4FTFModel.hh"
59 
60 #include "G4HadFinalState.hh"
62 
64  : G4HadronicInteraction("G4ElectroVDNuclearModel"),
65  leptonKE(0.0), photonEnergy(0.0), photonQ2(0.0)
66 {
67  SetMinEnergy(0.0);
68  SetMaxEnergy(1*PeV);
69  electroXS =
71  GetCrossSectionDataSet(G4ElectroNuclearCrossSection::Default_Name());
72  gammaXS =
74  GetCrossSectionDataSet(G4PhotoNuclearCrossSection::Default_Name());
75 
76  // reuse existing pre-compound model
77  G4GeneratorPrecompoundInterface* precoInterface
81  G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
82  if(!pre) { pre = new G4PreCompoundModel(); }
83  precoInterface->SetDeExcitation(pre);
84 
85  // string model
86  ftfp = new G4TheoFSGenerator();
87  ftfp->SetTransport(precoInterface);
88  theFragmentation = new G4LundStringFragmentation();
89  theStringDecay = new G4ExcitedStringDecay(theFragmentation);
90  G4FTFModel* theStringModel = new G4FTFModel();
91  theStringModel->SetFragmentationModel(theStringDecay);
92  ftfp->SetHighEnergyGenerator(theStringModel);
93 
94  // Build Bertini model
95  bert = new G4CascadeInterface();
96 }
97 
99 {
100  delete theFragmentation;
101  delete theStringDecay;
102 }
103 
104 void G4ElectroVDNuclearModel::ModelDescription(std::ostream& outFile) const
105 {
106  outFile << "G4ElectroVDNuclearModel handles the inelastic scattering\n"
107  << "of e- and e+ from nuclei using the equivalent photon\n"
108  << "approximation in which the incoming lepton generates a\n"
109  << "virtual photon at the electromagnetic vertex, and the\n"
110  << "virtual photon is converted to a real photon. At low\n"
111  << "energies, the photon interacts directly with the nucleus\n"
112  << "using the Bertini cascade. At high energies the photon\n"
113  << "is converted to a pi0 which interacts using the FTFP\n"
114  << "model. The electro- and gamma-nuclear cross sections of\n"
115  << "M. Kossov are used to generate the virtual photon spectrum\n";
116 }
117 
118 
121  G4Nucleus& targetNucleus)
122 {
123  // Set up default particle change (just returns initial state)
126  leptonKE = aTrack.GetKineticEnergy();
129 
130  // Set up sanity checks for real photon production
131  G4DynamicParticle lepton(aTrack.GetDefinition(), aTrack.Get4Momentum() );
132 
133  // Need to call GetElementCrossSection before calling GetEquivalentPhotonEnergy.
134  G4Material* mat = 0;
135  G4int targZ = targetNucleus.GetZ_asInt();
136  electroXS->GetElementCrossSection(&lepton, targZ, mat);
137 
138  photonEnergy = electroXS->GetEquivalentPhotonEnergy();
139  // Photon energy cannot exceed lepton energy
140  if (photonEnergy < leptonKE) {
141  photonQ2 = electroXS->GetEquivalentPhotonQ2(photonEnergy);
143  // Photon
144  if (photonEnergy > photonQ2/dM) {
145  // Produce recoil lepton and transferred photon
146  G4DynamicParticle* transferredPhoton = CalculateEMVertex(aTrack, targetNucleus);
147  // Interact gamma with nucleus
148  if (transferredPhoton) CalculateHadronicVertex(transferredPhoton, targetNucleus);
149  }
150  }
151  return &theParticleChange;
152 }
153 
154 
156 G4ElectroVDNuclearModel::CalculateEMVertex(const G4HadProjectile& aTrack,
157  G4Nucleus& targetNucleus)
158 {
159  G4DynamicParticle photon(G4Gamma::Gamma(), photonEnergy,
160  G4ThreeVector(0.,0.,1.) );
161 
162  // Get gamma cross section at Q**2 = 0 (real gamma)
163  G4int targZ = targetNucleus.GetZ_asInt();
164  G4Material* mat = 0;
165  G4double sigNu =
166  gammaXS->GetElementCrossSection(&photon, targZ, mat);
167 
168  // Change real gamma energy to equivalent energy and get cross section at that energy
170  photon.SetKineticEnergy(photonEnergy - photonQ2/dM);
171  G4double sigK =
172  gammaXS->GetElementCrossSection(&photon, targZ, mat);
173  G4double rndFraction = electroXS->GetVirtualFactor(photonEnergy, photonQ2);
174 
175  // No gamma produced, return null ptr
176  if (sigNu*G4UniformRand() > sigK*rndFraction) return 0;
177 
178  // Scatter the lepton
179  G4double mProj = aTrack.GetDefinition()->GetPDGMass();
180  G4double mProj2 = mProj*mProj;
181  G4double iniE = leptonKE + mProj; // Total energy of incident lepton
182  G4double finE = iniE - photonEnergy; // Total energy of scattered lepton
184  G4double iniP = std::sqrt(iniE*iniE-mProj2); // Incident lepton momentum
185  G4double finP = std::sqrt(finE*finE-mProj2); // Scattered lepton momentum
186  G4double cost = (iniE*finE - mProj2 - photonQ2/2.)/iniP/finP; // cos(theta) from Q**2
187  if (cost > 1.) cost= 1.;
188  if (cost < -1.) cost=-1.;
189  G4double sint = std::sqrt(1.-cost*cost);
190 
191  G4ThreeVector dir = aTrack.Get4Momentum().vect().unit();
192  G4ThreeVector ortx = dir.orthogonal().unit(); // Ortho-normal to scattering plane
193  G4ThreeVector orty = dir.cross(ortx); // Third unit vector
194  G4double phi = twopi*G4UniformRand();
195  G4double sinx = sint*std::sin(phi);
196  G4double siny = sint*std::cos(phi);
197  G4ThreeVector findir = cost*dir+sinx*ortx+siny*orty;
198  theParticleChange.SetMomentumChange(findir); // change lepton direction
199 
200  // Create a gamma with momentum equal to momentum transfer
201  G4ThreeVector photonMomentum = iniP*dir - finP*findir;
203  photonEnergy, photonMomentum);
204  return gamma;
205 }
206 
207 
208 void
209 G4ElectroVDNuclearModel::CalculateHadronicVertex(G4DynamicParticle* incident,
210  G4Nucleus& target)
211 {
212  G4HadFinalState* hfs = 0;
213  G4double gammaE = incident->GetTotalEnergy();
214 
215  if (gammaE < 10*GeV) {
216  G4HadProjectile projectile(*incident);
217  hfs = bert->ApplyYourself(projectile, target);
218  } else {
219  // At high energies convert incident gamma to a pion
221  G4double piMom = std::sqrt(gammaE*gammaE - piMass*piMass);
222  G4ThreeVector piMomentum(incident->GetMomentumDirection() );
223  piMomentum *= piMom;
224  G4DynamicParticle theHadron(G4PionZero::PionZero(), piMomentum);
225  G4HadProjectile projectile(theHadron);
226  hfs = ftfp->ApplyYourself(projectile, target);
227  }
228 
229  delete incident;
230 
231  // Copy secondaries from sub-model to model
233 }
234 
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
void AddSecondaries(const std::vector< G4HadSecondary > &addSecs)
virtual void ModelDescription(std::ostream &outFile) const
const XML_Char * target
Definition: expat.h:268
CLHEP::Hep3Vector G4ThreeVector
G4double GetTotalEnergy() const
static const G4double dM
void SetFragmentationModel(G4VStringFragmentation *aModel)
const char * p
Definition: xmltok.h:285
int G4int
Definition: G4Types.hh:78
void SetHighEnergyGenerator(G4VHighEnergyGenerator *const value)
static constexpr double twopi
Definition: G4SIunits.hh:76
void SetStatusChange(G4HadFinalStateStatus aS)
void SetMinEnergy(G4double anEnergy)
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
const G4ParticleDefinition * GetDefinition() const
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4CrossSectionDataSetRegistry * Instance()
static G4Proton * Proton()
Definition: G4Proton.cc:93
static constexpr double PeV
Definition: G4SIunits.hh:219
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
const G4LorentzVector & Get4Momentum() const
static G4PionZero * PionZero()
Definition: G4PionZero.cc:108
G4HadronicInteraction * FindModel(const G4String &name)
virtual G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *)
void SetEnergyChange(G4double anEnergy)
G4double GetPDGMass() const
G4double GetVirtualFactor(G4double nu, G4double Q2)
static G4HadronicInteractionRegistry * Instance()
virtual G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *mat)
Hep3Vector unit() const
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
Hep3Vector orthogonal() const
static constexpr double GeV
Definition: G4SIunits.hh:217
void SetMaxEnergy(const G4double anEnergy)
void SetDeExcitation(G4VPreCompoundModel *ptr)
void SetTransport(G4VIntraNuclearTransportModel *const value)
Hep3Vector cross(const Hep3Vector &) const
double G4double
Definition: G4Types.hh:76
void SetMomentumChange(const G4ThreeVector &aV)
G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)