PetrukhinModel.cc

Go to the documentation of this file.

// 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 <cmath>
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;
}