text: //////////////////////////////////////////////////////////
  // PetrukhinModel.cc Class:
  //
  // Improvements: Function of Muom Brem using  nuclear screening correction
  // Description: Muon bremsstrahlung using the Petrukhin's model in FASTSIM
  // Authors: Sandro Fonseca de Souza and Andre Sznajder (UERJ/Brazil)
  // Date: 23-Nov-2010
  ////////////////////////////////////////////////////////////////////
  
  #include <fstream>
  #include "TF1.h"
  #include "FastSimulation/MaterialEffects/interface/PetrukhinModel.h"
  #include <math.h>
  using namespace std;
  
  //=====================================================================
  
  ///////////////////////////////////////////////////
  //Function of Muom Brem using  nuclear-electron screening correction from G4 style
  //
  
   double PetrukhinFunc (double *x, double *p ){
   
     //Function independent variable 
     double nu = x[0]; //fraction of muon's energy transferred to the photon
     
    // Parameters
     double E=p[0]; //Muon Energy (in GeV)
     double A=p[1];// Atomic weight
     double Z=p[2];// Atomic number
  
   /*
  
  //////////////////////////////////////////////////
  //Function of Muom Brem using  nuclear screening correction
  //Ref: http://pdg.lbl.gov/2008/AtomicNuclearProperties/adndt.pdf
  
     //Physical constants
     double B = 182.7;
     double ee = sqrt(2.7181) ; // sqrt(e)
     double ZZ=  pow( Z,-1./3.); // Z^-1/3
     ///////////////////////////////////////////////////////////////////
    double emass = 0.0005109990615;  // electron mass (GeV/c^2)
    double mumass = 0.105658367;//mu mass  (GeV/c^2)
  
     double re = 2.817940285e-13;// Classical electron radius (Units: cm)
     double alpha = 1./137.03599976; // fine structure constant
     double Dn = 1.54* (pow(A,0.27));
     double constant =  pow((2.0 * Z * emass/mumass * re ),2.0);
     //////////////////////////////////////////////
     
     double delta = (mumass * mumass * nu) /(2.* E * (1.- nu)); 
      
     double Delta_n = TMath::Log(Dn / (1.+ delta *( Dn * ee -2.)/ mumass)); //nuclear screening correction 
      
     double Phi = TMath::Log((B * mumass * ZZ / emass)/ (1.+ delta * ee * B * ZZ  / emass)) - Delta_n;//phi(delta)
     
      //Diff. Cross Section for Muon Brem from a screened nuclear (Equation 16: REF: LBNL-44742)
     double f = alpha * constant *(4./3.-4./3.*nu + nu*nu)*Phi/nu;
  */
  
  //////////////////////////////////////////////////
  // Function for Muon Brem Xsec from G4
  //////////////////////////////////////////////////
  //Physical constants
     double B = 183.;
     double Bl= 1429.;
     double ee = 1.64872 ; // sqrt(e)
     double Z13=  pow( Z,-1./3.); // Z^-1/3
     double Z23=  pow( Z,-2./3.); // Z^-2/3
  
     //Original values of paper
     double emass = 0.0005109990615;  // electron mass (GeV/c^2)
     double mumass = 0.105658367;     // muon mass  (GeV/c^2)
     // double re = 2.817940285e-13;     // Classical electron radius (Units: cm)
     double alpha =  0.00729735;      // 1./137.03599976; // fine structure constant
     double constant = 1.85736e-30;   // pow( ( emass / mumass * re ) , 2.0);
  
     double Dn = 1.54*(pow(A,0.27));
     double Dnl= pow(Dn,(1.-1./Z));
  
     double delta = (mumass * mumass * nu)/(2.* E * (1.- nu)); 
      
      
     double Phi_n = TMath::Log( B * Z13 *( mumass + delta * ( Dnl * ee -2 ))
                    / ( Dnl * ( emass + delta * ee * B * Z13 ) ) );
  
      
     double Phi_e = TMath::Log( ( Bl * Z23 * mumass )
                    / ( 1.+ delta * mumass / ( emass*emass * ee ) )
                    / ( emass + delta *ee * Bl * Z23 ) );
  
  
  
  //Diff. Cross Section for Muon Brem from G4 ( without NA/A factor )
     double f =  16./3. * alpha * constant* Z * ( Z * Phi_n + Phi_e ) * (1./nu) * (1. - nu + 0.75*nu*nu) ;
  
     return f;
   
  
  
   }
  
  
  
  