Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // 27 // Author: D.H. Wright 27 // Author: D.H. Wright 28 // Date: 1 May 2012 28 // Date: 1 May 2012 29 // 29 // 30 // Description: model for electron and positro 30 // Description: model for electron and positron interaction with nuclei 31 // using the equivalent photon sp 31 // using the equivalent photon spectrum. A real gamma is 32 // produced from the virtual phot 32 // produced from the virtual photon spectrum and is then 33 // interacted hadronically by the 33 // interacted hadronically by the Bertini cascade at low 34 // energies. At high energies th 34 // energies. At high energies the gamma is treated as a 35 // pi0 and interacted with the nu 35 // pi0 and interacted with the nucleus using the FTFP model. 36 // The electro- and photo-nuclear 36 // The electro- and photo-nuclear cross sections of 37 // M. Kossov are used to generate 37 // M. Kossov are used to generate the virtual photon 38 // spectrum. 38 // spectrum. 39 // 39 // 40 40 41 #include "G4ElectroVDNuclearModel.hh" 41 #include "G4ElectroVDNuclearModel.hh" 42 42 43 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh" 44 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh" 45 45 46 #include "G4ElectroNuclearCrossSection.hh" 46 #include "G4ElectroNuclearCrossSection.hh" 47 #include "G4PhotoNuclearCrossSection.hh" 47 #include "G4PhotoNuclearCrossSection.hh" 48 #include "G4CrossSectionDataSetRegistry.hh" 48 #include "G4CrossSectionDataSetRegistry.hh" 49 49 50 #include "G4CascadeInterface.hh" 50 #include "G4CascadeInterface.hh" 51 #include "G4TheoFSGenerator.hh" 51 #include "G4TheoFSGenerator.hh" 52 #include "G4GeneratorPrecompoundInterface.hh" 52 #include "G4GeneratorPrecompoundInterface.hh" 53 #include "G4ExcitationHandler.hh" 53 #include "G4ExcitationHandler.hh" 54 #include "G4PreCompoundModel.hh" 54 #include "G4PreCompoundModel.hh" 55 #include "G4LundStringFragmentation.hh" 55 #include "G4LundStringFragmentation.hh" 56 #include "G4ExcitedStringDecay.hh" 56 #include "G4ExcitedStringDecay.hh" 57 #include "G4FTFModel.hh" 57 #include "G4FTFModel.hh" 58 58 59 #include "G4HadFinalState.hh" 59 #include "G4HadFinalState.hh" 60 #include "G4HadronicInteractionRegistry.hh" 60 #include "G4HadronicInteractionRegistry.hh" 61 #include "G4PhysicsModelCatalog.hh" 61 #include "G4PhysicsModelCatalog.hh" 62 62 63 #include "G4ElectroNuclearCrossSection.hh" 63 #include "G4ElectroNuclearCrossSection.hh" 64 #include "G4PhotoNuclearCrossSection.hh" 64 #include "G4PhotoNuclearCrossSection.hh" 65 #include "G4GammaNuclearXS.hh" 65 #include "G4GammaNuclearXS.hh" 66 66 67 67 68 G4ElectroVDNuclearModel::G4ElectroVDNuclearMod 68 G4ElectroVDNuclearModel::G4ElectroVDNuclearModel() 69 : G4HadronicInteraction("G4ElectroVDNuclearMo 69 : G4HadronicInteraction("G4ElectroVDNuclearModel"), 70 leptonKE(0.0), photonEnergy(0.0), photonQ2( 70 leptonKE(0.0), photonEnergy(0.0), photonQ2(0.0), secID(-1) 71 { 71 { 72 SetMinEnergy(0.0); 72 SetMinEnergy(0.0); 73 SetMaxEnergy(1*PeV); 73 SetMaxEnergy(1*PeV); 74 74 75 electroXS = 75 electroXS = 76 (G4ElectroNuclearCrossSection*)G4CrossSect 76 (G4ElectroNuclearCrossSection*)G4CrossSectionDataSetRegistry::Instance()-> 77 GetCrossSectionDataSet(G4ElectroNuclearCro 77 GetCrossSectionDataSet(G4ElectroNuclearCrossSection::Default_Name()); 78 if ( electroXS == nullptr ) { 78 if ( electroXS == nullptr ) { 79 electroXS = new G4ElectroNuclearCrossSecti 79 electroXS = new G4ElectroNuclearCrossSection; 80 } 80 } 81 81 82 gammaXS = 82 gammaXS = 83 (G4PhotoNuclearCrossSection*)G4CrossSectio 83 (G4PhotoNuclearCrossSection*)G4CrossSectionDataSetRegistry::Instance()-> 84 GetCrossSectionDataSet(G4PhotoNuclearCross 84 GetCrossSectionDataSet(G4PhotoNuclearCrossSection::Default_Name()); 85 if ( gammaXS == nullptr ) { 85 if ( gammaXS == nullptr ) { 86 gammaXS = 86 gammaXS = 87 (G4GammaNuclearXS*)G4CrossSectionDataSet 87 (G4GammaNuclearXS*)G4CrossSectionDataSetRegistry::Instance()-> 88 GetCrossSectionDataSet(G4GammaNuclearXS: 88 GetCrossSectionDataSet(G4GammaNuclearXS::Default_Name()); 89 if ( gammaXS == nullptr ) { 89 if ( gammaXS == nullptr ) { 90 gammaXS = new G4PhotoNuclearCrossSection 90 gammaXS = new G4PhotoNuclearCrossSection; 91 } 91 } 92 } 92 } 93 93 94 // reuse existing pre-compound model 94 // reuse existing pre-compound model 95 G4GeneratorPrecompoundInterface* precoInterf 95 G4GeneratorPrecompoundInterface* precoInterface 96 = new G4GeneratorPrecompoundInterface(); 96 = new G4GeneratorPrecompoundInterface(); 97 G4HadronicInteraction* p = 97 G4HadronicInteraction* p = 98 G4HadronicInteractionRegistry::Instance()- 98 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO"); 99 G4VPreCompoundModel* pre = static_cast<G4VPr 99 G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p); 100 if(!pre) { pre = new G4PreCompoundModel(); } 100 if(!pre) { pre = new G4PreCompoundModel(); } 101 precoInterface->SetDeExcitation(pre); 101 precoInterface->SetDeExcitation(pre); 102 102 103 // string model 103 // string model 104 ftfp = new G4TheoFSGenerator(); 104 ftfp = new G4TheoFSGenerator(); 105 ftfp->SetTransport(precoInterface); 105 ftfp->SetTransport(precoInterface); 106 theFragmentation = new G4LundStringFragmenta 106 theFragmentation = new G4LundStringFragmentation(); 107 theStringDecay = new G4ExcitedStringDecay(th 107 theStringDecay = new G4ExcitedStringDecay(theFragmentation); 108 G4FTFModel* theStringModel = new G4FTFModel( 108 G4FTFModel* theStringModel = new G4FTFModel(); 109 theStringModel->SetFragmentationModel(theStr 109 theStringModel->SetFragmentationModel(theStringDecay); 110 ftfp->SetHighEnergyGenerator(theStringModel) 110 ftfp->SetHighEnergyGenerator(theStringModel); 111 111 112 // Build Bertini model 112 // Build Bertini model 113 bert = new G4CascadeInterface(); 113 bert = new G4CascadeInterface(); 114 114 115 // Creator model ID 115 // Creator model ID 116 secID = G4PhysicsModelCatalog::GetModelID( " 116 secID = G4PhysicsModelCatalog::GetModelID( "model_" + GetModelName() ); 117 } 117 } 118 118 119 G4ElectroVDNuclearModel::~G4ElectroVDNuclearMo 119 G4ElectroVDNuclearModel::~G4ElectroVDNuclearModel() 120 { 120 { 121 delete theFragmentation; 121 delete theFragmentation; 122 delete theStringDecay; 122 delete theStringDecay; 123 } 123 } 124 124 125 void G4ElectroVDNuclearModel::ModelDescription 125 void G4ElectroVDNuclearModel::ModelDescription(std::ostream& outFile) const 126 { 126 { 127 outFile << "G4ElectroVDNuclearModel handles 127 outFile << "G4ElectroVDNuclearModel handles the inelastic scattering\n" 128 << "of e- and e+ from nuclei using t 128 << "of e- and e+ from nuclei using the equivalent photon\n" 129 << "approximation in which the incom 129 << "approximation in which the incoming lepton generates a\n" 130 << "virtual photon at the electromag 130 << "virtual photon at the electromagnetic vertex, and the\n" 131 << "virtual photon is converted to a 131 << "virtual photon is converted to a real photon. At low\n" 132 << "energies, the photon interacts d 132 << "energies, the photon interacts directly with the nucleus\n" 133 << "using the Bertini cascade. At h 133 << "using the Bertini cascade. At high energies the photon\n" 134 << "is converted to a pi0 which inte 134 << "is converted to a pi0 which interacts using the FTFP\n" 135 << "model. The electro- and gamma-n 135 << "model. The electro- and gamma-nuclear cross sections of\n" 136 << "M. Kossov are used to generate t 136 << "M. Kossov are used to generate the virtual photon spectrum\n"; 137 } 137 } 138 138 139 139 140 G4HadFinalState* 140 G4HadFinalState* 141 G4ElectroVDNuclearModel::ApplyYourself(const G 141 G4ElectroVDNuclearModel::ApplyYourself(const G4HadProjectile& aTrack, 142 G4Nucle 142 G4Nucleus& targetNucleus) 143 { 143 { 144 // Set up default particle change (just re 144 // Set up default particle change (just returns initial state) 145 theParticleChange.Clear(); 145 theParticleChange.Clear(); 146 theParticleChange.SetStatusChange(isAlive) 146 theParticleChange.SetStatusChange(isAlive); 147 leptonKE = aTrack.GetKineticEnergy(); 147 leptonKE = aTrack.GetKineticEnergy(); 148 theParticleChange.SetEnergyChange(leptonKE 148 theParticleChange.SetEnergyChange(leptonKE); 149 theParticleChange.SetMomentumChange(aTrack 149 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit() ); 150 150 151 // Set up sanity checks for real photon pr 151 // Set up sanity checks for real photon production 152 G4DynamicParticle lepton(aTrack.GetDefinit 152 G4DynamicParticle lepton(aTrack.GetDefinition(), aTrack.Get4Momentum() ); 153 153 154 // Need to call GetElementCrossSection bef 154 // Need to call GetElementCrossSection before calling GetEquivalentPhotonEnergy. 155 const G4Material* mat = aTrack.GetMaterial 155 const G4Material* mat = aTrack.GetMaterial(); 156 G4int targZ = targetNucleus.GetZ_asInt(); 156 G4int targZ = targetNucleus.GetZ_asInt(); 157 electroXS->GetElementCrossSection(&lepton, 157 electroXS->GetElementCrossSection(&lepton, targZ, mat); 158 158 159 photonEnergy = electroXS->GetEquivalentPho 159 photonEnergy = electroXS->GetEquivalentPhotonEnergy(); 160 // Photon energy cannot exceed lepton ener 160 // Photon energy cannot exceed lepton energy 161 if (photonEnergy < leptonKE) { 161 if (photonEnergy < leptonKE) { 162 photonQ2 = electroXS->GetEquivalentPho 162 photonQ2 = electroXS->GetEquivalentPhotonQ2(photonEnergy); 163 G4double dM = G4Proton::Proton()->GetP 163 G4double dM = G4Proton::Proton()->GetPDGMass() + G4Neutron::Neutron()->GetPDGMass(); 164 // Photon 164 // Photon 165 if (photonEnergy > photonQ2/dM) { 165 if (photonEnergy > photonQ2/dM) { 166 // Produce recoil lepton and trans 166 // Produce recoil lepton and transferred photon 167 G4DynamicParticle* transferredPhot 167 G4DynamicParticle* transferredPhoton = CalculateEMVertex(aTrack, targetNucleus); 168 // Interact gamma with nucleus 168 // Interact gamma with nucleus 169 if (transferredPhoton) CalculateHa 169 if (transferredPhoton) CalculateHadronicVertex(transferredPhoton, targetNucleus); 170 } 170 } 171 } 171 } 172 return &theParticleChange; 172 return &theParticleChange; 173 } 173 } 174 174 175 175 176 G4DynamicParticle* 176 G4DynamicParticle* 177 G4ElectroVDNuclearModel::CalculateEMVertex(con 177 G4ElectroVDNuclearModel::CalculateEMVertex(const G4HadProjectile& aTrack, 178 G4N 178 G4Nucleus& targetNucleus) 179 { 179 { 180 G4DynamicParticle photon(G4Gamma::Gamma(), p 180 G4DynamicParticle photon(G4Gamma::Gamma(), photonEnergy, 181 G4ThreeVector(0.,0. 181 G4ThreeVector(0.,0.,1.) ); 182 182 183 // Get gamma cross section at Q**2 = 0 (real 183 // Get gamma cross section at Q**2 = 0 (real gamma) 184 G4int targZ = targetNucleus.GetZ_asInt(); 184 G4int targZ = targetNucleus.GetZ_asInt(); 185 const G4Material* mat = aTrack.GetMaterial() 185 const G4Material* mat = aTrack.GetMaterial(); 186 G4double sigNu = 186 G4double sigNu = 187 gammaXS->GetElementCrossSection(&photon, t 187 gammaXS->GetElementCrossSection(&photon, targZ, mat); 188 188 189 // Change real gamma energy to equivalent en 189 // Change real gamma energy to equivalent energy and get cross section at that energy 190 G4double dM = G4Proton::Proton()->GetPDGMass 190 G4double dM = G4Proton::Proton()->GetPDGMass() + G4Neutron::Neutron()->GetPDGMass(); 191 photon.SetKineticEnergy(photonEnergy - photo 191 photon.SetKineticEnergy(photonEnergy - photonQ2/dM); 192 G4double sigK = 192 G4double sigK = 193 gammaXS->GetElementCrossSection(&photon, t 193 gammaXS->GetElementCrossSection(&photon, targZ, mat); 194 G4double rndFraction = electroXS->GetVirtual 194 G4double rndFraction = electroXS->GetVirtualFactor(photonEnergy, photonQ2); 195 195 196 // No gamma produced, return null ptr 196 // No gamma produced, return null ptr 197 if (sigNu*G4UniformRand() > sigK*rndFraction 197 if (sigNu*G4UniformRand() > sigK*rndFraction) return 0; 198 198 199 // Scatter the lepton 199 // Scatter the lepton 200 G4double mProj = aTrack.GetDefinition()->Get 200 G4double mProj = aTrack.GetDefinition()->GetPDGMass(); 201 G4double mProj2 = mProj*mProj; 201 G4double mProj2 = mProj*mProj; 202 G4double iniE = leptonKE + mProj; 202 G4double iniE = leptonKE + mProj; // Total energy of incident lepton 203 G4double finE = iniE - photonEnergy; 203 G4double finE = iniE - photonEnergy; // Total energy of scattered lepton 204 theParticleChange.SetEnergyChange(finE-mProj 204 theParticleChange.SetEnergyChange(finE-mProj); 205 G4double iniP = std::sqrt(iniE*iniE-mProj2); 205 G4double iniP = std::sqrt(iniE*iniE-mProj2); // Incident lepton momentum 206 G4double finP = std::sqrt(finE*finE-mProj2); 206 G4double finP = std::sqrt(finE*finE-mProj2); // Scattered lepton momentum 207 G4double cost = (iniE*finE - mProj2 - photon 207 G4double cost = (iniE*finE - mProj2 - photonQ2/2.)/iniP/finP; // cos(theta) from Q**2 208 if (cost > 1.) cost= 1.; 208 if (cost > 1.) cost= 1.; 209 if (cost < -1.) cost=-1.; 209 if (cost < -1.) cost=-1.; 210 G4double sint = std::sqrt(1.-cost*cost); 210 G4double sint = std::sqrt(1.-cost*cost); 211 211 212 G4ThreeVector dir = aTrack.Get4Momentum().ve 212 G4ThreeVector dir = aTrack.Get4Momentum().vect().unit(); 213 G4ThreeVector ortx = dir.orthogonal().unit() 213 G4ThreeVector ortx = dir.orthogonal().unit(); // Ortho-normal to scattering plane 214 G4ThreeVector orty = dir.cross(ortx); 214 G4ThreeVector orty = dir.cross(ortx); // Third unit vector 215 G4double phi = twopi*G4UniformRand(); 215 G4double phi = twopi*G4UniformRand(); 216 G4double sinx = sint*std::sin(phi); 216 G4double sinx = sint*std::sin(phi); 217 G4double siny = sint*std::cos(phi); 217 G4double siny = sint*std::cos(phi); 218 G4ThreeVector findir = cost*dir+sinx*ortx+si 218 G4ThreeVector findir = cost*dir+sinx*ortx+siny*orty; 219 theParticleChange.SetMomentumChange(findir); 219 theParticleChange.SetMomentumChange(findir); // change lepton direction 220 220 221 // Create a gamma with momentum equal to mom 221 // Create a gamma with momentum equal to momentum transfer 222 G4ThreeVector photonMomentum = iniP*dir - fi 222 G4ThreeVector photonMomentum = iniP*dir - finP*findir; 223 G4DynamicParticle* gamma = new G4DynamicPart 223 G4DynamicParticle* gamma = new G4DynamicParticle(G4Gamma::Gamma(), 224 224 photonEnergy, photonMomentum); 225 return gamma; 225 return gamma; 226 } 226 } 227 227 228 228 229 void 229 void 230 G4ElectroVDNuclearModel::CalculateHadronicVert 230 G4ElectroVDNuclearModel::CalculateHadronicVertex(G4DynamicParticle* incident, 231 231 G4Nucleus& target) 232 { 232 { 233 G4HadFinalState* hfs = 0; 233 G4HadFinalState* hfs = 0; 234 G4double gammaE = incident->GetTotalEnergy() 234 G4double gammaE = incident->GetTotalEnergy(); 235 235 236 if (gammaE < 10*GeV) { 236 if (gammaE < 10*GeV) { 237 G4HadProjectile projectile(*incident); 237 G4HadProjectile projectile(*incident); 238 hfs = bert->ApplyYourself(projectile, targ 238 hfs = bert->ApplyYourself(projectile, target); 239 } else { 239 } else { 240 // At high energies convert incident gamma 240 // At high energies convert incident gamma to a pion 241 G4double piMass = G4PionZero::PionZero()-> 241 G4double piMass = G4PionZero::PionZero()->GetPDGMass(); 242 G4double piMom = std::sqrt(gammaE*gammaE - 242 G4double piMom = std::sqrt(gammaE*gammaE - piMass*piMass); 243 G4ThreeVector piMomentum(incident->GetMome 243 G4ThreeVector piMomentum(incident->GetMomentumDirection() ); 244 piMomentum *= piMom; 244 piMomentum *= piMom; 245 G4DynamicParticle theHadron(G4PionZero::Pi 245 G4DynamicParticle theHadron(G4PionZero::PionZero(), piMomentum); 246 G4HadProjectile projectile(theHadron); 246 G4HadProjectile projectile(theHadron); 247 hfs = ftfp->ApplyYourself(projectile, targ 247 hfs = ftfp->ApplyYourself(projectile, target); 248 } 248 } 249 249 250 delete incident; 250 delete incident; 251 251 252 // Assign the creator model ID to the second 252 // Assign the creator model ID to the secondaries 253 for ( size_t i = 0; i < hfs->GetNumberOfSeco 253 for ( size_t i = 0; i < hfs->GetNumberOfSecondaries(); ++i ) { 254 hfs->GetSecondary( i )->SetCreatorModelID( 254 hfs->GetSecondary( i )->SetCreatorModelID( secID ); 255 } 255 } 256 256 257 // Copy secondaries from sub-model to model 257 // Copy secondaries from sub-model to model 258 theParticleChange.AddSecondaries(hfs); 258 theParticleChange.AddSecondaries(hfs); 259 } 259 } 260 260 261 261