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