Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lepto_nuclear/src/G4MuonVDNuclearModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/lepto_nuclear/src/G4MuonVDNuclearModel.cc (Version 11.3.0) and /processes/hadronic/models/lepto_nuclear/src/G4MuonVDNuclearModel.cc (Version 10.0)


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