const double Pi     = 3.14159265359;
const double c      = 299792.458; // en km s^-1
const double c2     = c*c;
const double hbar   = 6.58211928e-25; // en GeV s
const double N0     = 6.02214179e23;

/* en [kg^-1 day^-1 keV^-1] */
Double_t dmFlux(const Double_t *X /* X[0]=ER_keV en keV */, const Double_t *par /* par[0]=M en GeV c^-2, par[1]=sigma en cm^2 */){

  const double ER_keV = X[0]; // en keV
  const double ER     = ER_keV/1e6; // en GeV
  const double M      = par[0]/c2;// en GeV c^-2 
  const double sigma  = par[1];// en cm^2
  
  const double t      = 0   ; //days 
  const double v0     = 220 ; //km s^-1 
  const double vesc   = 650 ; //km s^-1 
  const double vE     = 244  + 15*cos(2*Pi* (t-152.5)/365.25); // km s^-1
  const double rho    = 0.3/c2 ; //GeV c^-2 cm^-3
  const double A      = 28  ; //silicon
  const double MN     = 26.0603162/c2  ; //GeV c^-2 : Nucleo
  const double Mn     = 0.938272/c2    ; //GeV c^-2 : nucleon
  

  // calcula el radio del nuclo que se va a usar en el factor de forma
  double fc = 1.23*pow(A,1./3.)-0.60;
  double fa = 0.52;
  double fs = 0.9;
  double rn = sqrt( fc*fc + 7./3.*Pi*Pi*fa*fa - 5*fs*fs); // en fm
  
  
  // calcula el factor de forma para corregir la seccion eficaz por perdida de coherencia
  double qrn = 6.92e-3*sqrt(A*ER_keV)*(rn); //dimensionless
  double F   = (qrn>0)? 3. * ( sin(qrn) - qrn*cos(qrn) )/pow(qrn,3) * exp( -pow(qrn*fs/rn,2)/2.) : 1.;

  
  // calcula los parametros necesarios para obtener el rate
  double m_N = M * MN / (M+MN); //masa reducida DM-Nucleo
  double m_n = M * Mn / (M+Mn); //masa reducida DM-nucleon
  double sigma0 = sigma * m_N*m_N/m_n/m_n * A*A ;   //cm^2
  double sigmaF = sigma0 * F*F ;   // seccion eficaz corregida por el factor de forma (longitud de onda h/q no tan grande comparada con el radio nuclear)
  double sigmaAux = sigmaF;
  const double v02  = v0*v0;
  double k0 = pow( Pi*v02, 3./2.);
  double k1 = k0*( TMath::Erf(vesc/v0) - 2/sqrt(Pi) * vesc/v0 * exp(-vesc/v0 * vesc/v0) );
  double n0 = rho/M;
  double R0 = 2./sqrt(Pi) * N0/A * n0 * v0 * sigmaAux;
  double E0 = M * v02/2.;
  double r  = 4 * M * MN / pow(M+MN,2);  
  double vmin = v0 * sqrt(ER/E0/r);
  

  double conFactor = 8640*1000; // para ir de [km cm GeV^-1 s^-1 g^-1] a [kg^-1 day^-1 keV^-1] 
  
  double dRdER = 0;
  if( vmin <= vesc-vE )
    dRdER = conFactor * k0/k1 * R0/E0/r * (sqrt(Pi)/4. * v0/vE * ( TMath::Erf( (vmin+vE)/v0 )-TMath::Erf( (vmin-vE)/v0 ) ) - exp(-vesc*vesc/v02) );
  else if( vmin <= vesc+vE )
    dRdER = conFactor * k0/k1 * R0/E0/r * (sqrt(Pi)/4. * v0/vE * ( TMath::Erf( (vesc)/v0 )-TMath::Erf( (vmin-vE)/v0 ) ) - (vE+vesc-vmin)/vE/2.*exp(-vesc*vesc/v02) );
  
  return dRdER; /* en [kg^-1 day^-1 keV^-1] */
}


/* En este ejempolo se crea un objeto TF1
   para el rate diferencial dR/dER */
void rateDM(){
  
  const float minE = 1e-3; // en keV
  const float maxE = 100; // en keV
  
  TF1 fRateDM("fRateDM",dmFlux, minE,maxE,2);

  const Double_t par[]={ 10, // Mass: GeV c^-2
                         1e-39   // sigma: cm^  
                       };
  
  fRateDM.SetParameters(par);
  fRateDM.SetNpx(10000);
  TCanvas canvas;
  // canvas.SetLogx();
  // canvas.SetLogy();
  fRateDM.Draw();
  canvas.WaitPrimitive();
}