Float_t ProtonMomentumLoss(Float_t p,Float_t theta){

  //*********************************************************************
  //* Returns the momentum lost (in [GeV]) by a proton in passing from  *
  //* the target cell to the exit plane of the outer-most wire chamber  *
  //*                                                                   *
  //* Arguments:                                                        *
  //*   p     : reconstructed proton momentum in [GeV]                  *
  //*   theta : polar angle of reconstructed proton momentum in the     *
  //*           BLAST coordinate system in [deg]                        *
  //*********************************************************************

  const Float_t PLOW=0.135;
  const Float_t PUPP=1.000;

  const Int_t NUMDIVSTHETA=7;
  const Float_t THETALOW=15.0;
  const Float_t THETAUPP=85.0;

  const Float_t C[NUMDIVSTHETA][9]=
  {{+2.525887e-01,-2.463591e+00,+1.109465e+01,-2.902014e+01,+4.737158e+01,
    -4.889279e+01,+3.092569e+01,-1.088934e+01,+1.622587e+00},
   {+2.467519e-01,-2.586977e+00,+1.255536e+01,-3.550835e+01,+6.294246e+01,
    -7.098161e+01,+4.949543e+01,-1.946127e+01,+3.299227e+00},
   {+2.470371e-01,-2.782622e+00,+1.457769e+01,-4.467991e+01,+8.609354e+01,
    -1.057935e+02,+8.054456e+01,-3.464091e+01,+6.435032e+00},
   {+2.377869e-01,-2.742439e+00,+1.467808e+01,-4.586923e+01,+8.994726e+01,
    -1.122878e+02,+8.671446e+01,-3.777704e+01,+7.099781e+00},
   {+2.250083e-01,-2.584070e+00,+1.372053e+01,-4.242290e+01,+8.215929e+01,
    -1.011799e+02,+7.703141e+01,-3.307514e+01,+6.126511e+00},
   {+2.226852e-01,-2.578033e+00,+1.379630e+01,-4.299399e+01,+8.392735e+01,
    -1.041823e+02,+7.994897e+01,-3.459866e+01,+6.458445e+00},
   {+2.219288e-01,-2.600811e+00,+1.412805e+01,-4.477074e+01,+8.894891e+01,
    -1.124014e+02,+8.778393e+01,-3.864007e+01,+7.330892e+00}};

  Int_t thetaIndice=
    Int_t(NUMDIVSTHETA*(theta-THETALOW)/((THETAUPP-THETALOW)));

  if(p<PLOW) p=PLOW;
  else if(p>PUPP) p=PUPP;

  if(thetaIndice<0) thetaIndice=0;
  else if(thetaIndice>=NUMDIVSTHETA) thetaIndice=NUMDIVSTHETA-1;

  return(C[thetaIndice][0]+C[thetaIndice][1]*p+C[thetaIndice][2]*p*p
    +C[thetaIndice][3]*p*p*p+C[thetaIndice][4]*p*p*p*p
    +C[thetaIndice][5]*p*p*p*p*p+C[thetaIndice][6]*p*p*p*p*p*p
    +C[thetaIndice][7]*p*p*p*p*p*p*p+C[thetaIndice][8]*p*p*p*p*p*p*p*p);
}

Float_t DeuteronMomentumLoss(Float_t p,Float_t theta){

  //*********************************************************************
  //* Returns the momentum lost (in [GeV]) by a deuteron in passing     *
  //* from the target cell to the exit plane of the outer-most wire     *
  //* chamber                                                           *
  //*                                                                   *
  //* Arguments:                                                        *
  //*   p     : reconstructed proton momentum in [GeV]                  *
  //*   theta : polar angle of reconstructed proton momentum in the     *
  //*           BLAST coordinate system in [deg]                        *
  //*********************************************************************

  const Float_t PLOW=0.195;
  const Float_t PUPP=1.000;

  const Int_t NUMDIVSTHETA=7;
  const Float_t THETALOW=15.0;
  const Float_t THETAUPP=85.0;

  const Float_t C[NUMDIVSTHETA][9]=
  {{-1.809151e-01,+6.310025e+00,-4.794185e+01,+1.786171e+02,-3.884445e+02,
    +5.186335e+02,-4.198597e+02,+1.895051e+02,-3.663460e+01},
   {+4.858613e-01,-3.317892e+00,+9.801995e+00,-1.313674e+01,-2.502832e-01,
    +2.650557e+01,-3.757053e+01,+2.288091e+01,-5.395364e+00},
   {+7.957955e-01,-8.346428e+00,+4.296252e+01,-1.331760e+02,+2.630077e+02,
    -3.329085e+02,+2.613282e+02,-1.157758e+02,+2.211583e+01},
   {+2.973238e-01,-1.308478e+00,-1.976850e-01,+1.517645e+01,-4.863401e+01,
    +7.664054e+01,-6.768753e+01,+3.210503e+01,-6.388907e+00},
   {+1.227123e-01,+1.082067e+00,-1.457410e+01,+6.379820e+01,-1.492508e+02,
    +2.069895e+02,-1.709908e+02,+7.794567e+01,-1.511996e+01},
   {+2.906743e-01,-1.360304e+00,+3.116820e-01,+1.351816e+01,-4.580120e+01,
    +7.386021e+01,-6.613844e+01,+3.166455e+01,-6.342882e+00},
   {+2.983323e-01,-1.480130e+00,+9.941900e-01,+1.149800e+01,-4.235965e+01,
    +7.045492e+01,-6.431991e+01,+3.126104e+01,-6.344395e+00}};

  Int_t thetaIndice=
    Int_t(NUMDIVSTHETA*(theta-THETALOW)/((THETAUPP-THETALOW)));

  if(p<PLOW) p=PLOW;
  else if(p>PUPP) p=PUPP;

  if(thetaIndice<0) thetaIndice=0;
  else if(thetaIndice>=NUMDIVSTHETA) thetaIndice=NUMDIVSTHETA-1;

  return(C[thetaIndice][0]+C[thetaIndice][1]*p+C[thetaIndice][2]*p*p
    +C[thetaIndice][3]*p*p*p+C[thetaIndice][4]*p*p*p*p
    +C[thetaIndice][5]*p*p*p*p*p+C[thetaIndice][6]*p*p*p*p*p*p
    +C[thetaIndice][7]*p*p*p*p*p*p*p+C[thetaIndice][8]*p*p*p*p*p*p*p*p);
}

const float E=0.850, M=0.93827231, EM2=(E+M)*(E+M), deg=57.2957795131;

float pkp_(float twp)			// proton momentum
{
  if (twp>89.9) twp = 89.9;
  float xcp = cos(twp/deg); xcp*=xcp;	// cos(th_p)
  float xsp = 1-xcp;			// sin(th_p)
  float xa  = E*E*xcp-EM2;
  float xb  = 2*M*E*E*xcp + 2*E*xsp*EM2;// Ax^2+Bx+c
  float xc  = E*E*(M*M*xcp - xsp*EM2);
  float xd  = -xb/2/xa;
  float x   = xd - sqrt(xd*xd-xc/xa);	// E'*cos(th_e)
  if (x!=x) return 0;
  float Ep  = E*(M+x)/(M+E);		// E'
  return sqrt(E*E + Ep*Ep - 2*E*x);	// q = p_p
}

Double_t MomCor(Double_t* x, Double_t* p)
{
  return ProtonMomentumLoss(pkp_(*x),*x);
}


/// from Eugene


// E loss as a function of theta_p
Double_t ElossTH(Double_t *x, Double_t *p) {
  Double_t THp=*x;
  return = 0.0597905 - 0.00529057*THp + 1.73551e-04*THp*THp - 
    2.49832e-06*THp*THp*THp + 1.33493e-08*THp*THp*THp*THp;
}

//E loss as a function of Pp
Float_t ElossP(Float_t Pp) {
  Float_t p_lost = 0.0305683 - 0.164857*Pp + 0.334062*Pp*Pp - 
    0.299825*Pp*Pp*Pp + 0.100794*Pp*Pp*Pp*Pp;
  return p_lost;
}

void eloss_aaron()
{
  TF1* fa=new TF1("cor1",MomCor,38,68,0);
  fa->Draw();
  TF1* fe=new TF1("cor2",ElossTH,38,68,0);
  fe->SetLineColor(kRed);
  fe->Draw("same");
}