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

/*! @file TrkDedxHeCut.h  The TrkDedxHeCut.h class definition file. */

#ifndef TRKDEDXHECUT_H_
#define TRKDEDXHECUT_H_

#include "../../PamCutBase/PamCutBase.h"

#include <TH2F.h>
#include <TF1.h>
#include <TStyle.h>
#include <TMath.h>


/*! @brief Cutting a band DEDX vs. Rigidity for Helium.
 * This cut discards events where DEDX vs. Rigidity by Tracker measurements are outside the helium "corridor".
 */
class TrkDedxHeCut: public PamCut {

public:
  /*! @brief Constructor.
   *
   * @param cutName The cut's name.
   */
  TrkDedxHeCut(const char *cutName) :
    PamCut(cutName){
    _cuthededxrig = BBbetagamma(4*.938, 2, 14, 28.0855, 173, 0, 500, 4.44, .1492, 3.3, .2014, 2.87,1.); //Bethe-Bloch for He in Si


#ifdef DEBUGPAMCUT

#define RIG_XBIN 600
#define DEDX_YBIN 600

  Double_t RIB[RIG_XBIN+1];
  for(Int_t i = 0; i<RIG_XBIN+1; i++) RIB[i]=(0.01)*TMath::Power(10.,i/(Double_t(RIG_XBIN))*5.);
  h_dedx_rig_before = new TH2F("DEDX_RIG_before","DEDX(Trk) vs. Rigidity distribution",
			       RIG_XBIN,RIB, DEDX_YBIN, 0., 30.);
  h_dedx_rig_before->GetYaxis()->SetTitle("Trk DE/DX (MIP)");
  h_dedx_rig_before->GetYaxis()->SetLabelSize(0.025);
  h_dedx_rig_before->GetXaxis()->SetTitle("Rigidity (GV)");
  h_dedx_rig_before->GetXaxis()->SetLabelSize(0.025);
  h_dedx_rig_before->SetContour(50);
  h_dedx_rig_before->SetStats(kFALSE);
  gStyle->SetPalette(1);
  Double_t sum = h_dedx_rig_before->GetBinContent(h_dedx_rig_before->GetMaximumBin());
  Double_t lev[50];
  for(Int_t i = 0; i<50; i++){
    lev[i]=TMath::Power(10.,i/(50.)*TMath::Log10(sum));
    h_dedx_rig_before->SetContourLevel(i,lev[i]);
  }
  h_dedx_rig_after = (TH2F*)h_dedx_rig_before->Clone("DEDX_RIG_after");
#endif

  }
  /*! @brief Destructor. */
  ~TrkDedxHeCut() {

    delete _cuthededxrig;

#ifdef DEBUGPAMCUT
    h_dedx_rig_before->Write();
    h_dedx_rig_after->Write();
#endif
  }

  /*! @brief The DEDX vs. beta Helium check.
   *
   * @param event The event to analyze.
   * @return #CUTOK if event is inside helium "corridor".
   * @return 0 otherwise.
   */
  int Check(PamLevel2 *event);

private:

  TF1* _cuthededxrig;

#ifdef DEBUGPAMCUT
  TH2F* h_dedx_rig_before;
  TH2F* h_dedx_rig_after;
#endif


  TF1* 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);
  Bool_t CutHeDEDXrig(Float_t rig, Float_t dedx, Float_t factor, TF1* _cuthededxrig);

};
#endif /* TRKDEDXHECUT_H_  */