source: trunk/source/processes/hadronic/models/photolepton_hadron/muon_nuclear/src/G4MuonNucleusInteractionModel.cc @ 1340

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

update ti head
File size: 13.4 KB
Line 
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// $Id: G4MuonNucleusInteractionModel.cc,v 1.6 2006/06/29 20:57:36 gunter Exp $
28// GEANT4 tag $Name: geant4-09-03-ref-09 $
29//
30// --------------------------------------------------------------
31// G4MuonNucleusInteractionModel.cc
32//
33//     M.Takahata (Makoto.Takahata@cern.ch)
34
35#include "G4MuonNucleusInteractionModel.hh"
36
37
38//-----------------------------------------------------------------------------
39  G4MuonNucleusInteractionModel::G4MuonNucleusInteractionModel()
40    : G4LeptonHadronInteractionModel()
41//-----------------------------------------------------------------------------
42  {
43    // build the physics vector
44    Nbin = 90;
45    kEmin = 1.0e-5*GeV;
46    kEmax = 1.0e+4*GeV;
47    cascadeModelMarginalEnergy = 25.0*GeV;
48    theCoefficientVector = new G4PhysicsLogVector(kEmin, kEmax, Nbin);
49    makePhysicsVector();
50
51    // construct variables
52    LEPionMinusInelastic  = new G4LEPionMinusInelastic;
53    LEPionPlusInelastic   = new G4LEPionPlusInelastic;
54    HEPionMinusInelastic  = new G4HEPionMinusInelastic;
55    HEPionPlusInelastic   = new G4HEPionPlusInelastic;
56  }
57
58
59//-----------------------------------------------------------------------------
60  G4MuonNucleusInteractionModel::~G4MuonNucleusInteractionModel()
61//-----------------------------------------------------------------------------
62  {
63    delete LEPionMinusInelastic;
64    delete LEPionPlusInelastic;
65    delete HEPionMinusInelastic;
66    delete HEPionPlusInelastic;
67
68    delete theCoefficientVector;
69  }
70
71
72//-----------------------------------------------------------------------------
73  G4double G4MuonNucleusInteractionModel::tetal[35] = {
74//-----------------------------------------------------------------------------
75    1.0000000,  0.9999995,  0.9999990,  0.9999981,  0.9999962,
76    0.9999943,  0.9999905,  0.9999847,  0.9999752,  0.9999599,
77    0.9999352,  0.9998951,  0.9998302,  0.9997253,  0.9995556,
78    0.9992810,  0.9988368,  0.9981183,  0.9969561,  0.9950773,
79    0.9920409,  0.9871377,  0.9792297,  0.9665010,  0.9460785,
80    0.9134827,  0.8618938,  0.7813507,  0.6583430,  0.4770452,
81    0.2247237, -0.0955139, -0.4461272, -0.7495149, -0.9900000
82  };
83
84
85//-----------------------------------------------------------------------------
86  G4double G4MuonNucleusInteractionModel::xeml[23] = {
87//-----------------------------------------------------------------------------
88    1.000,  0.998,  0.997,  0.996,  0.995,  0.994, 0.992, 0.990,
89    0.970,  0.950,  0.920,  0.890,  0.850,  0.800, 0.750, 0.700,
90    0.600,  0.500,  0.400,  0.300,  0.200,  0.100, 0.050
91  };
92
93
94//-----------------------------------------------------------------------------
95  G4double G4MuonNucleusInteractionModel::computeMicroscopicCrossSection
96    (const G4Track &muonTrack)
97//-----------------------------------------------------------------------------
98  {
99    const G4DynamicParticle *muonDynamics = muonTrack.GetDynamicParticle();
100    G4double kineticEnergy = muonDynamics->GetKineticEnergy();
101    G4double muonMass      = muonDynamics->GetDefinition()->GetPDGMass();
102
103    G4double totalEnergy = kineticEnergy + muonMass;
104
105    G4double microscopicCrossSection;
106    if(totalEnergy <= 30.0*GeV) {
107      microscopicCrossSection
108        = 0.0003*millibarn;
109    } else {
110      microscopicCrossSection
111        = 0.0003*std::pow((totalEnergy/(30.0*GeV)), 0.25)*millibarn;
112    }
113
114    return microscopicCrossSection;
115  }
116
117
118//-----------------------------------------------------------------------------
119  void G4MuonNucleusInteractionModel::makePhysicsVector()
120//-----------------------------------------------------------------------------
121  {
122    G4double Ei, Ef;  // initial and final energy of incident muon;
123    G4double muonMass = G4MuonMinus::MuonMinus()->GetPDGMass();
124
125    for (G4int i=0; i<=(Nbin-1); i++)
126    {
127      G4double totalCrossSection = 0.0;
128      Ei = theCoefficientVector->GetLowEdgeEnergy(i) + muonMass;
129      for (G4int j=1; j<=34; j++)
130      {
131        cosTheta = 0.5 * (tetal[j] + tetal[j-1]);
132        for (G4int k=1; k<=22; k++)
133        {
134          Ef = 0.5 * Ei * (xeml[k]+xeml[k-1]);
135          G4double dsigma = computeDifferentialCrossSection(Ei,Ef,cosTheta);
136          totalCrossSection = totalCrossSection
137            + Ei * (tetal[j-1]-tetal[j])*(xeml[k-1]-xeml[k]) * dsigma;
138        }
139      }
140      theCoefficientVector->PutValue(i, totalCrossSection);
141    }
142  }
143
144
145//-----------------------------------------------------------------------------
146  G4VParticleChange* G4MuonNucleusInteractionModel::applyInteractionModel
147    (const G4Track &muonTrack, G4Nucleus &targetNucleus )
148//-----------------------------------------------------------------------------
149  {
150    G4int icos=0, ie1=0;
151    G4double E1=0., P1=0.;
152    G4double rndm[3];
153    G4bool isOutRange;
154
155    // Initialization
156    aParticleChange.Initialize(muonTrack);
157
158    const G4DynamicParticle *muonDynamics = muonTrack.GetDynamicParticle();
159    G4double kineticEnergy = muonDynamics->GetKineticEnergy();
160    G4double totalMomentum = muonDynamics->GetTotalMomentum();
161    G4double totalEnergy   = muonDynamics->GetTotalEnergy();
162    G4double muonMass      = muonDynamics->GetDefinition()->GetPDGMass();
163
164
165    G4double W2 = 0.0;  G4int W2try = 0;
166    while (W2 <= 0.0)
167    {
168      G4double totalCrossSection = 0.0;
169      G4bool interpolated = false;
170      G4double fRndm = G4UniformRand();
171      G4double Hmax
172        = theCoefficientVector->GetValue(kineticEnergy, isOutRange);
173      for (G4int i=1; i<=34; i++)
174      {
175        cosTheta = 0.5 * (tetal[i] + tetal[i-1]);
176        for (G4int j=1; j<=22; j++)
177        {
178          E1 = 0.5 * totalEnergy * (xeml[j]+xeml[j-1]);
179          G4double dsigma
180            = computeDifferentialCrossSection(totalEnergy, E1, cosTheta);
181          totalCrossSection = totalCrossSection
182             + totalEnergy*(tetal[i-1]-tetal[i])*(xeml[j-1]-xeml[j])*dsigma; 
183
184          if((fRndm*Hmax)<totalCrossSection) {
185            interpolated = true;
186            icos = i;  ie1 = j;
187            break;
188          }
189        }
190
191        if(interpolated) {
192          // calculate energy, momentum and angle of outgoing muon
193          CLHEP::RandFlat::shootArray(3, rndm);
194          G4double theta = std::acos(tetal[icos-1]) 
195                          + rndm[0]*(std::acos(tetal[icos])-std::acos(tetal[icos-1]));
196          cosTheta = std::cos(theta);
197          E1 = (xeml[ie1] + rndm[1]*(xeml[ie1-1]-xeml[ie1])) * totalEnergy;
198            if(E1<muonMass) E1 = muonMass + 0.0001*GeV;
199          P1 = std::sqrt(std::abs(E1*E1-muonMass*muonMass));
200
201          // invariant mass of final hadron state must be greater than zero
202          W2 = proton_mass_c2*proton_mass_c2
203              +2.0*proton_mass_c2*(totalEnergy-E1)
204              -2.0*(totalEnergy*E1-totalMomentum*P1*cosTheta-muonMass*muonMass);
205
206          break;
207        }
208      }
209
210      W2try++;
211      if (W2try>100) return &aParticleChange;
212    }
213
214    // calculate momentum of outgoing muon / pion(photon)
215    G4double sinTheta = std::sqrt(std::abs(1.0 - cosTheta*cosTheta));
216    G4double phi = rndm[2]*twopi;
217    G4ThreeVector muonDirection(sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
218    G4ThreeVector muonDirectionInit = muonTrack.GetMomentumDirection();
219    muonDirection.rotateUz(muonDirectionInit);
220
221    G4ParticleMomentum pionMomentum
222      = muonDynamics->GetMomentum() - P1*muonDirection;
223    G4double muonKineticEnergy
224      = std::sqrt(P1*P1 + muonMass*muonMass) - muonMass;
225    aParticleChange.ProposeMomentumDirection(muonDirection);
226    aParticleChange.ProposeEnergy(muonKineticEnergy);
227    aParticleChange.ProposeTrackStatus(fAlive);
228
229
230    // virtual photon is exchanged with a pion of same Q2
231    // select pi+/pi- randomly and generate pion track
232    G4ParticleDefinition* pdPion;
233    if(CLHEP::RandBit::shootBit())
234      pdPion = G4PionMinus::PionMinusDefinition();
235    else
236      pdPion = G4PionPlus::PionPlusDefinition();
237
238    G4DynamicParticle* pionDynamics
239        = new G4DynamicParticle(pdPion, pionMomentum);
240
241    G4Track* pionTrack = new G4Track(pionDynamics,
242                                     muonTrack.GetGlobalTime(), 
243                                     muonTrack.GetPosition() );
244    pionTrack->SetStep(muonTrack.GetStep());
245
246    // Invoke pion-nucleus inelastic process
247    invokePionNucleus(*pionTrack, targetNucleus);
248
249    // Termination
250    delete pionTrack;
251
252    return &aParticleChange;
253  }
254
255
256//-----------------------------------------------------------------------------
257  void G4MuonNucleusInteractionModel::invokePionNucleus
258    (const G4Track &pionTrack, G4Nucleus &targetNucleus )
259//-----------------------------------------------------------------------------
260  {
261    // force interaction of pion with target nucleus
262    G4double pionKineticEnergy = pionTrack.GetKineticEnergy();
263    if(pionTrack.GetDefinition()->GetParticleName() == "pi-") {
264      if(pionKineticEnergy <= cascadeModelMarginalEnergy)
265        pionChange
266          = LEPionMinusInelastic->ApplyYourself(pionTrack, targetNucleus);
267      else
268        pionChange
269          = HEPionMinusInelastic->ApplyYourself(pionTrack, targetNucleus);
270    } else if(pionTrack.GetDefinition()->GetParticleName() == "pi+") {
271      if(pionKineticEnergy <= cascadeModelMarginalEnergy)
272        pionChange
273          = LEPionPlusInelastic->ApplyYourself(pionTrack, targetNucleus);
274      else
275        pionChange
276          = HEPionPlusInelastic->ApplyYourself(pionTrack, targetNucleus);
277    }
278
279
280    // add local energy deposit
281    G4double localEnergyDeposited = 0.0;
282    localEnergyDeposited = pionChange->GetLocalEnergyDeposit();
283    aParticleChange.ProposeLocalEnergyDeposit(localEnergyDeposited);
284
285
286    // register secondary particles
287    G4int numSecondaries = pionChange->GetNumberOfSecondaries();
288    aParticleChange.SetNumberOfSecondaries(numSecondaries);
289
290    G4ParticleMomentum secondaryMomentum = G4ThreeVector(0.,0.,0.);
291    for(G4int iS=0; iS<=(numSecondaries-1); iS++) {
292      secondaryMomentum
293        = secondaryMomentum + pionChange->GetSecondary(iS)->GetParticle()->GetMomentum();
294      aParticleChange.AddSecondary(pionChange->GetSecondary(iS)->GetParticle());
295    }
296    pionChange->Clear();
297
298    return;
299  }
300
301
302//-----------------------------------------------------------------------------
303  G4double G4MuonNucleusInteractionModel::computeDifferentialCrossSection
304    (G4double initialEnergy, G4double finalEnergy, G4double aCosTheta)
305//-----------------------------------------------------------------------------
306  {
307    G4double muonMass = G4MuonMinus::MuonMinus()->GetPDGMass();
308
309    if(finalEnergy < muonMass) return(0.0);
310    if(aCosTheta >= 1.0) return DBL_MAX;
311
312    G4double initialMomentum
313      = std::sqrt(initialEnergy*initialEnergy - muonMass*muonMass);
314    G4double finalMomentum
315      = std::sqrt(finalEnergy*finalEnergy - muonMass*muonMass);
316
317
318    // calculate momentum transfer (Q2)
319    // and invariant mass of final state of hadrons (W2)
320    G4double Q2
321      = 2.0*(initialEnergy*finalEnergy
322             -initialMomentum*finalMomentum*aCosTheta-muonMass*muonMass);
323    if(Q2 < 0.0) return(0.0);
324
325    G4double W2
326      = proton_mass_c2*proton_mass_c2
327       +2.0*proton_mass_c2*(initialEnergy-finalEnergy)-Q2;
328    if(W2 < 0.0) return(0.0);
329
330
331    // calculate factors
332    //   Nu      : energy transfer
333    //   K       : incident flux of photon
334    //   Epsilon : virtual photon polarization
335    G4double fNu = initialEnergy-finalEnergy;
336    G4double fK = fNu+Q2/(2.0*fNu);
337    G4double fEpsilon
338      = 1.0/(1.0+2.0*((1.0-aCosTheta)/(1.0+aCosTheta))*(Q2+fNu*fNu)/Q2);
339    if(fEpsilon > 1.0) return DBL_MAX;
340
341
342    // calculate photoabsorption cross sections
343    //   fGamma  : flux of transverse photons
344    //   sigma_t : for transverse photons
345    //   sigma_l : for longitudinal photons
346    G4double fGamma
347      = fine_structure_const*fK*finalEnergy
348        / (pi*Q2*initialEnergy*(1.0 - fEpsilon));
349    G4double sigma_t = 0.12*millibarn;
350    G4double sigma_l = 0.3*(1.0-Q2/(1.868*GeV*fNu))*sigma_t;
351    if(sigma_l < 0.) sigma_l = 0.;
352
353    return fGamma*(sigma_t+fEpsilon*sigma_l);
354  }
Note: See TracBrowser for help on using the repository browser.