G4ElectroVDNuclearModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: $
00027 //
00028 // Author:  D.H. Wright
00029 // Date:    1 May 2012
00030 //
00031 // Description: model for electron and positron interaction with nuclei
00032 //              using the equivalent photon spectrum.  A real gamma is 
00033 //              produced from the virtual photon spectrum and is then 
00034 //              interacted hadronically by the Bertini cascade at low
00035 //              energies.  At high energies the gamma is treated as a 
00036 //              pi0 and interacted with the nucleus using the FTFP model.
00037 //              The electro- and photo-nuclear cross sections of
00038 //              M. Kossov are used to generate the virtual photon
00039 //              spectrum.
00040 //
00041 
00042 #include "G4ElectroVDNuclearModel.hh"
00043 
00044 #include "G4PhysicalConstants.hh"
00045 #include "G4SystemOfUnits.hh"
00046 
00047 #include "G4ElectroNuclearCrossSection.hh"
00048 #include "G4PhotoNuclearCrossSection.hh"
00049 
00050 #include "G4CascadeInterface.hh"
00051 #include "G4TheoFSGenerator.hh"
00052 #include "G4GeneratorPrecompoundInterface.hh"
00053 #include "G4ExcitationHandler.hh"
00054 #include "G4PreCompoundModel.hh"
00055 #include "G4LundStringFragmentation.hh"
00056 #include "G4ExcitedStringDecay.hh"
00057 #include "G4FTFModel.hh"
00058 
00059 #include "G4HadFinalState.hh"
00060 
00061 
00062 G4ElectroVDNuclearModel::G4ElectroVDNuclearModel()
00063  : G4HadronicInteraction("G4ElectroVDNuclearModel"),
00064    leptonKE(0.0), photonEnergy(0.0), photonQ2(0.0)
00065 {
00066   SetMinEnergy(0.0);
00067   SetMaxEnergy(1*PeV);
00068 
00069   electroXS = new G4ElectroNuclearCrossSection();
00070   gammaXS = new G4PhotoNuclearCrossSection();      
00071   ftfp = new G4TheoFSGenerator();
00072   precoInterface = new G4GeneratorPrecompoundInterface();
00073   theHandler = new G4ExcitationHandler();
00074   preEquilib = new G4PreCompoundModel(theHandler);
00075   precoInterface->SetDeExcitation(preEquilib);
00076   ftfp->SetTransport(precoInterface);
00077   theFragmentation = new G4LundStringFragmentation();
00078   theStringDecay = new G4ExcitedStringDecay(theFragmentation);
00079   theStringModel = new G4FTFModel();
00080   theStringModel->SetFragmentationModel(theStringDecay);
00081   ftfp->SetHighEnergyGenerator(theStringModel);
00082 
00083   // Build Bertini model
00084   bert = new G4CascadeInterface();
00085 }
00086 
00087 
00088 G4ElectroVDNuclearModel::~G4ElectroVDNuclearModel()
00089 {
00090   delete electroXS;
00091   delete gammaXS;
00092   delete ftfp;
00093   delete preEquilib;
00094   delete theFragmentation;
00095   delete theStringDecay;
00096   delete theStringModel;
00097   delete bert;
00098 }
00099 
00100     
00101 void G4ElectroVDNuclearModel::ModelDescription(std::ostream& outFile) const 
00102 {
00103   outFile << "G4ElectroVDNuclearModel handles the inelastic scattering\n"
00104           << "of e- and e+ from nuclei using the equivalent photon\n"
00105           << "approximation in which the incoming lepton generates a\n"
00106           << "virtual photon at the electromagnetic vertex, and the\n"
00107           << "virtual photon is converted to a real photon.  At low\n"
00108           << "energies, the photon interacts directly with the nucleus\n"
00109           << "using the Bertini cascade.  At high energies the photon\n"
00110           << "is converted to a pi0 which interacts using the FTFP\n"
00111           << "model.  The electro- and gamma-nuclear cross sections of\n"
00112           << "M. Kossov are used to generate the virtual photon spectrum\n";
00113 }
00114 
00115 
00116 G4HadFinalState*
00117 G4ElectroVDNuclearModel::ApplyYourself(const G4HadProjectile& aTrack, 
00118                                        G4Nucleus& targetNucleus)
00119 {
00120   // Set up default particle change (just returns initial state)
00121   theParticleChange.Clear();
00122   theParticleChange.SetStatusChange(isAlive);
00123   leptonKE = aTrack.GetKineticEnergy();
00124   theParticleChange.SetEnergyChange(leptonKE);
00125   theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit() );
00126  
00127   // Set up sanity checks for real photon production 
00128   G4DynamicParticle lepton(aTrack.GetDefinition(), aTrack.Get4Momentum() );
00129   G4int targZ = targetNucleus.GetZ_asInt();
00130   G4int targA = targetNucleus.GetA_asInt();
00131   G4Isotope* iso = 0;
00132   G4Element* ele = 0;
00133   G4Material* mat = 0; 
00134   G4double eXS = electroXS->GetIsoCrossSection(&lepton, targZ, targA, iso, ele, mat);
00135 
00136   // If electronuclear cross section is negative, return initial track
00137   if (eXS > 0.0) {
00138     photonEnergy = electroXS->GetEquivalentPhotonEnergy();
00139     // Photon energy cannot exceed lepton energy
00140     if (photonEnergy < leptonKE) {
00141       photonQ2 = electroXS->GetEquivalentPhotonQ2(photonEnergy);
00142       G4double dM = G4Proton::Proton()->GetPDGMass() + G4Neutron::Neutron()->GetPDGMass();
00143       // Photon
00144       if (photonEnergy > photonQ2/dM) {
00145         // Produce recoil lepton and transferred photon
00146         G4DynamicParticle* transferredPhoton = CalculateEMVertex(aTrack, targetNucleus);
00147         // Interact gamma with nucleus
00148         if (transferredPhoton) CalculateHadronicVertex(transferredPhoton, targetNucleus);
00149       }
00150     }
00151   }
00152   return &theParticleChange;
00153 }
00154 
00155 
00156 G4DynamicParticle*
00157 G4ElectroVDNuclearModel::CalculateEMVertex(const G4HadProjectile& aTrack,
00158                                            G4Nucleus& targetNucleus)
00159 {
00160   G4DynamicParticle photon(G4Gamma::Gamma(), photonEnergy,
00161                            G4ThreeVector(0.,0.,1.) );
00162 
00163   // Get gamma cross section at Q**2 = 0 (real gamma)
00164   G4int targZ = targetNucleus.GetZ_asInt();
00165   G4int targA = targetNucleus.GetA_asInt();
00166   G4Isotope* iso = 0;
00167   G4Element* ele = 0;
00168   G4Material* mat = 0;
00169   G4double sigNu =
00170     gammaXS->GetIsoCrossSection(&photon, targZ, targA, iso, ele, mat);
00171 
00172   // Change real gamma energy to equivalent energy and get cross section at that energy 
00173   G4double dM = G4Proton::Proton()->GetPDGMass() + G4Neutron::Neutron()->GetPDGMass();
00174   photon.SetKineticEnergy(photonEnergy - photonQ2/dM);      
00175   G4double sigK =
00176     gammaXS->GetIsoCrossSection(&photon, targZ, targA, iso, ele, mat);
00177   G4double rndFraction = electroXS->GetVirtualFactor(photonEnergy, photonQ2);
00178 
00179   // No gamma produced, return null ptr
00180   if (sigNu*G4UniformRand() > sigK*rndFraction) return 0;
00181 
00182   // Scatter the lepton
00183   G4double mProj = aTrack.GetDefinition()->GetPDGMass();
00184   G4double mProj2 = mProj*mProj;
00185   G4double iniE = leptonKE + mProj;               // Total energy of incident lepton
00186   G4double finE = iniE - photonEnergy;            // Total energy of scattered lepton
00187   theParticleChange.SetEnergyChange(finE-mProj);
00188   G4double iniP = std::sqrt(iniE*iniE-mProj2);    // Incident lepton momentum
00189   G4double finP = std::sqrt(finE*finE-mProj2);    // Scattered lepton momentum
00190   G4double cost = (iniE*finE - mProj2 - photonQ2/2.)/iniP/finP;  // cos(theta) from Q**2
00191   if (cost > 1.) cost= 1.;
00192   if (cost < -1.) cost=-1.;
00193   G4double sint = std::sqrt(1.-cost*cost);
00194 
00195   G4ThreeVector dir = aTrack.Get4Momentum().vect().unit();
00196   G4ThreeVector ortx = dir.orthogonal().unit();   // Ortho-normal to scattering plane
00197   G4ThreeVector orty = dir.cross(ortx);           // Third unit vector
00198   G4double phi = twopi*G4UniformRand();
00199   G4double sinx = sint*std::sin(phi);
00200   G4double siny = sint*std::cos(phi);
00201   G4ThreeVector findir = cost*dir+sinx*ortx+siny*orty;
00202   theParticleChange.SetMomentumChange(findir);    // change lepton direction
00203 
00204   // Create a gamma with momentum equal to momentum transfer
00205   G4ThreeVector photonMomentum = iniP*dir - finP*findir;
00206   G4DynamicParticle* gamma = new G4DynamicParticle(G4Gamma::Gamma(),
00207                                                    photonEnergy, photonMomentum);
00208   return gamma;
00209 }
00210 
00211 
00212 void
00213 G4ElectroVDNuclearModel::CalculateHadronicVertex(G4DynamicParticle* incident,
00214                                                  G4Nucleus& target)
00215 {
00216   G4HadFinalState* hfs = 0;
00217   G4double gammaE = incident->GetTotalEnergy();
00218 
00219   if (gammaE < 10*GeV) {
00220     G4HadProjectile projectile(*incident);
00221     hfs = bert->ApplyYourself(projectile, target);
00222   } else {
00223     // At high energies convert incident gamma to a pion
00224     G4double piMass = G4PionZero::PionZero()->GetPDGMass();
00225     G4double piMom = std::sqrt(gammaE*gammaE - piMass*piMass);
00226     G4ThreeVector piMomentum(incident->GetMomentumDirection() );
00227     piMomentum *= piMom;
00228     G4DynamicParticle theHadron(G4PionZero::PionZero(), piMomentum);
00229     G4HadProjectile projectile(theHadron);
00230     hfs = ftfp->ApplyYourself(projectile, target);
00231   }
00232 
00233   delete incident;
00234 
00235   // Copy secondaries from sub-model to model
00236   theParticleChange.AddSecondaries(hfs);
00237 }
00238 

Generated on Mon May 27 17:48:07 2013 for Geant4 by  doxygen 1.4.7