G4MuonVDNuclearModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: $
00027 //
00028 // Author:      D.H. Wright
00029 // Date:        2 February 2011
00030 //
00031 // Description: model of muon nuclear interaction in which a gamma from
00032 //              the virtual photon spectrum interacts in the nucleus as
00033 //              a real gamma at low energies and as a pi0 at high energies.
00034 //              Kokoulin's muon cross section and equivalent gamma spectrum
00035 //              are used.
00036 //
00037 
00038 #include "G4MuonVDNuclearModel.hh"
00039 
00040 #include "Randomize.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4CascadeInterface.hh"
00044 #include "G4TheoFSGenerator.hh" 
00045 #include "G4GeneratorPrecompoundInterface.hh"
00046 #include "G4ExcitationHandler.hh"
00047 #include "G4PreCompoundModel.hh"
00048 #include "G4LundStringFragmentation.hh"
00049 #include "G4ExcitedStringDecay.hh"
00050 #include "G4FTFModel.hh"
00051 
00052 G4MuonVDNuclearModel::G4MuonVDNuclearModel()
00053  : G4HadronicInteraction("G4MuonVDNuclearModel")
00054 {
00055   SetMinEnergy(0.0);
00056   SetMaxEnergy(1*PeV);
00057   CutFixed = 0.2*GeV;
00058   NBIN = 1000;
00059 
00060   for (G4int k = 0; k < 5; k++) {
00061     for (G4int j = 0; j < 8; j++) {
00062       for (G4int i = 0; i < 1001; i++) {
00063         proba[k][j][i] = 0.0;
00064         ya[i] = 0.0;
00065       }
00066     }
00067   }
00068 
00069   MakeSamplingTable();
00070 
00071   // Build FTFP model
00072   ftfp = new G4TheoFSGenerator();
00073   precoInterface = new G4GeneratorPrecompoundInterface();
00074   theHandler = new G4ExcitationHandler();
00075   preEquilib = new G4PreCompoundModel(theHandler);
00076   precoInterface->SetDeExcitation(preEquilib);
00077   ftfp->SetTransport(precoInterface);
00078   theFragmentation = new G4LundStringFragmentation();
00079   theStringDecay = new G4ExcitedStringDecay(theFragmentation);    
00080   theStringModel = new G4FTFModel;
00081   theStringModel->SetFragmentationModel(theStringDecay);
00082   ftfp->SetHighEnergyGenerator(theStringModel);
00083 
00084   // Build Bertini cascade
00085   bert = new G4CascadeInterface();
00086 }
00087 
00088 
00089 G4MuonVDNuclearModel::~G4MuonVDNuclearModel()
00090 {
00091   delete ftfp;
00092   delete preEquilib;
00093   delete theFragmentation;
00094   delete theStringDecay;
00095   delete theStringModel;
00096   delete bert;
00097 }
00098 
00099   
00100 G4HadFinalState*
00101 G4MuonVDNuclearModel::ApplyYourself(const G4HadProjectile& aTrack,
00102                                     G4Nucleus& targetNucleus)
00103 {
00104   theParticleChange.Clear();
00105 
00106   // For very low energy, return initial track
00107   G4double epmax = aTrack.GetTotalEnergy() - 0.5*proton_mass_c2;
00108   if (epmax <= CutFixed) {
00109     theParticleChange.SetStatusChange(isAlive);
00110     theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00111     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00112     return &theParticleChange;
00113   }
00114 
00115   // Produce recoil muon and transferred photon
00116   G4DynamicParticle* transferredPhoton = CalculateEMVertex(aTrack, targetNucleus);
00117 
00118   // Interact the gamma with the nucleus
00119   CalculateHadronicVertex(transferredPhoton, targetNucleus);
00120   return &theParticleChange;
00121 }
00122 
00123 
00124 G4DynamicParticle*
00125 G4MuonVDNuclearModel::CalculateEMVertex(const G4HadProjectile& aTrack,
00126                                         G4Nucleus& targetNucleus)
00127 {
00128   // Select sampling table
00129   G4double KineticEnergy = aTrack.GetKineticEnergy();
00130   G4double TotalEnergy = aTrack.GetTotalEnergy();
00131   G4double Mass = G4MuonMinus::MuonMinus()->GetPDGMass();
00132   G4double lnZ = std::log(G4double(targetNucleus.GetZ_asInt() ) );
00133 
00134   G4double epmin = CutFixed;
00135   G4double epmax = TotalEnergy - 0.5*proton_mass_c2;
00136   G4double m0 = 0.2*GeV;
00137 
00138   G4double delmin = 1.e10;
00139   G4double del;
00140   G4int izz = 0;
00141   G4int itt = 0;
00142   G4int NBINminus1 = NBIN - 1;
00143 
00144   G4int nzdat = 5;
00145   G4double zdat[] = {1.,4.,13.,29.,92.};
00146   for (G4int iz = 0; iz < nzdat; iz++) {
00147     del = std::abs(lnZ-std::log(zdat[iz]));
00148     if (del < delmin) {
00149       delmin = del;
00150       izz = iz;
00151     }
00152   }
00153  
00154   G4int ntdat = 8;
00155   G4double tdat[] = {1.e3,1.e4,1.e5,1.e6,1.e7,1.e8,1.e9,1.e10};
00156   delmin = 1.e10;
00157   for (G4int it = 0; it < ntdat; it++) {
00158     del = std::abs(std::log(KineticEnergy)-std::log(tdat[it]) );
00159     if (del < delmin) {
00160       delmin = del;
00161       itt = it;
00162     }
00163   }
00164 
00165   // Sample the energy transfer according to the probability table
00166   G4double r = G4UniformRand();
00167 
00168   G4int iy = -1;
00169   do {
00170        iy += 1 ;
00171      } while (((proba[izz][itt][iy]) < r)&&(iy < NBINminus1)) ;
00172 
00173   // Sampling is done uniformly in y in the bin
00174 
00175   G4double y; 
00176   if (iy < NBIN)
00177     y = ya[iy] + G4UniformRand() * (ya[iy+1] - ya[iy]);
00178   else
00179     y = ya[iy];
00180 
00181   G4double x = std::exp(y);
00182   G4double ep = epmin*std::exp(x*std::log(epmax/epmin) );
00183 
00184   // Sample scattering angle of mu, but first t should be sampled.
00185   G4double yy = ep/TotalEnergy;
00186   G4double tmin = Mass*Mass*yy*yy/(1.-yy);
00187   G4double tmax = 2.*proton_mass_c2*ep;
00188   G4double t1;
00189   G4double t2;
00190   if (m0 < ep) {
00191     t1 = m0*m0;
00192     t2 = ep*ep;
00193   } else {
00194     t1 = ep*ep;
00195     t2 = m0*m0;
00196   }
00197 
00198   G4double w1 = tmax*t1;
00199   G4double w2 = tmax+t1;
00200   G4double w3 = tmax*(tmin+t1)/(tmin*w2);
00201   G4double y1 = 1.-yy;
00202   G4double y2 = 0.5*yy*yy;
00203   G4double y3 = y1+y2;
00204 
00205   G4double t;
00206   G4double rej;
00207 
00208   // Now sample t
00209   G4int ntry = 0;
00210   do
00211   {
00212     ntry += 1;
00213     t = w1/(w2*std::exp(G4UniformRand()*std::log(w3))-tmax);
00214     rej = (1.-t/tmax)*(y1*(1.-tmin/t)+y2)/(y3*(1.-t/t2)); 
00215   } while (G4UniformRand() > rej) ;
00216 
00217   // compute angle from t
00218   G4double sinth2 =
00219              0.5*(t-tmin)/(2.*(TotalEnergy*(TotalEnergy-ep)-Mass*Mass)-tmin);
00220   G4double theta = std::acos(1. - 2.*sinth2);
00221 
00222   G4double phi = twopi*G4UniformRand();
00223   G4double sinth = std::sin(theta);
00224   G4double dirx = sinth*std::cos(phi);
00225   G4double diry = sinth*std::sin(phi);
00226   G4double dirz = std::cos(theta);
00227   G4ThreeVector finalDirection(dirx,diry,dirz);
00228   G4ThreeVector ParticleDirection(aTrack.Get4Momentum().vect().unit() );
00229   finalDirection.rotateUz(ParticleDirection);
00230 
00231   G4double NewKinEnergy = KineticEnergy - ep;
00232   G4double finalMomentum = std::sqrt(NewKinEnergy*(NewKinEnergy+2.*Mass) );
00233   G4double Ef = NewKinEnergy + Mass;
00234   G4double initMomentum = std::sqrt(KineticEnergy*(TotalEnergy+Mass) );
00235 
00236   // Set energy and direction of scattered primary in theParticleChange
00237   theParticleChange.SetStatusChange(isAlive);
00238   theParticleChange.SetEnergyChange(NewKinEnergy);
00239   theParticleChange.SetMomentumChange(finalDirection);
00240 
00241   // Now create the emitted gamma 
00242   G4LorentzVector primaryMomentum(initMomentum*ParticleDirection, TotalEnergy);
00243   G4LorentzVector fsMomentum(finalMomentum*finalDirection, Ef);
00244   G4LorentzVector momentumTransfer = primaryMomentum - fsMomentum;
00245 
00246   G4DynamicParticle* gamma = 
00247            new G4DynamicParticle(G4Gamma::Gamma(), momentumTransfer);
00248  
00249   return gamma;
00250 }
00251 
00252 
00253 void
00254 G4MuonVDNuclearModel::CalculateHadronicVertex(G4DynamicParticle* incident,
00255                                               G4Nucleus& target)
00256 {
00257   G4HadFinalState* hfs = 0;
00258   G4double gammaE = incident->GetTotalEnergy();
00259 
00260   if (gammaE < 10*GeV) {
00261     G4HadProjectile projectile(*incident);
00262     hfs = bert->ApplyYourself(projectile, target);
00263   } else {
00264     // convert incident gamma to a pi0
00265     G4double piMass = G4PionZero::PionZero()->GetPDGMass();
00266     G4double piKE = incident->GetTotalEnergy() - piMass;
00267     G4double piMom = std::sqrt(piKE*(piKE + 2*piMass) );
00268     G4ThreeVector piMomentum(incident->GetMomentumDirection() );
00269     piMomentum *= piMom;
00270     G4DynamicParticle theHadron(G4PionZero::PionZero(), piMomentum);
00271     G4HadProjectile projectile(theHadron);
00272     hfs = ftfp->ApplyYourself(projectile, target);
00273   }
00274 
00275   delete incident;
00276 
00277   // Copy secondaries from sub-model to model
00278   theParticleChange.AddSecondaries(hfs);
00279 } 
00280 
00281 
00282 void G4MuonVDNuclearModel::MakeSamplingTable()
00283 {
00284   G4double adat[] = {1.01,9.01,26.98,63.55,238.03};
00285   G4double zdat[] = {1.,4.,13.,29.,92.};
00286   G4int nzdat = 5;
00287 
00288   G4double tdat[] = {1.e3,1.e4,1.e5,1.e6,1.e7,1.e8,1.e9,1.e10};
00289   G4int ntdat = 8;
00290 
00291   G4int nbin;
00292   G4double KineticEnergy;
00293   G4double TotalEnergy;
00294   G4double Maxep;
00295   G4double CrossSection;
00296 
00297   G4double c;
00298   G4double y;
00299   G4double ymin,ymax;
00300   G4double dy,yy;
00301   G4double dx,x;
00302   G4double ep;
00303 
00304   G4double AtomicNumber;
00305   G4double AtomicWeight;
00306 
00307   for (G4int iz = 0; iz < nzdat; iz++) {
00308     AtomicNumber = zdat[iz];
00309     AtomicWeight = adat[iz]*(g/mole);  
00310 
00311     for (G4int it = 0; it < ntdat; it++) {
00312       KineticEnergy = tdat[it];
00313       TotalEnergy = KineticEnergy + G4MuonMinus::MuonMinus()->GetPDGMass();
00314       Maxep = TotalEnergy - 0.5*proton_mass_c2;
00315 
00316       CrossSection = 0.0;
00317 
00318       // Calculate the differential cross section
00319       // numerical integration in log .........
00320       c = std::log(Maxep/CutFixed);
00321       ymin = -5.0;
00322       ymax = 0.0;
00323       dy = (ymax-ymin)/NBIN; 
00324 
00325       nbin=-1;              
00326 
00327       y = ymin - 0.5*dy;
00328       yy = ymin - dy;
00329       for (G4int i = 0; i < NBIN; i++) {
00330         y += dy;
00331         x = std::exp(y);
00332         yy += dy;
00333         dx = std::exp(yy+dy)-std::exp(yy);
00334       
00335         ep = CutFixed*std::exp(c*x);
00336 
00337         CrossSection +=
00338            ep*dx*muNucXS.ComputeDDMicroscopicCrossSection(KineticEnergy,
00339                                                           AtomicNumber,
00340                                                           AtomicWeight, ep);
00341         if (nbin < NBIN) {
00342           nbin += 1;
00343           ya[nbin] = y;
00344           proba[iz][it][nbin] = CrossSection;
00345         }
00346       }
00347       ya[NBIN] = 0.;
00348  
00349       if (CrossSection > 0.0) {
00350         for (G4int ib = 0; ib <= nbin; ib++) proba[iz][it][ib] /= CrossSection;
00351       }
00352     } // loop on it
00353   } // loop on iz
00354 
00355   // G4cout << " Kokoulin XS = "
00356   //       << muNucXS.ComputeDDMicroscopicCrossSection(1*GeV, 20.0, 40.0*g/mole, 0.3*GeV)/millibarn
00357   //       << G4endl; 
00358 }
00359 

Generated on Mon May 27 17:48:55 2013 for Geant4 by  doxygen 1.4.7