ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MuonVDNuclearModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MuonVDNuclearModel.cc
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 //
27 // Author: D.H. Wright
28 // Date: 2 February 2011
29 //
30 // Description: model of muon nuclear interaction in which a gamma from
31 // the virtual photon spectrum interacts in the nucleus as
32 // a real gamma at low energies and as a pi0 at high energies.
33 // Kokoulin's muon cross section and equivalent gamma spectrum
34 // are used.
35 //
36 
37 #include "G4MuonVDNuclearModel.hh"
38 
39 #include "Randomize.hh"
40 #include "G4Log.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4CascadeInterface.hh"
44 #include "G4TheoFSGenerator.hh"
45 #include "G4GeneratorPrecompoundInterface.hh"
46 #include "G4PreCompoundModel.hh"
47 #include "G4LundStringFragmentation.hh"
48 #include "G4ExcitedStringDecay.hh"
49 #include "G4FTFModel.hh"
50 #include "G4HadronicInteractionRegistry.hh"
51 #include "G4KokoulinMuonNuclearXS.hh"
52 #include "G4CrossSectionDataSetRegistry.hh"
53 #include "G4ElementData.hh"
54 #include "G4Physics2DVector.hh"
55 #include "G4Pow.hh"
56 
57 const G4int G4MuonVDNuclearModel::zdat[] = {1, 4, 13, 29, 92};
58 const G4double G4MuonVDNuclearModel::adat[] = {1.01,9.01,26.98,63.55,238.03};
59 const G4double G4MuonVDNuclearModel::tdat[] = {
60  1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3,
61  1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4,
62  1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5,
63  1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6,
64  1.e7,2.e7,3.e7,4.e7,5.e7,6.e7,7.e7,8.e7,9.e7,
65  1.e8,2.e8,3.e8,4.e8,5.e8,6.e8,7.e8,8.e8,9.e8,
66  1.e9,2.e9,3.e9,4.e9,5.e9,6.e9,7.e9,8.e9,9.e9,
67  1.e10,2.e10,3.e10,4.e10,5.e10,6.e10,7.e10,8.e10,9.e10,1.e11};
68 
69 G4ElementData* G4MuonVDNuclearModel::fElementData = nullptr;
70 
71 G4MuonVDNuclearModel::G4MuonVDNuclearModel()
72  : G4HadronicInteraction("G4MuonVDNuclearModel"),isMaster(false)
73 {
74  muNucXS = (G4KokoulinMuonNuclearXS*)G4CrossSectionDataSetRegistry::Instance()->
75  GetCrossSectionDataSet(G4KokoulinMuonNuclearXS::Default_Name());
76 
77  SetMinEnergy(0.0);
78  SetMaxEnergy(1*CLHEP::PeV);
79  CutFixed = 0.2*CLHEP::GeV;
80 
81  if(!fElementData && G4Threading::IsMasterThread()) {
82  fElementData = new G4ElementData();
83  MakeSamplingTable();
84  isMaster = true;
85  }
86 
87  // reuse existing pre-compound model
88  G4GeneratorPrecompoundInterface* precoInterface
89  = new G4GeneratorPrecompoundInterface();
90  G4HadronicInteraction* p =
91  G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
92  G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
93  if(!pre) { pre = new G4PreCompoundModel(); }
94  precoInterface->SetDeExcitation(pre);
95 
96  // Build FTFP model
97  ftfp = new G4TheoFSGenerator();
98  ftfp->SetTransport(precoInterface);
99  theFragmentation = new G4LundStringFragmentation();
100  theStringDecay = new G4ExcitedStringDecay(theFragmentation);
101  G4FTFModel* theStringModel = new G4FTFModel;
102  theStringModel->SetFragmentationModel(theStringDecay);
103  ftfp->SetHighEnergyGenerator(theStringModel);
104 
105  // Build Bertini cascade
106  bert = new G4CascadeInterface();
107 }
108 
109 G4MuonVDNuclearModel::~G4MuonVDNuclearModel()
110 {
111  delete theFragmentation;
112  delete theStringDecay;
113 
114  if(isMaster) {
115  delete fElementData;
116  fElementData = nullptr;
117  }
118 }
119 
120 G4HadFinalState*
121 G4MuonVDNuclearModel::ApplyYourself(const G4HadProjectile& aTrack,
122  G4Nucleus& targetNucleus)
123 {
124  theParticleChange.Clear();
125 
126  // For very low energy, return initial track
127  G4double epmax = aTrack.GetTotalEnergy() - 0.5*proton_mass_c2;
128  if (epmax <= CutFixed) {
129  theParticleChange.SetStatusChange(isAlive);
130  theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
131  theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
132  return &theParticleChange;
133  }
134 
135  // Produce recoil muon and transferred photon
136  G4DynamicParticle* transferredPhoton = CalculateEMVertex(aTrack, targetNucleus);
137 
138  // Interact the gamma with the nucleus
139  CalculateHadronicVertex(transferredPhoton, targetNucleus);
140  return &theParticleChange;
141 }
142 
143 G4DynamicParticle*
144 G4MuonVDNuclearModel::CalculateEMVertex(const G4HadProjectile& aTrack,
145  G4Nucleus& targetNucleus)
146 {
147  // Select sampling table
148  G4double KineticEnergy = aTrack.GetKineticEnergy();
149  G4double TotalEnergy = aTrack.GetTotalEnergy();
150  G4double Mass = G4MuonMinus::MuonMinus()->GetPDGMass();
151  G4Pow* g4calc = G4Pow::GetInstance();
152  G4double lnZ = g4calc->logZ(targetNucleus.GetZ_asInt());
153 
154  G4double epmin = CutFixed;
155  G4double epmax = TotalEnergy - 0.5*proton_mass_c2;
156  G4double m0 = CutFixed;
157 
158  G4double delmin = 1.e10;
159  G4double del;
160  G4int izz = 0;
161  G4int itt = 0;
162 
163  for (G4int iz = 0; iz < nzdat; ++iz) {
164  del = std::abs(lnZ - g4calc->logZ(zdat[iz]));
165  if (del < delmin) {
166  delmin = del;
167  izz = iz;
168  }
169  }
170 
171  delmin = 1.e10;
172  for (G4int it = 0; it < ntdat; ++it) {
173  del = std::abs(G4Log(KineticEnergy)-G4Log(tdat[it]) );
174  if (del < delmin) {
175  delmin = del;
176  itt = it;
177  }
178  }
179 
180  // Sample the energy transfer according to the probability table
181  G4double r = G4UniformRand();
182 
183  G4int iy;
184 
185  G4int Z = zdat[izz];
186 
187  for(iy = 0; iy<NBIN; ++iy) {
188 
189  G4double pvv = fElementData->GetElement2DData(Z)->GetValue(iy, itt);
190  if(pvv >= r) { break; }
191  }
192 
193  // Sampling is done uniformly in y in the bin
194  G4double pvx = fElementData->GetElement2DData(Z)->GetX(iy);
195  G4double pvx1 = fElementData->GetElement2DData(Z)->GetX(iy+1);
196  G4double y = pvx + G4UniformRand() * (pvx1 - pvx);
197 
198  G4double x = G4Exp(y);
199  G4double ep = epmin*G4Exp(x*G4Log(epmax/epmin) );
200 
201  // Sample scattering angle of mu, but first t should be sampled.
202  G4double yy = ep/TotalEnergy;
203  G4double tmin = Mass*Mass*yy*yy/(1.-yy);
204  G4double tmax = 2.*proton_mass_c2*ep;
205  G4double t1;
206  G4double t2;
207  if (m0 < ep) {
208  t1 = m0*m0;
209  t2 = ep*ep;
210  } else {
211  t1 = ep*ep;
212  t2 = m0*m0;
213  }
214 
215  G4double w1 = tmax*t1;
216  G4double w2 = tmax+t1;
217  G4double w3 = tmax*(tmin+t1)/(tmin*w2);
218  G4double y1 = 1.-yy;
219  G4double y2 = 0.5*yy*yy;
220  G4double y3 = y1+y2;
221 
222  G4double t;
223  G4double rej;
224 
225  // Now sample t
226  G4int ntry = 0;
227  do
228  {
229  ntry += 1;
230  if (ntry > 10000) {
231  G4ExceptionDescription eda;
232  eda << " While count exceeded " << G4endl;
233  G4Exception("G4MuonVDNuclearModel::CalculateEMVertex()", "HAD_RPG_100", JustWarning, eda);
234  break;
235  }
236 
237  t = w1/(w2*G4Exp(G4UniformRand()*G4Log(w3))-tmax);
238  rej = (1.-t/tmax)*(y1*(1.-tmin/t)+y2)/(y3*(1.-t/t2));
239  } while (G4UniformRand() > rej) ; /* Loop checking, 01.09.2015, D.Wright */
240 
241  // compute angle from t
242  G4double sinth2 =
243  0.5*(t-tmin)/(2.*(TotalEnergy*(TotalEnergy-ep)-Mass*Mass)-tmin);
244  G4double theta = std::acos(1. - 2.*sinth2);
245 
246  G4double phi = twopi*G4UniformRand();
247  G4double sinth = std::sin(theta);
248  G4double dirx = sinth*std::cos(phi);
249  G4double diry = sinth*std::sin(phi);
250  G4double dirz = std::cos(theta);
251  G4ThreeVector finalDirection(dirx,diry,dirz);
252  G4ThreeVector ParticleDirection(aTrack.Get4Momentum().vect().unit() );
253  finalDirection.rotateUz(ParticleDirection);
254 
255  G4double NewKinEnergy = KineticEnergy - ep;
256  G4double finalMomentum = std::sqrt(NewKinEnergy*(NewKinEnergy+2.*Mass) );
257  G4double Ef = NewKinEnergy + Mass;
258  G4double initMomentum = std::sqrt(KineticEnergy*(TotalEnergy+Mass) );
259 
260  // Set energy and direction of scattered primary in theParticleChange
261  theParticleChange.SetStatusChange(isAlive);
262  theParticleChange.SetEnergyChange(NewKinEnergy);
263  theParticleChange.SetMomentumChange(finalDirection);
264 
265  // Now create the emitted gamma
266  G4LorentzVector primaryMomentum(initMomentum*ParticleDirection, TotalEnergy);
267  G4LorentzVector fsMomentum(finalMomentum*finalDirection, Ef);
268  G4LorentzVector momentumTransfer = primaryMomentum - fsMomentum;
269 
270  G4DynamicParticle* gamma =
271  new G4DynamicParticle(G4Gamma::Gamma(), momentumTransfer);
272 
273  return gamma;
274 }
275 
276 void
277 G4MuonVDNuclearModel::CalculateHadronicVertex(G4DynamicParticle* incident,
278  G4Nucleus& target)
279 {
280  G4HadFinalState* hfs = 0;
281  G4double gammaE = incident->GetTotalEnergy();
282 
283  if (gammaE < 10*GeV) {
284  G4HadProjectile projectile(*incident);
285  hfs = bert->ApplyYourself(projectile, target);
286  } else {
287  // convert incident gamma to a pi0
288  G4double piMass = G4PionZero::PionZero()->GetPDGMass();
289  G4double piKE = incident->GetTotalEnergy() - piMass;
290  G4double piMom = std::sqrt(piKE*(piKE + 2*piMass) );
291  G4ThreeVector piMomentum(incident->GetMomentumDirection() );
292  piMomentum *= piMom;
293  G4DynamicParticle theHadron(G4PionZero::PionZero(), piMomentum);
294  G4HadProjectile projectile(theHadron);
295  hfs = ftfp->ApplyYourself(projectile, target);
296  }
297 
298  delete incident;
299 
300  // Copy secondaries from sub-model to model
301  theParticleChange.AddSecondaries(hfs);
302 }
303 
304 
305 void G4MuonVDNuclearModel::MakeSamplingTable()
306 {
307  G4int nbin;
308  G4double KineticEnergy;
309  G4double TotalEnergy;
310  G4double Maxep;
311  G4double CrossSection;
312 
313  G4double c;
314  G4double y;
315  G4double ymin,ymax;
316  G4double dy,yy;
317  G4double dx,x;
318  G4double ep;
319 
320  G4double AtomicNumber;
321  G4double AtomicWeight;
322 
323  G4double mumass = G4MuonMinus::MuonMinus()->GetPDGMass();
324 
325  for (G4int iz = 0; iz < nzdat; ++iz) {
326  AtomicNumber = zdat[iz];
327  AtomicWeight = adat[iz]*(g/mole);
328 
329  G4Physics2DVector* pv = new G4Physics2DVector(NBIN+1,ntdat+1);
330  G4double pvv;
331 
332  for (G4int it = 0; it < ntdat; ++it) {
333  KineticEnergy = tdat[it];
334  TotalEnergy = KineticEnergy + mumass;
335  Maxep = TotalEnergy - 0.5*proton_mass_c2;
336 
337  CrossSection = 0.0;
338 
339  // Calculate the differential cross section
340  // numerical integration in log .........
341  c = G4Log(Maxep/CutFixed);
342  ymin = -5.0;
343  ymax = 0.0;
344  dy = (ymax-ymin)/NBIN;
345 
346  nbin=-1;
347 
348  y = ymin - 0.5*dy;
349  yy = ymin - dy;
350  for (G4int i = 0; i < NBIN; ++i) {
351  y += dy;
352  x = G4Exp(y);
353  yy += dy;
354  dx = G4Exp(yy+dy)-G4Exp(yy);
355 
356  ep = CutFixed*G4Exp(c*x);
357 
358  CrossSection +=
359  ep*dx*muNucXS->ComputeDDMicroscopicCrossSection(KineticEnergy,
360  AtomicNumber,
361  AtomicWeight, ep);
362  if (nbin < NBIN) {
363  ++nbin;
364  pv->PutValue(nbin, it, CrossSection);
365  pv->PutX(nbin, y);
366  }
367  }
368  pv->PutX(NBIN, 0.);
369 
370  if (CrossSection > 0.0) {
371  for (G4int ib = 0; ib <= nbin; ++ib) {
372  pvv = pv->GetValue(ib, it);
373  pvv = pvv/CrossSection;
374  pv->PutValue(ib, it, pvv);
375  }
376  }
377  } // loop on it
378 
379  fElementData->InitialiseForElement(zdat[iz], pv);
380  } // loop on iz
381 
382  // G4cout << " Kokoulin XS = "
383  // << muNucXS->ComputeDDMicroscopicCrossSection(1*GeV, 20.0,
384  // 40.0*g/mole, 0.3*GeV)/millibarn
385  // << G4endl;
386 }
387 
388 void G4MuonVDNuclearModel::ModelDescription(std::ostream& outFile) const
389 {
390  outFile << "G4MuonVDNuclearModel handles the inelastic scattering\n"
391  << "of mu- and mu+ from nuclei using the equivalent photon\n"
392  << "approximation in which the incoming lepton generates a\n"
393  << "virtual photon at the electromagnetic vertex, and the\n"
394  << "virtual photon is converted to a real photon. At low\n"
395  << "energies, the photon interacts directly with the nucleus\n"
396  << "using the Bertini cascade. At high energies the photon\n"
397  << "is converted to a pi0 which interacts using the FTFP\n"
398  << "model. The muon-nuclear cross sections of R. Kokoulin \n"
399  << "are used to generate the virtual photon spectrum\n";
400 }
401