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