// Idea behind this macro is to implement cuts which are functions of
// several ntuple variables.
// 
// Tracy Usher (July 18, 2001)

#ifndef JoseCutClass_c
#define JoseCutClass_c
#ifndef nTupleVars_c
#include "nTupleVars.cxx"
#endif
#ifndef PlotGraphs_c
#include "PlotGraphs.cxx"
#endif
#include "iostream.h"
#include "math.h"
#endif

class JoseCutClass
{
private:
  PlotGraphs* pGraph;

  TH1F*       pNoTrks;
  TH1F*       pRawPSF;
  TH1F*       pCutPSF;
  TH1F*       pTAng;
  TH1F*       pKink;
  TH1F*       pCutVal;

  int         nEvents;

  //Private functions for this class
  Float_t JoseSelCut(bool top, Float_t kalEne, Float_t zDir);
  Float_t fixCorrEnergy(Float_t CalEnergy, Float_t Conv_Twr_Dist,
                        Float_t gamZdir,   Float_t nFitGaps,
                        Float_t nFitHits   Float_t ShwrHits1,
                        Float_t ShwrHits2);

public:
  JoseCutPlots();
  ~JoseCutPlots();

  void fillHists(nTupleVars* pTupleVals);
  void plotHists();
};

JoseCutClass::JoseCutClass()
{
  printf("--In JoseCutClass constructor\n");

  pGraph = new PlotGraphs("ComplexCuts", "Complex Cuts", "Test of multivariable Cuts", 6);

  //********************************************************************
  //Set up logarithmic binning in energy
  int     decades = 5;
  int     iBin    = 0;
  Float_t binEdge = 0.001;
  Float_t binInc  = 0.001;
  Float_t xLowEdge[51];

  xLowEdge[iBin++] = binEdge;

  while(decades--)
  {
    int nIncrements = 9;
    
    while(nIncrements--)
    {
      binEdge += binInc;
      xLowEdge[iBin++] = binEdge;
    }

    binInc *= 10.;
  }
  xLowEdge[iBin] = binEdge + 1.;
  
  //********************************************************************
  //Define and fill the histogram of the # of events reconstructing a track
  pNoTrks = new TH1F("cntEm",  "Tracks",  25,  0., 25.00);
  pRawPSF = new TH1F("mcGam",  "PSF",    100,  0.,  0.25);
  pCutPSF = new TH1F("zRes",   "PSF",    100,  0.,  0.25);
  pCutVal = new TH1F("cutV",   "CutVal", 100,  0.,  0.2);
  pTAng   = new TH1F("tAng",   "tAngle", 100,  0.,  0.1);
  pKink   = new TH1F("Kink",   "Kink",   200, -.2,  0.2);


  nEvents = 0;

  printf("--Done with JoseCutClass constructor\n");

  return;
}

void JoseCutClass::fillHists(nTupleVars* pTupleVals)
{
  if (pTupleVals)
  {
    int     numTrks = pTupleVals->getNumTrks();
    int     fitType = pTupleVals->getFitType();
    Float_t tAngle  = pTupleVals->getTAngle();
    Float_t fitKink = pTupleVals->getFitKink();
    Float_t calEne  = pTupleVals->getCalEne();

    if (numTrks>0 && fitType>0 && fabs(fitKink) < 1. && tAngle < .2 && calEne > 20.) 
    {
      bool    tempVar = 1;
      Float_t logEne  = 0.00001;
      Float_t gamErr  = pTupleVals->getMcGamErr();
      Float_t zDir    = pTupleVals->getZDir();
      Float_t corEne  = pTupleVals->getCorEne();

      //First, fix the corrected energy (or try to)
      //kalEne = fixCorrEnergy(0.001*calEne, Rec_Conv_Twr_Dist, zDir, Tkr_Fit_Gaps, Tkr_No_Hits, Rec_showerHits1, Rec_showerHits2);
      //printf("  rawEne=%f, kalEne=%f\n", 0.001*calEne, corEne);
      //kalEne = 0.001* calEne;
      //kalEne = mcGamEne;

      if (corEne > 0.) {logEne = log10(corEne);}

      pNoTrks->Fill(numTrks);
      pRawPSF->Fill(gamErr);

      Float_t cutVal = JoseSelCut(tempVar, corEne, zDir);

      pCutVal->Fill(cutVal);
      pTAng->Fill(tAngle);
      pKink->Fill(fitKink);

      if (fabs(fitKink) <= cutVal && tAngle <= 0.1*cutVal)
      {
        pCutPSF->Fill(gamErr);
      }
    }

    nEvents++;
  }
  return;
}

void JoseCutClass::plotHists()
{
  int numPassCuts = pNoTrks->GetEntries();
  printf("# events: %i, # nonzero entries: %i\n", nEvents, numPassCuts);

  int nRawPSF     = pRawPSF->GetEntries();
  int nCutPSF     = pCutPSF->GetEntries();
  printf("# in raw PSF distribution: %i, # in cut PSF distribution: %i\n", nRawPSF, nCutPSF);
   
  //********************************
  //Ok, time to start drawing stuff
  pGraph->setPlot();
  pGraph->drawPlots(0, pNoTrks, 0);
  pGraph->drawPlots(1, pRawPSF, 0);
  pGraph->drawPlots(2, pCutPSF, 0);
  pGraph->drawPlots(3, pCutVal, 0);
  pGraph->drawPlots(4, pTAng,   0);
  pGraph->drawPlots(5, pKink,   0);

  return;
}

Float_t JoseCutClass::JoseSelCut(bool top, Float_t kalEne, Float_t zDir)
{
  Float_t tolerance = -1;

  if (kalEne > 0)
  {
    //Determine the energy and angle dependent tolerance
    Float_t kalFit;

    if (!top) {kalFit = (2.8 /kalEne) + (2.1 /sqrt(kalEne));}
    else      {kalFit = (6.35/kalEne) + (3.75/sqrt(kalEne));}

    kalFit = - (0.001*kalFit) / zDir;

    tolerance = 3.5 * kalFit;
  }

  return tolerance;
}

Float_t JoseCutClass::fixCorrEnergy(Float_t CalEnergy,
                                    Float_t Conv_Twr_Dist,
                                    Float_t gamZdir,
                                    Float_t nFitGaps,
                                    Float_t nFitHits,
                                    Float_t ShwrHits1,
                                    Float_t ShwrHits2)
{

    Float_t Corr_Energy = CalEnergy;

    if (Conv_Twr_Dist > 0)
    {
        //Energy in Hits
        Float_t hit_energy = 0.0003*ShwrHits1 + .0013*ShwrHits2;  
        //Number of layers crossed
        Float_t num_Layers_Crossed = (nFitGaps + nFitHits) / 2.;
        // Take out observed layer correlation
        Float_t x_ing_energy = .0016*num_Layers_Crossed; 
        // Edge of Tower correction 
        Float_t edg_corr     =     (Conv_Twr_Dist < 10.) ? 1 : .68 + (.069  - .0036* Conv_Twr_Dist )*Conv_Twr_Dist; 
        // Total Tracker energy 
        Float_t trk_energy = hit_energy + x_ing_energy;  

        //printf("--in fixCorrEnergy:\n");
        //printf("  Conv_Twr_Dist=%f, edg_corr=%f, gamZdir=%f\n",Conv_Twr_Dist,edg_corr,gamZdir);
        //printf("  CalEnergy=%f, trk_energy=%f\n", CalEnergy, trk_energy);
    
        trk_energy /= (fabs(gamZdir) > .2) ? fabs(gamZdir) : .2;
        if (edg_corr != 0.) {Corr_Energy = (CalEnergy + trk_energy)/edg_corr;}
    }

    return Corr_Energy;
}