// 
// Class for analyzing the PSF function
// 
// Tracy Usher (July 8, 2001)

class PSFplot
{
 private:
  Double_t    contain68;
  Double_t    contain95;
  Double_t    cutEffic;

  TString*    pHistName;
  TString*    pSumName;
  TString*    pHistTitle;
  int         nBins;
  Axis_t      lowBinEdge;
  Axis_t      hiBinEdge;

  TH1F*       pHist;
  TH1F*       pSum;

  TGaxis*     axis;
  TPaveLabel* pLabel68;
  TPaveLabel* pLabel95;
  int         numTotal;

 public:
  PSFplot();
  PSFplot(int histId, char* pHstTtle, int numBins, Double_t loEdge, Double_t hiEdge, int nEntries);
  ~PSFplot();

  void fillHists(TTree* pTree, TCut cutVal);
  void fillHists(TTree* pTree, char* pFillVal, TCut cutVal);
  void plotHists(TCanvas* pCan, TPad* pPad); 

  Double_t getContain68()   {return contain68;} 
  Double_t getContain95()   {return contain95;}
  Double_t getEfficiency()  {return cutEffic;}
};

PSFplot::PSFplot()
{
  contain68  = 0.;
  contain95  = 0.;
  cutEffic   = -1.;

  pHistName  = 0;
  pSumName   = 0;
  pHistTitle = 0;
  nBins      = 0;
  lowBinEdge = 0.;
  hiBinEdge  = 0.;

  pHist      = 0;
  pSum       = 0;

  axis       = 0;
  pLabel68   = 0;
  pLabel95   = 0;
  numTotal   = 0;

  return;
}

PSFplot::PSFplot(int histId, char* pHstTtle, int numBins, Double_t loEdge, Double_t hiEdge, int nEntries)
{
  //Take care of histogram name first
  char histName[] = "gamErr01";
  sprintf(histName,"gamErr%-2i",histId);
  pHistName  = new TString(histName);

  sprintf(histName,"gamSum%-2i",histId);
  pSumName   = new TString(histName);

  //Copy over the histogram title
  pHistTitle = new TString(pHstTtle);

  //Copy over the rest of the stuff
  nBins      = numBins;
  lowBinEdge = loEdge;
  hiBinEdge  = hiEdge;

  //Zero out everything else
  contain68  = 0.;
  contain95  = 0.;
  cutEffic   = -1.;
  pHist      = 0;
  pSum       = 0;

  axis       = 0;
  pLabel68   = 0;
  pLabel95   = 0;
  numTotal   = nEntries;

  return;
}

PSFplot::~PSFplot()
{
  if (pHistName)  delete pHistName;
  if (PSumName)   delete pSumName;
  if (pHistTitle) delete pHistTitle;
  if (pHist)      delete pHist;
  if (pSum)       delete pSum;
  if (axis)       delete axis;
  if (pLabel68)   delete pLabel68;
  if (pLabel95)   delete pLabel95;

  return;
}

void PSFplot::fillHists(TTree* pTree, TCut cutVal)
{
  //Define a set of temporary histograms for finding the binning
  int      idx       =  0;
  int      numBins   = 500;
  Double_t xBinLow   =  0.;
  Double_t xBinHigh  =  2.;
  TH1F*    pTmpHis   = new TH1F("tmpHis","temp",numBins,xBinLow,xBinHigh);
  TH1F*    pTmpInteg = new TH1F("tmpInteg","temp",numBins,xBinLow,xBinHigh);

  //Project will fill but not draw the histogram
  pTree->Project("tmpHis", "MC_Gamma_Err", cutVal);

  //Now integrate this histogram
  Double_t intSum    =  0.;
  while(idx <= numBins+1)
  {
    intSum += pTmpHis->GetBinContent(idx);
    pTmpInteg->SetBinContent(idx,intSum);
    idx++;
  }

  //Make sure nothing bad is happening
  if (intSum == 0.)
  {
    printf("zero integral value found \n");
    delete pTmpHis;
    delete pTmpInteg;

    return;
  }

  //Now find 95% containment
  Double_t intScale = 1.0/intSum;
  pTmpInteg->Scale(intScale);
  idx =  0;
  while(idx <= numBins)
  {
    Double_t normVal = pTmpInteg->GetBinContent(idx);
    if (normVal > 0.95) break;
    idx++;
  }

  //Determine the new upper edge of the hist
  contain95 = pTmpInteg->GetBinCenter(idx+1);

  //Do we want to rescale the histogram?
  if (hiBinEdge <= lowBinEdge)
  {
    int nBinVal = (100*contain95+0.5);
    hiBinEdge   = 0.01 * nBinVal;
  }

  //Get rid of the temporary histograms
  delete pTmpHis;
  delete pTmpInteg;

  printf("*************** %s ***************\n",pHistTitle->Data());
  printf("--> Find the 95%% containment in bin %3i, center %f \n",idx,contain95);

  //Define the histogram to hold the PSF 
  pHist = new TH1F(pHistName->Data(),pHistTitle->Data(),nBins,lowBinEdge,hiBinEdge);

  pTree->Project(pHistName->Data(), "MC_Gamma_Err", cutVal);

  pSum = new TH1F(pSumName->Data(),pHistTitle->Data(),nBins,lowBinEdge,hiBinEdge);

  //Go through integrating again...
  idx    =  0;
  intSum =  0.;
  while(idx <= nBins+1)
  {
    intSum += pHist->GetBinContent(idx);
    pSum->SetBinContent(idx,intSum);
    idx++;
  }

  //If non zero integral then scale the integral hist
  if (intSum > 0.)
  {
    Double_t intScale = 1.0/intSum;
    pSum->Scale(intScale);
    pSum->SetLabelSize(0.07,"X"); 
    pSum->SetLabelSize(0.07,"Y");

    idx = nBins+1;
    while(idx--)
    {
      if (pSum->GetBinContent(idx) < 0.68) break;
    }
    //make sure we found something...
    if (idx < 1) idx = 1;

    contain68 = pSum->GetBinCenter(idx);
    printf("    Find the 68%% containment in bin %3i, center %f\n", idx, contain68); 
  }
  else
  {
    printf("attempt to divide by zero when setting integrals\n");
  }

  //Finally, deal with the cut efficiency
  int nEntries = pHist->GetEntries();

  if (numTotal < 1) {numTotal = pTree->GetEntries();}

  if (numTotal) {cutEffic = ((Double_t) nEntries) / ((Double_t) numTotal);}
  else          {cutEffic = 0.;}

  printf("    -- Total events:%i, events in hist:%i\n",numTotal,nEntries);

  printf("    Efficiency with this cut: %5.2f%%\n", 100.*cutEffic);

  return;
}


void PSFplot::fillHists(TTree* pTree, char* pFillVal, TCut cutVal)
{
  //Define the histogram for this variable
  pHist = new TH1F(pHistName->Data(),pHistTitle->Data(),nBins,lowBinEdge,hiBinEdge);

  //Fill the histogram
  pTree->Project(pHistName->Data(), pFillVal, cutVal);

  return;
}

//**********************************************************************
// Display the plot(s)
//
void PSFplot::plotHists(TCanvas* pCan, TPad* pPad) 
{
  //Make the first plot
  pPad->Draw();
  pPad->cd();
  pPad->SetLogy(0);
  pPad->GetFrame()->SetFillColor(10);
  pPad->SetGrid();   
  pHist->SetFillColor(9);
  pHist->SetLabelSize(0.07,"X");
  pHist->SetLabelSize(0.07,"Y");
  //****************************USER**********************************
  //range is set on each layer, this is a default
  //pHist[i]->GetXaxis()->SetRange(0,15);
  //****************************USER**********************************
  //SetMaximum sets the Y-axis max range upon # of events selected;
  // i.e., 50 if Go(10), 500 if Go(1000)
  //pHist[i]->SetMaximum(500);
  //    pHist[i]->SetXTitle("Nhits");
  //    pHist[i]->SetYTitle("Nevents");
  pHist->SetTitleOffset(1.3,"X");
  pHist->SetTitleOffset(1.3,"Y");
  pHist->SetTitleSize(0.09,"X");
  pHist->SetTitleSize(0.09,"Y");
  pHist->Draw();

  pCan->Update();

  //If a second plot to overlay, then do here
  if (pSum)
  {
    Double_t yAxisMax = pPad->GetUymax();
    axis = new TGaxis(pPad->GetUxmax(), pPad->GetUymin(),
                      pPad->GetUxmax(), pPad->GetUymax(),
                      0, 1.1, 510, "+L");
    axis->SetLineColor(2); //Red
    axis->SetLabelSize(0.07);
    axis->Draw();
    pSum->Scale(yAxisMax/1.1);
    pSum->SetFillColor(0);
    pSum->SetLineColor(2); //Red
    pSum->Draw("same");

    //Extract information for sizing the labelling of the plot
    TAxis*   pAxisX  = pHist->GetXaxis();
    Axis_t   loXEdge = pAxisX->GetXmin();
    Axis_t   hiXEdge = pAxisX->GetXmax();
    Double_t xDiff   = hiXEdge - loXEdge;

    loXEdge          = hiXEdge - 0.55 * xDiff;
    hiXEdge          = hiXEdge - 0.05 * xDiff;

    Double_t loYEdge = pHist->GetMinimum();
    Double_t hiYEdge = pHist->GetMaximum();
    Double_t yDiff   = hiYEdge - loYEdge;

    loYEdge          = loYEdge + 0.55 * yDiff;
    hiYEdge          = loYEdge + 0.20 * yDiff;
    
    //Now create the 68% label to output
    char labelText[30];

    sprintf(labelText,"68%% point: %6.4f",contain68);
    
    //Label the overlaid plot
    pLabel68 = new TPaveLabel(loXEdge,loYEdge,hiXEdge,hiYEdge,labelText);
    pLabel68->SetBorderSize(0);
    pLabel68->SetFillColor(19);
    pLabel68->SetTextAlign(32);
    pLabel68->SetTextSize(0.60);
    pLabel68->Draw();
    
    //Finally, create the 95% label to output
    sprintf(labelText,"95%% point: %6.4f",contain95);

    hiYEdge = loYEdge;
    loYEdge = loYEdge - 0.20 * yDiff;
    
    //Label the overlaid plot
    pLabel95 = new TPaveLabel(loXEdge,loYEdge,hiXEdge,hiYEdge,labelText);
    pLabel95->SetBorderSize(0);
    pLabel95->SetFillColor(19);
    pLabel95->SetTextAlign(32);
    pLabel95->SetTextSize(0.60);
    pLabel95->Draw();
  }

  pCan->cd();
}