MaterialEffects/src/PetrukhinModel.cc

0001 //////////////////////////////////////////////////////////
0002 // PetrukhinModel.cc Class:
0003 //
0004 // Improvements: Function of Muom Brem using  nuclear screening correction
0005 // Description: Muon bremsstrahlung using the Petrukhin's model in FASTSIM
0006 // Authors: Sandro Fonseca de Souza and Andre Sznajder (UERJ/Brazil)
0007 // Date: 23-Nov-2010
0008 ////////////////////////////////////////////////////////////////////
0009
0010 #include <fstream>
0011 #include "TF1.h"
0012 #include "FastSimulation/MaterialEffects/interface/PetrukhinModel.h"
0013 #include <cmath>
0014 using namespace std;
0015
0016 //=====================================================================
0017
0018 ///////////////////////////////////////////////////
0019 //Function of Muom Brem using  nuclear-electron screening correction from G4 style
0020 //
0021
0022 double PetrukhinFunc(double *x, double *p) {
0023   //Function independent variable
0024   double nu = x[0];  //fraction of muon's energy transferred to the photon
0025
0026   // Parameters
0027   double E = p[0];  //Muon Energy (in GeV)
0028   double A = p[1];  // Atomic weight
0029   double Z = p[2];  // Atomic number
0030
0031   /*
0032
0033 //////////////////////////////////////////////////
0034 //Function of Muom Brem using  nuclear screening correction
0035 //Ref: http://pdg.lbl.gov/2008/AtomicNuclearProperties/adndt.pdf
0036
0037    //Physical constants
0038    double B = 182.7;
0039    double ee = sqrt(2.7181) ; // sqrt(e)
0040    double ZZ=  pow( Z,-1./3.); // Z^-1/3
0041    ///////////////////////////////////////////////////////////////////
0042   double emass = 0.0005109990615;  // electron mass (GeV/c^2)
0043   double mumass = 0.105658367;//mu mass  (GeV/c^2)
0044
0045    double re = 2.817940285e-13;// Classical electron radius (Units: cm)
0046    double alpha = 1./137.03599976; // fine structure constant
0047    double Dn = 1.54* (pow(A,0.27));
0048    double constant =  pow((2.0 * Z * emass/mumass * re ),2.0);
0049    //////////////////////////////////////////////
0050
0051    double delta = (mumass * mumass * nu) /(2.* E * (1.- nu));
0052
0053    double Delta_n = TMath::Log(Dn / (1.+ delta *( Dn * ee -2.)/ mumass)); //nuclear screening correction
0054
0055    double Phi = TMath::Log((B * mumass * ZZ / emass)/ (1.+ delta * ee * B * ZZ  / emass)) - Delta_n;//phi(delta)
0056
0057     //Diff. Cross Section for Muon Brem from a screened nuclear (Equation 16: REF: LBNL-44742)
0058    double f = alpha * constant *(4./3.-4./3.*nu + nu*nu)*Phi/nu;
0059 */
0060
0061   //////////////////////////////////////////////////
0062   // Function for Muon Brem Xsec from G4
0063   //////////////////////////////////////////////////
0064   //Physical constants
0065   double B = 183.;
0066   double Bl = 1429.;
0067   double ee = 1.64872;            // sqrt(e)
0068   double Z13 = pow(Z, -1. / 3.);  // Z^-1/3
0069   double Z23 = pow(Z, -2. / 3.);  // Z^-2/3
0070
0071   //Original values of paper
0072   double emass = 0.0005109990615;  // electron mass (GeV/c^2)
0073   double mumass = 0.105658367;     // muon mass  (GeV/c^2)
0074   // double re = 2.817940285e-13;     // Classical electron radius (Units: cm)
0075   double alpha = 0.00729735;      // 1./137.03599976; // fine structure constant
0076   double constant = 1.85736e-30;  // pow( ( emass / mumass * re ) , 2.0);
0077
0078   double Dn = 1.54 * (pow(A, 0.27));
0079   double Dnl = pow(Dn, (1. - 1. / Z));
0080
0081   double delta = (mumass * mumass * nu) / (2. * E * (1. - nu));
0082
0083   double Phi_n = TMath::Log(B * Z13 * (mumass + delta * (Dnl * ee - 2)) / (Dnl * (emass + delta * ee * B * Z13)));
0084
0085   double Phi_e =
0086       TMath::Log((Bl * Z23 * mumass) / (1. + delta * mumass / (emass * emass * ee)) / (emass + delta * ee * Bl * Z23));
0087
0088   //Diff. Cross Section for Muon Brem from G4 ( without NA/A factor )
0089   double f = 16. / 3. * alpha * constant * Z * (Z * Phi_n + Phi_e) * (1. / nu) * (1. - nu + 0.75 * nu * nu);
0090
0091   return f;
0092 }