//////////////////////////////////////////////////////////////////////////
//                                                                      //
// GFDstarCount                                                         //
//
//   Author:      Gero Flucke
//   Date:        October 8th, 2003
//   last update: $Date: 2005/03/09 12:45:31 $
//   by:          $Author: flucke $
//                                                                      //
// Encapsulate Gero's routines for counting D*                          // 
//   (prior part of GFHistManip)                                        //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include "GFAnalysis/GFDstarCount.h"

#include <iostream>
#include <string.h> 

#include "TMath.h"
#include "TObjArray.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TAxis.h"
#include "GFUtils/GFHistArray.h"

#include "/afs/desy.de/user/f/flucke/h1/root/DstarDmFitter/DstarDmFitter.h"

using namespace std;

ClassImp(GFDstarCount)

//__________________________________________________________________
 TArrayD GFDstarCount::NDstarBackgrSubtracted(TH1* dstarHist, TH1* wcHist, 
					    Double_t dmLow, Double_t dmHigh, 
					    Double_t startBack, Double_t endBack)
{
  // N(D*) from scaling up wc due to factor between startBack and endBack
  // and subtracting background in region dmLow to dmHigh
  // 0: N(D*), 1: sigma(N), 2: signal to background ratio 3: scale for wcHist  

  if(!dstarHist || ! wcHist){
    cerr << "<GFDstarCount::NDstarBackgrSubtracted>: missing a hist!" << endl;
    return TArrayD(4);
  }
  const Int_t binStartBack = dstarHist->FindBin(startBack);
  const Int_t binEndBack = dstarHist->FindBin(endBack);
  const Double_t dstarBack = dstarHist->Integral(binStartBack, binEndBack);
  const Double_t wcBack = wcHist->Integral(binStartBack, binEndBack);

  const Double_t scale = wcBack ? dstarBack/wcBack : 0.;

  const Int_t binDmLow = dstarHist->FindBin(dmLow);
  const Int_t binDmHigh = dstarHist->FindBin(dmHigh);
  const Double_t numDstar = dstarHist->Integral(binDmLow, binDmHigh);
  const Double_t numWc = wcHist->Integral(binDmLow, binDmHigh);

  TArrayD result(4);
  result[0] = numDstar - numWc * scale;
  result[1] = TMath::Sqrt(numDstar + numWc * scale * 3.); // cf. logbook 18.4.2004
  result[2] = (numWc ? numDstar/numWc : 0.);
  result[3] = scale;

  return result;
}


//__________________________________________________________________
 TArrayD GFDstarCount::NDstarBackgrFuncSubtracted(TH1* dstarHist, TH1* wcHist, 
						Double_t dmLow, Double_t dmHigh, 
						Double_t startBack, Double_t endBack)
{
  // N(D*) from scaling up wc due to factor between startBack and endBack
  // and subtracting background in region dmLow to dmHigh (but result from fit!)
  // 0: N(D*), 1: sigma(N), 2: signal to background ratio 3: scale for wcHist

  if(!dstarHist || ! wcHist){
    cerr << "<GFDstarCount::NDstarBackgrFuncSubtracted>: missing a hist!" << endl;
    return TArrayD(4);
  }
  
  Int_t binDmLow = dstarHist->FindBin(dmLow);
  if(binDmLow != dstarHist->FindBin(dmLow+0.00001)){
    binDmLow = dstarHist->FindBin(dmLow+0.00001);
  }
  Int_t binDmHigh = dstarHist->FindBin(dmHigh);
  if(binDmHigh != dstarHist->FindBin(dmHigh-0.00001)){
    binDmHigh = dstarHist->FindBin(dmHigh-0.00001);
  }
  const Double_t numDstar = dstarHist->Integral(binDmLow, binDmHigh);

  DstarDmFitter fitter(wcHist);
  fitter.Fit(-1);
  // adjust for same edges of fit and func->Integral:
  dmLow = dstarHist->GetBinLowEdge(binDmLow);
  dmHigh = dstarHist->GetBinLowEdge(binDmHigh) + dstarHist->GetBinWidth(binDmHigh);
  const Double_t numWc = 
    fitter.GetFittedFunc()->Integral(dmLow, dmHigh) / fitter.GetHistBinSize();

  const Double_t scale = GetScale(dstarHist, wcHist, startBack, endBack);
//   wcHist->Scale(scale);
//cout << "GFDstarCount::NDstarBackgrFuncSubtracted: Scaling up background by "<< scale << endl;

  TArrayD result(4);
  result[0] = numDstar - numWc*scale;
  result[1] = TMath::Sqrt(numDstar + numWc * scale * 3.);//cf.logbook 18.4.2004 with fit ~ count
  result[2] = (numWc*scale ? numDstar / numWc*scale : 0.);
  result[3] = scale;

  return result;
}

//__________________________________________________________________
 TArrayD GFDstarCount::NDstarBackgrFuncSubtracted(TH1* dstarHist, TH1* wcHist, 
						Double_t dmLow, Double_t dmHigh) 
{
  // N(D*) from N(D*) in between dmLow and dmHigh
  // and subtracting background fit in region dmLow to dmHigh (adjusted to next bin border)
  // 0: N(D*), 1: sigma(N), 2: just == 1.
  // sigma is not correct,at least not if wcHist as been scaled before!

  if(!dstarHist || ! wcHist){
    cerr << "<GFDstarCount::NDstarBackgrFuncSubtracted>: missing a hist!" << endl;
    return TArrayD(4);
  }
  Int_t binDmLow = dstarHist->FindBin(dmLow);
  if(binDmLow != dstarHist->FindBin(dmLow+0.00001)){
    binDmLow = dstarHist->FindBin(dmLow+0.00001);
  }
  Int_t binDmHigh = dstarHist->FindBin(dmHigh);
  if(binDmHigh != dstarHist->FindBin(dmHigh-0.00001)){
    binDmHigh = dstarHist->FindBin(dmHigh-0.00001);
  }
  const Double_t numDstar = dstarHist->Integral(binDmLow, binDmHigh);

  DstarDmFitter fitter(wcHist);
  fitter.Fit(-1);
  // adjust for same edges of fit and func->Integral:
  dmLow = dstarHist->GetBinLowEdge(binDmLow);
  dmHigh = dstarHist->GetBinLowEdge(binDmHigh) + dstarHist->GetBinWidth(binDmHigh);
  const Double_t numWc = 
    fitter.GetFittedFunc()->Integral(dmLow, dmHigh) / fitter.GetHistBinSize();

  TArrayD result(4);
  result[0] = numDstar - numWc;
  result[1] = TMath::Sqrt(numDstar + numWc);
  result[2] = (numWc ? numDstar/numWc : 0.);
  result[3] = 1.;

  return result;
}


//___________________________________________________________
 Double_t GFDstarCount::GetScale(TH1* dstarHist, TH1* wcHist, 
			       Double_t startBack, Double_t endBack)
{
  // returns factor that has to be applied to wcHist to get same number of entries 
  // as dstarHist from startback to endBack
  // assuming same binning of both hists!!!!!
  // 0. if wcHist is empty in that range

  const Int_t binStartBack = dstarHist->FindBin(startBack);
  const Int_t binEndBack = dstarHist->FindBin(endBack);
  const Double_t dstarBack = dstarHist->Integral(binStartBack, binEndBack);
  const Double_t wcBack = wcHist->Integral(binStartBack, binEndBack);
  
  return wcBack ? dstarBack/wcBack : 0.;
}