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