////////////////////////////////////////////////////////////
// Class H1JetProfile2D 
//
// Author     : Ger Flucke
// Created    : 2003/09/23
// Last update: $Date: 2004/11/24 08:45:06 $
//          by: $Author: flucke $
// 
// Fill Jet eta/phi profile histograms with an additional quantity on the y-axis.
// The idea is to make jet profiles where the bin content is determined by e.g.
// a fit in the y-variable. The fit has to be dne by the user.
//
////////////////////////////////////////////////////////////
//
//
// TH2 *histeta = new TH2F("etaprofile","etaprofile",20,-5,5, nBinsy, startY, endY);
// TH2 *histphi = new TH2F("phiprofile","phiprofile",20,-TMath::Pi(),TMath::Pi(), 
//                         nBinsy, startY, endY);
//
// H1PartInclJetArrayPtr jets;
// H1PartCandArrayPtr    parts;
// H1PartDstarArrayPtr   dstars;
//
// H1JetProfile2D profiler;
//
// profiler.SetEtaHistogram(histeta);
// profiler.SetPhiHistogram(histphi);
// 
// H1Tree::Instance()->SelectHat("NumDstar > 0 && NumInclKtJets > 0");
// while(H1Tree::Instance()->Next()){
//     for(Int_t i = 0; i < dstars.GetEntries(); ++i){
//         Double_t dm = dstars[i]->GetDm();
//         TObjArray *list = parts.GetArray(); // works even with any list of H1Part!
//         profiler.FillProfile(jets[0], list, dm);
//     }
// }
//
// profiler.NormaliseHistograms();
//
// TCanvas *profile_eta = new TCanvas;
// histeta->Draw();
// TCanvas *profile_phi = new TCanvas;
// histphi->Draw();
// 
// (Better project the hists along y-axis for each delta eta/phi bin, do a D*-fit 
//  and store the results in a new 1-D profile hist...)

#include <TH2.h>
#include <TCollection.h>

#include "GFData/H1JetProfile2D.h"
#include "H1Mods/H1PartJet.h"

#include "GFUtils/GFHistManip.h"

#include <iostream>
using namespace std;

ClassImp(H1JetProfile2D)

//______________________________________________________
 H1JetProfile2D::H1JetProfile2D(Bool_t relPt, Float_t distEta, Float_t distPhi)
  : fHistEta(NULL), fHistPhi(NULL), fHistNumJets(NULL), fFlagRelPt(relPt)
{
  // both distances set via SetDistances(...), set fFlagRelPt (cf. FillProfile(...))
  this->SetDistances(distEta, distPhi);
}

//______________________________________________________
 H1JetProfile2D::~H1JetProfile2D()
{
  delete fHistNumJets;
}

//______________________________________________________
 void H1JetProfile2D::SetDistances(Float_t distEta, Float_t distPhi)
{
  // distances below 0. are interpreted as default, i.e. 1.
  fDistEta = (distEta > 0. ? distEta : 1.);
  fDistPhi = (distPhi > 0. ? distPhi : 1.);
}

//______________________________________________________
 void H1JetProfile2D::SetHistsEtaPhi(TH2* histEta, TH2* histPhi)
{
  // setting hists and creating fHistNumJets with y-bins
  // the name of the latter is taken as the name of the given hists where the first occurence
  // of 'Eta', 'eta', 'Phi' or 'phi' will be replaced by 'Norm'
  // setting flag fFlagAbsPhi according to histPhi 

  fHistPhi = histPhi; 
  fHistEta = histEta;
  delete fHistNumJets;
  fHistNumJets = NULL;

  TH2* hRef = NULL;
  if(!fHistPhi && fHistEta) hRef = fHistEta;
  else if(fHistPhi && !fHistEta) hRef = fHistPhi;
  else if(fHistPhi && fHistEta){
    if(fHistPhi->GetNbinsY() != fHistEta->GetNbinsY() 
       || fHistPhi->GetYaxis()->GetXmax() != fHistEta->GetYaxis()->GetXmax()
       || fHistPhi->GetYaxis()->GetXmin() != fHistEta->GetYaxis()->GetXmin()){
      this->Warning("SetHistsEtaPhi", "y-bins differ, take phi hist to create fHistNumJets");
    }
//     hRef = fHistPhi;
    hRef = fHistEta;
  }

  if(hRef) {
    fHistNumJets = GFHistManip::CreateHistYBins(hRef, "y-axis");
    if(hRef->GetSumw2N()) fHistNumJets->Sumw2();
    TString name(hRef->GetName());
    // searching for first occurence of 'Eta', 'eta', 'Phi' or 'phi' to replace it by 'Norm'
    const char* phiToReplace = "";
    Ssiz_t firstPhi = name.Index("Phi");
    if(firstPhi == kNPOS) {
      firstPhi = name.Index("phi");
      if(firstPhi != kNPOS) phiToReplace = "phi";
    } else {
      phiToReplace = "Phi";
    }
    const char* etaToReplace = "";
    Ssiz_t firstEta = name.Index("Eta");
    if(firstEta == kNPOS) {
      firstEta = name.Index("eta");
      if(firstEta != kNPOS) etaToReplace = "eta";
    } else {
      etaToReplace = "Eta";
    }

    if(firstPhi != kNPOS && (firstEta == kNPOS || firstPhi < firstEta)){
      name.Replace(firstPhi, strlen(phiToReplace), "Norm");
    } else if(firstEta != kNPOS){
      name.Replace(firstEta, strlen(etaToReplace), "Norm");
    } else {
      name += "Norm";
    }
    fHistNumJets->SetName(name);
  }

  fFlagAbsPhi = (fHistPhi && fHistPhi->GetXaxis()->GetXmin() >= 0.);
}

//______________________________________________________
 Bool_t H1JetProfile2D::FillProfile(const H1PartJet *jet, const TCollection *parts, 
				 Double_t yVar, Double_t weight)
{
  // Extracts eta and phi of the jet, then loops over all particles i in the array 
  // and for
  //
  // - delta phi(i,jet) < fDistPhi calls 
  //                         fHistEta->Fill(delta eta(i,jet), yVar, pt*weight)
  // - delta eta(i,jet) < fDistEta calls 
  //                         fHistPhi->Fill(delta phi(i,jet), yVar, pt*weight)
  //
  //     where pt = (p_t,i/p_t,jet)  if fFlagRelPt = kTRUE (default)
  //                 p_t,i           otherwise.
  //
  // False if something went wrong.

  if(!jet || !parts) return kFALSE;
  if(fHistEta == NULL && fHistPhi == NULL){
    this->Warning("FillProfile","Pointer to histograms are zero! No calculation!");
    return kFALSE;
  }
  const Double_t jetPt  = jet->GetPt();
  if(fFlagRelPt && jetPt == 0.) {
    this->Warning("FillProfile","Cannot calculate relative pt for pt(jet)=0.!");
    return kFALSE;
  }
  const Double_t jetEta = jet->GetEta();
  const Double_t jetPhi = jet->GetPhi();

  // increase the number of jets that are in the histograms
  // is double rather than int because of the weight that is necessary
  // for proper normalisation
  fHistNumJets->Fill(yVar, weight);
    
  TIter next(parts);

  // loop over all particles
  while (H1Part *part = static_cast<H1Part*>(next())) {
    const Double_t partEta = part->GetEta();
    const Double_t partPhi = part->GetPhi();
    const Double_t partRelPt = (fFlagRelPt ? part->GetPt()/jetPt : part->GetPt());

    // this routine returnes the correct phi in [-pi,+pi]
    const Double_t deltaPhi = H1PartJet::GetRealPhi(partPhi - jetPhi);
    const Double_t absDeltaPhi = TMath::Abs(deltaPhi);
    const Double_t deltaEta = partEta - jetEta;

    if(TMath::Abs(deltaEta) < fDistEta){
      if(fHistPhi) fHistPhi->Fill((fFlagAbsPhi?absDeltaPhi:deltaPhi), yVar, partRelPt*weight);
    }
    if(absDeltaPhi < fDistPhi){
      if(fHistEta) fHistEta->Fill(deltaEta, yVar, partRelPt*weight);
    }
  }

  return kTRUE;
}

//______________________________________________________
 Double_t H1JetProfile2D::GetNumJetsInHistograms() const
{
  if(fHistNumJets) return fHistNumJets->Integral();
  else return 0.;
}

//______________________________________________________
 Bool_t H1JetProfile2D::NormaliseHistograms(Double_t n)
{
  // Scales the histograms by 1/n where n should the number of jets that entered
  // the histograms, best defined by the same fit that is applied to all bins of the 
  // profile hists.
  // If n is not given or 0, just devide through the (weighted) number of calls to 
  // FillProfile()
  //
  // Calling this routine more than once has no effect. 
  // True if called the first time and all OK, false otherwise.

  if(this->GetUniqueID() == 0){ // NormaliseHistograms() not called yet
    if(fHistEta == NULL && fHistPhi == NULL){
      this->Warning("NormaliseHistograms","Pointer to histograms are zero! Do nothing");
      return kFALSE;
    }
    Double_t scale = 0.;
    if(n) scale = 1./n;
    else if(fHistNumJets && fHistNumJets->Integral()) scale = 1./fHistNumJets->Integral();
//     else if(fNumJets) scale = 1./fNumJets;
    else {
      this->Warning("NormaliseHistograms","No jets filled? Do nothing");
      return kFALSE;
    } 
    this->SetUniqueID(1);
    if(fHistEta) fHistEta->Scale(scale);
    if(fHistPhi) fHistPhi->Scale(scale);

    return kTRUE;
  }

  return kFALSE;
}

//______________________________________________________
 void H1JetProfile2D::Print(Option_t *opt) const
{
  this->TObject::Print(opt);
  if(fHistEta) fHistEta->Print(opt);
  if(fHistPhi) fHistPhi->Print(opt);
}