Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ElectroNuclearReaction.hh
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 // Class description:
29 // G4ElectroNuclearReaction = eA interface to use CHIPS in the G4Hadronic frame
30 // Created: J.P. Wellisch, following M. Kossov's algorithm. 12/11/2001
31 // The last update: J.P. Wellisch, 06-June-02
32 // 17.02.2009 M.Kossov, now it is recommended to use the G4QCollision process
33 // 10.11.2010 V.Ivanchenko use cross sections by pointer and not by value
34 // 07.09.2011 M.Kelsey, follow changes to G4HadFinalState interface
35 
36 #ifndef G4ElectroNuclearReaction_h
37 #define G4ElectroNuclearReaction_h 1
38 
39 #include <iostream>
40 #include <CLHEP/Units/PhysicalConstants.h>
41 
42 #include "globals.hh"
43 #include "G4HadronicInteraction.hh"
44 #include "G4ChiralInvariantPhaseSpace.hh"
45 #include "G4ElectroNuclearCrossSection.hh"
46 #include "G4PhotoNuclearCrossSection.hh"
47 #include "G4GammaParticipants.hh"
48 #include "G4QGSModel.hh"
49 #include "G4QGSMFragmentation.hh"
50 #include "G4Nucleus.hh"
51 #include "G4HadFinalState.hh"
52 #include "G4HadProjectile.hh"
53 #include "G4Electron.hh"
54 #include "G4Positron.hh"
55 #include "G4Gamma.hh"
56 #include "G4TheoFSGenerator.hh"
57 #include "G4GeneratorPrecompoundInterface.hh"
58 #include "G4ExcitedStringDecay.hh"
59 
60 
61 class G4ElectroNuclearReaction : public G4HadronicInteraction
62 {
63 public:
64  G4ElectroNuclearReaction():G4HadronicInteraction("CHIPSElectroNuclear")
65  {
66  SetMinEnergy(0*CLHEP::GeV);
67  SetMaxEnergy(100*CLHEP::TeV);
68 
69  theHEModel = new G4TheoFSGenerator;
70  theCascade = new G4GeneratorPrecompoundInterface;
71  theHEModel->SetTransport(theCascade);
72  theHEModel->SetHighEnergyGenerator(&theStringModel);
73  theStringDecay = new G4ExcitedStringDecay(&theFragmentation);
74  theStringModel.SetFragmentationModel(theStringDecay);
75  theHEModel->SetMinEnergy(2.5*CLHEP::GeV);
76  theHEModel->SetMaxEnergy(100*CLHEP::TeV);
77  theElectronData = new G4ElectroNuclearCrossSection;
78  thePhotonData = new G4PhotoNuclearCrossSection;
79  Description();
80  }
81 
82  ~G4ElectroNuclearReaction()
83  {
84  delete theHEModel;
85  delete theCascade;
86  delete theStringDecay;
87  delete theElectronData;
88  delete thePhotonData;
89  }
90 
91  G4HadFinalState*
92  ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aTargetNucleus);
93 
94  void Description() const
95  {
96  char* dirName = getenv("G4PhysListDocDir");
97  if (dirName) {
98  std::ofstream outFile;
99  G4String outFileName = GetModelName() + ".html";
100  G4String pathName = G4String(dirName) + "/" + outFileName;
101 
102  outFile.open(pathName);
103  outFile << "<html>\n";
104  outFile << "<head>\n";
105 
106  outFile << "<title>Description of CHIPS ElectroNuclear Model</title>\n";
107  outFile << "</head>\n";
108  outFile << "<body>\n";
109 
110  outFile << "G4ElectroNuclearReaction handles the inelastic scattering\n"
111  << "of e- and e+ from nuclei using the Chiral Invariant Phase\n"
112  << "Space (CHIPS) model of M. Kossov. This model uses the\n"
113  << "Equivalent Photon Approximation in which the incoming\n"
114  << "electron generates a virtual photon at the electromagnetic\n"
115  << "vertex, and the virtual photon is converted to a real photon\n"
116  << "before it interacts with the nucleus. The real photon\n"
117  << "interacts with the hadrons in the target by producing\n"
118  << "quasmons (or generalized excited hadrons) which then decay\n"
119  << "into final state hadrons. This model is valid for e- and\n"
120  << "e+ of all incident energies.\n";
121 
122  outFile << "</body>\n";
123  outFile << "</html>\n";
124  outFile.close();
125  }
126  }
127 
128 private:
129 
130  G4ChiralInvariantPhaseSpace theLEModel;
131  G4TheoFSGenerator* theHEModel;
132  G4GeneratorPrecompoundInterface* theCascade;
133  G4QGSModel<G4GammaParticipants> theStringModel;
134  G4QGSMFragmentation theFragmentation;
135  G4ExcitedStringDecay* theStringDecay;
136  G4ElectroNuclearCrossSection* theElectronData;
137  G4PhotoNuclearCrossSection* thePhotonData;
138  G4HadFinalState theResult;
139 };
140 
141 inline
142 G4HadFinalState* G4ElectroNuclearReaction::ApplyYourself(const G4HadProjectile& aTrack,
143  G4Nucleus& aTargetNucleus)
144 {
145  theResult.Clear();
146  static const G4double dM=G4Proton::Proton()->GetPDGMass() +
147  G4Neutron::Neutron()->GetPDGMass(); // MeanDoubleNucleon Mass = m_n+m_p (@@ no binding)
148  static const G4double me=G4Electron::Electron()->GetPDGMass(); // electron mass
149  static const G4double me2=me*me; // squared electron mass
150  G4DynamicParticle theTempEl(const_cast<G4ParticleDefinition *>(aTrack.GetDefinition()),
151  aTrack.Get4Momentum().vect());
152  const G4DynamicParticle* theElectron=&theTempEl;
153  const G4ParticleDefinition* aD = theElectron->GetDefinition();
154  if((aD != G4Electron::ElectronDefinition()) && (aD != G4Positron::PositronDefinition()))
155  throw G4HadronicException(__FILE__, __LINE__,
156  "G4ElectroNuclearReaction::ApplyYourself called for neither electron or positron");
157  const G4ElementTable* aTab = G4Element::GetElementTable();
158  G4Element * anElement = 0;
159  G4int aZ = static_cast<G4int>(aTargetNucleus.GetZ_asInt()+.1);
160  for(size_t ii=0; ii<aTab->size(); ++ii) if ( std::abs((*aTab)[ii]->GetZ()-aZ) < .1)
161  {
162  anElement = (*aTab)[ii];
163  break;
164  }
165  if(0==anElement)
166  {
167  G4cerr<<"***G4ElectroNuclearReaction::ApplyYourself: element with Z="
168  <<aTargetNucleus.GetZ_asInt()<<" is not in the element table"<<G4endl;
169  throw G4HadronicException(__FILE__, __LINE__, "Anomalous element error.");
170  }
171 
172  // Note: high energy gamma nuclear now implemented.
173  G4double xSec = theElectronData->GetCrossSection(theElectron, anElement);// Check out XS
174  if(xSec<=0.)
175  {
176  theResult.SetStatusChange(isAlive);
177  theResult.SetEnergyChange(theElectron->GetKineticEnergy());
178  // new direction for the electron
179  theResult.SetMomentumChange(theElectron->GetMomentumDirection());
180  return &theResult; // DO-NOTHING condition
181  }
182  G4double photonEnergy = theElectronData->GetEquivalentPhotonEnergy();
183  G4double theElectronKinEnergy=theElectron->GetKineticEnergy();
184  if( theElectronKinEnergy < photonEnergy )
185  {
186  G4cout<<"G4ElectroNuclearReaction::ApplyYourself: photonEnergy is very high"<<G4endl;
187  G4cout<<">>> If this condition persists, please contact Geant4 group"<<G4endl;
188  theResult.SetStatusChange(isAlive);
189  theResult.SetEnergyChange(theElectron->GetKineticEnergy());
190  // new direction for the electron
191  theResult.SetMomentumChange(theElectron->GetMomentumDirection());
192  return &theResult; // DO-NOTHING condition
193  }
194  G4double photonQ2 = theElectronData->GetEquivalentPhotonQ2(photonEnergy);
195  G4double W=photonEnergy-photonQ2/dM; // Hadronic energy flow from the virtual photon
196  if(getenv("debug_G4ElectroNuclearReaction") )
197  {
198  G4cout << "G4ElectroNuclearReaction: Equivalent Energy = "<<W<<G4endl;
199  }
200  if(W<0.)
201  {
202  theResult.SetStatusChange(isAlive);
203  theResult.SetEnergyChange(theElectron->GetKineticEnergy());
204  // new direction for the electron
205  theResult.SetMomentumChange(theElectron->GetMomentumDirection());
206  return &theResult; // DO-NOTHING condition
207  // throw G4HadronicException(__FILE__, __LINE__,
208  // "G4ElectroNuclearReaction::ApplyYourself: negative equivalent energy");
209  }
210  G4DynamicParticle* theDynamicPhoton = new
211  G4DynamicParticle(G4Gamma::GammaDefinition(),
212  G4ParticleMomentum(1.,0.,0.), photonEnergy*CLHEP::MeV); //----->-*
213  G4double sigNu=thePhotonData->GetCrossSection(theDynamicPhoton, anElement); // |
214  theDynamicPhoton->SetKineticEnergy(W); // Redefine photon with equivalent energy |
215  G4double sigK =thePhotonData->GetCrossSection(theDynamicPhoton, anElement); // |
216  delete theDynamicPhoton; // <--------------------------------------------------------*
217  G4double rndFraction = theElectronData->GetVirtualFactor(photonEnergy, photonQ2);
218  if(sigNu*G4UniformRand()>sigK*rndFraction)
219  {
220  theResult.SetStatusChange(isAlive);
221  theResult.SetEnergyChange(theElectron->GetKineticEnergy());
222  // new direction for the electron
223  theResult.SetMomentumChange(theElectron->GetMomentumDirection());
224  return &theResult; // DO-NOTHING condition
225  }
226  theResult.SetStatusChange(isAlive);
227  // Scatter an electron and make gamma+A reaction
228  G4double iniE=theElectronKinEnergy+me; // Initial total energy of electron
229  G4double finE=iniE-photonEnergy; // Final total energy of electron
230  theResult.SetEnergyChange(std::max(0.,finE-me)); // Modifies the KINETIC ENERGY
231  G4double EEm=iniE*finE-me2; // Just an intermediate value to avoid "2*"
232  G4double iniP=std::sqrt(iniE*iniE-me2); // Initial momentum of the electron
233  G4double finP=std::sqrt(finE*finE-me2); // Final momentum of the electron
234  G4double cost=(EEm+EEm-photonQ2)/iniP/finP;// std::cos(theta) for the electron scattering
235  if(cost> 1.) cost= 1.;
236  if(cost<-1.) cost=-1.;
237  G4ThreeVector dir=theElectron->GetMomentumDirection(); // Direction of primary electron
238  G4ThreeVector ort=dir.orthogonal(); // Not normed orthogonal vector (!) (to dir)
239  G4ThreeVector ortx = ort.unit(); // First unit vector orthogonal to the direction
240  G4ThreeVector orty = dir.cross(ortx); // Second unit vector orthoganal to the direction
241  G4double sint=std::sqrt(1.-cost*cost); // Perpendicular component
242  G4double phi=CLHEP::twopi*G4UniformRand(); // phi of scattered electron
243  G4double sinx=sint*std::sin(phi); // x-component
244  G4double siny=sint*std::cos(phi); // y-component
245  G4ThreeVector findir=cost*dir+sinx*ortx+siny*orty;
246  theResult.SetMomentumChange(findir); // new direction for the electron
247  G4ThreeVector photonMomentum=iniP*dir-finP*findir;
248  G4DynamicParticle localGamma(G4Gamma::GammaDefinition(), photonEnergy, photonMomentum);
249  //G4DynamicParticle localGamma(G4Gamma::GammaDefinition(),photonDirection, photonEnergy);
250  //G4DynamicParticle localGamma(G4Gamma::GammaDefinition(), photonLorentzVector);
251  G4ThreeVector position(0,0,0);
252  G4HadProjectile localTrack(localGamma);
253  G4HadFinalState * result;
254  if (photonEnergy < 3*CLHEP::GeV)
255  result = theLEModel.ApplyYourself(localTrack, aTargetNucleus, &theResult);
256  else {
257  // G4cout << "0) Getting a high energy electro-nuclear reaction"<<G4endl;
258  G4HadFinalState * aResult = theHEModel->ApplyYourself(localTrack, aTargetNucleus);
259  theResult.AddSecondaries(aResult);
260  aResult->Clear();
261  result = &theResult;
262  }
263  return result;
264 }
265 
266 #endif
Generated on Sat May 25 2013 14:33:50 for Geant4 by   doxygen 1.8.4