/*
 * TrkDedxHeCut.cpp
 *
 *  Created on: 19-may-2009
 *      Author: N. Nikonov
 */

/*! @file TrkDedxHeCut.cpp The TrkDedxHeCut class implementation file. */

#include "TrkDedxHeCut.h"

int TrkDedxHeCut::Check(PamLevel2 *event) {

  TrkTrack* trk2 = event->GetTrack(0)->GetTrkTrack();
  Double_t rig = trk2->GetRigidity();
  Double_t dedxtr = trk2->GetDEDX();

#ifdef DEBUGPAMCUT
  h_dedx_rig_before->Fill(rig, dedxtr);
#endif

  if ( !CutHeDEDXrig(rig,dedxtr,1.15,_cuthededxrig) ||
       CutHeDEDXrig(rig,dedxtr,0.5,_cuthededxrig)) return 0;

#ifdef DEBUGPAMCUT
  h_dedx_rig_after->Fill(rig, dedxtr);
#endif

  return CUTOK;
}


TF1* TrkDedxHeCut::BBbetagamma(Float_t m_GeV, Float_t z, Int_t Z, Float_t A, Float_t I_eV, Double_t xmin, Double_t xmax, Float_t C_positive, Float_t a, Float_t m, Float_t x0, Float_t x1, Float_t scale)
{
  TF1* f;
  if(C_positive != 0.) {
    // do density correction
    // x<x_0 ( delta = 0)
    // I don't need a function in this range
    // x0<x<x1
    // I build the function as a string.
    TString delta2 = "(4.6052*log([1]*x/[0])/log(10)-[5]+[6]*pow(([9]-log([1]*x/[0])/log(10)),[7])) * (([8] <= log([1]*x/[0])) && ( log([1]*x/[0])/log(10)  < [9] ) )";
    // x>x1
    TString delta3 = "(4.6052*log([1]*x/[0])/log(10)-[5]) * (log([1]*x/[0])/log(10) >= [9])";

    TString fstring = "[10]*(.3071 * ([2]/[3]) * [1]^2 / ([0]^2/([1]*x)^2+1)^-1 * ( 1./2 * log( (2 * 511e3 * ([1]*x/[0])^2 )^2 ) - log([4]) - ([0]^2/([1]*x)^2+1)^-1 - (";
    fstring += delta2.Data();
    fstring += "+";
    fstring += delta3.Data();
    fstring += ")/2))";
    f = new TF1("bethebloch", fstring.Data());
  }
  else
    f = new TF1("bethebloch", ".3071 * ([2]/[3]) * [1]^2 / ([0]^2/([1]*x)^2+1)^-1 * ( 1./2 * log( (2 * 511e3 * ([1]*x/[0])^2 )^2 ) - log([4]) - ([0]^2/([1]*x)^2+1)^-1)");
  f->SetParameters(m_GeV, z, Z, A, I_eV, C_positive, a, m, x0, x1,scale);

  return f;
}


Bool_t TrkDedxHeCut::CutHeDEDXrig(Float_t rig, Float_t dedx,
					 Float_t factor, TF1* _cuthededxrig){
  if(factor != 1.) dedx = dedx*factor;
  if(rig < 0.2)  return kFALSE;
  if(dedx > 0.55*_cuthededxrig->Eval(rig)) return kTRUE;
  else return kFALSE;
}