src/cpp/jail/resonanceDecay/functions.cpp
// This function interpolates the needed spectra for a given y, pt and phi.

This function interpolates the needed spectra for a given y, pt and phi.
This function interpolates the needed spectra for a given y pt and phi .
double Edndp3(double yr, double ptr, double phirin, int res_num)
{
  if (phirin < 0.0) { printf("ERROR: phir %15.8le < 0 !!! \n", phirin); exit(0); }
  if (phirin > 2.0*PI) { printf("ERROR: phir %15.8le > 2PI !!! \n", phirin); exit(0); }
  double phir = phirin;
  int pn = partid[MHALF + res_num];
  if (ptr > particle[pn].pt[particle[pn].npT - 1]) return 0.;
  if (fabs(yr) > particle[pn].ymax && !BOOST_INV) return 0.;
  int nphi = 1;
  while ( (phir > PHI[nphi]) && ( nphi < (particle[pn].nphi-1) ) ) nphi++;
  int npT = 1;
  while (ptr > particle[pn].pt[npT] && npT<(particle[pn].npT - 1)) npT++;
  int ny = 1;
  while ( yr > particle[pn].y[ny] && ny < (particle[pn].ny - 1) ) ny++;
  double f1 = lin_int(PHI[nphi-1], PHI[nphi], particle[pn].dNdypTdpTdphi[ny-1][npT-1][nphi-1], particle[pn].dNdypTdpTdphi[ny-1][npT-1][nphi], phir);
  double f2 = lin_int(PHI[nphi-1], PHI[nphi], particle[pn].dNdypTdpTdphi[ny-1][npT][nphi-1], particle[pn].dNdypTdpTdphi[ny-1][npT][nphi], phir);
  if (f1 < 0.) f1 = 1e-30;
  if (f2 < 0.) f2 = 1e-30;
  double f1s = f1;
  double f2s = f2;
  double val1 = lin_int(particle[pn].pt[npT-1], particle[pn].pt[npT], f1, f2, ptr);
  f1 = lin_int(PHI[nphi-1], PHI[nphi], particle[pn].dNdypTdpTdphi[ny][npT-1][nphi-1], particle[pn].dNdypTdpTdphi[ny][npT-1][nphi], phir);
  f2 = lin_int(PHI[nphi-1], PHI[nphi], particle[pn].dNdypTdpTdphi[ny][npT][nphi-1], particle[pn].dNdypTdpTdphi[ny][npT][nphi], phir);
  if (f1 < 0.) f1 = 1e-30;
  if (f2 < 0.) f2 = 1e-30;
  f1s = f1;
  f2s = f2;
  double val2 = lin_int(particle[pn].pt[npT-1], particle[pn].pt[npT], f1, f2, ptr);
  double val = lin_int(particle[pn].y[ny-1], particle[pn].y[ny], val1, val2, yr);
  return val;
}
double Edndp3 ( double yr , double ptr , double phirin , int res_num ) { if ( phirin < 0.0 ) { printf ( " \n " , phirin ) ; exit ( 0 ) ; } if ( phirin > 2.0 * PI ) { printf ( " \n " , phirin ) ; exit ( 0 ) ; } double phir = phirin ; int pn = partid [ MHALF + res_num ] ; if ( ptr > particle [ pn ] . pt [ particle [ pn ] . npT - 1 ] ) return 0. ; if ( fabs ( yr ) > particle [ pn ] . ymax && ! BOOST_INV ) return 0. ; int nphi = 1 ; while ( ( phir > PHI [ nphi ] ) && ( nphi < ( particle [ pn ] . nphi - 1 ) ) ) nphi ++ ; int npT = 1 ; while ( ptr > particle [ pn ] . pt [ npT ] && npT < ( particle [ pn ] . npT - 1 ) ) npT ++ ; int ny = 1 ; while ( yr > particle [ pn ] . y [ ny ] && ny < ( particle [ pn ] . ny - 1 ) ) ny ++ ; double f1 = lin_int ( PHI [ nphi - 1 ] , PHI [ nphi ] , particle [ pn ] . dNdypTdpTdphi [ ny - 1 ] [ npT - 1 ] [ nphi - 1 ] , particle [ pn ] . dNdypTdpTdphi [ ny - 1 ] [ npT - 1 ] [ nphi ] , phir ) ; double f2 = lin_int ( PHI [ nphi - 1 ] , PHI [ nphi ] , particle [ pn ] . dNdypTdpTdphi [ ny - 1 ] [ npT ] [ nphi - 1 ] , particle [ pn ] . dNdypTdpTdphi [ ny - 1 ] [ npT ] [ nphi ] , phir ) ; if ( f1 < 0. ) f1 = 1e-30 ; if ( f2 < 0. ) f2 = 1e-30 ; double f1s = f1 ; double f2s = f2 ; double val1 = lin_int ( particle [ pn ] . pt [ npT - 1 ] , particle [ pn ] . pt [ npT ] , f1 , f2 , ptr ) ; f1 = lin_int ( PHI [ nphi - 1 ] , PHI [ nphi ] , particle [ pn ] . dNdypTdpTdphi [ ny ] [ npT - 1 ] [ nphi - 1 ] , particle [ pn ] . dNdypTdpTdphi [ ny ] [ npT - 1 ] [ nphi ] , phir ) ; f2 = lin_int ( PHI [ nphi - 1 ] , PHI [ nphi ] , particle [ pn ] . dNdypTdpTdphi [ ny ] [ npT ] [ nphi - 1 ] , particle [ pn ] . dNdypTdpTdphi [ ny ] [ npT ] [ nphi ] , phir ) ; if ( f1 < 0. ) f1 = 1e-30 ; if ( f2 < 0. ) f2 = 1e-30 ; f1s = f1 ; f2s = f2 ; double val2 = lin_int ( particle [ pn ] . pt [ npT - 1 ] , particle [ pn ] . pt [ npT ] , f1 , f2 , ptr ) ; double val = lin_int ( particle [ pn ] . y [ ny - 1 ] , particle [ pn ] . y [ ny ] , val1 , val2 , yr ) ; return val ; }
This function interpolates the needed spectra for a given y, pt and phi.
This function interpolates the needed spectra for a given y pt and phi .
// double   yr;     /* y  of resonance */
 // double   ptr;        /* pt of resonance */
 // double   phirin;     /* phi angle  of resonance */
 // int  res_num;    /* Montecarlo number of resonance   */
 // if pseudofreeze flag is set, yr is the *pseudorapidity* of the resonance
 /* phi interpolation */ // security: if for some reason we got a negative number of particles
 // (happened in the viscous code at large eta sometimes)
 /*
  if (ptr > PTCHANGE && f1s > 0 && f2s > 0) {
    f1 = log(f1);
    f2 = log(f2);
  }
  */ //if (ptr > PTCHANGE && f1s > 0 && f2s > 0) val1 = exp(val1);
 // security: if for some reason we got a negative number of particles
 // (happened in the viscous code at large eta sometimes)
 /*
  if (ptr > PTCHANGE && f1s > 0 && f2s > 0)
  {
    f1 = log(f1);
    f2 = log(f2);
  }
  */ //if (ptr > PTCHANGE && f1s > 0 && f2s > 0) val2 = exp(val2);
