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