/***************************************************************************
 *
 * $Id$
 *
 * Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
 ***************************************************************************
 *
 * Description: part of STAR HBT Framework: AliFemtoMaker package
 *   3D Bertsch-Pratt decomposition in the LCMS frame
 *
 ***************************************************************************
 *
 * $Log$
 * Revision 1.1.1.1  2007/04/25 15:38:41  panos
 * Importing the HBT code dir
 *
 * Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
 * First version on CVS
 *
 * Revision 1.7  2003/01/31 19:20:54  magestro
 * Cleared up simple compiler warnings on i386_linux24
 *
 * Revision 1.6  2002/06/07 22:51:38  lisa
 * Widely used AliFemtoBPLCMS3DCorrFctn class now accumulates UNcorrected denominator and has a WriteOutHistos method
 *
 * Revision 1.5  2001/05/23 00:19:04  lisa
 * Add in Smearing classes and methods needed for momentum resolution studies and correction
 *
 * Revision 1.4  2000/10/26 19:48:49  rcwells
 * Added functionality for Coulomb correction of <qInv> in 3D correltions
 *
 * Revision 1.3  2000/09/14 18:36:53  lisa
 * Added Qinv and ExitSep pair cuts and AliFemtoBPLCMS3DCorrFctn_SIM CorrFctn
 *
 * Revision 1.2  2000/08/23 19:43:43  lisa
 * added alternate normalization algorithm to 3d CorrFctns in case normal one fails
 *
 * Revision 1.1  2000/08/17 20:48:39  lisa
 * Adding correlationfunction in LCMS frame
 *
 *
 **************************************************************************/

#include "CorrFctn/AliFemtoBPLCMS3DCorrFctn.h"
//#include "Infrastructure/AliFemtoHisto.h"
#include <cstdio>

#ifdef __ROOT__ 
ClassImp(AliFemtoBPLCMS3DCorrFctn)
#endif

//____________________________
AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi)
  :
  fIDNumHisto(0),
  fIDDenHisto(0),
  fIDRatHisto(0),
  fSMNumHisto(0),
  fSMDenHisto(0),
  fSMRatHisto(0),
  fCorrectionHisto(0),
  fCorrCFHisto(0),
  fNumerator(0),
  fDenominator(0),
  fRatio(0),
  fQinvHisto(0),
  fLambda(0),
  fRout2(0),
  fRside2(0),
  fRlong2(0),
  fPairCut(0), 
  fQinvNormLo(0),
  fQinvNormHi(0),
  fNumRealsNorm(0),
  fNumMixedNorm(0)
{

  // set some stuff...
  fQinvNormLo = 0.15;
  fQinvNormHi = 0.18;
  fNumRealsNorm = 0;
  fNumMixedNorm = 0;
  //  fCorrection = 0;  // pointer to Coulomb Correction object

  fPairCut = 0; // added Sept2000 - CorrFctn-specific PairCut

  //  fSmearPair = 0; // no resolution correction unless user sets SmearPair

  // set up numerator
  char TitNum[100] = "Num";
  strcat(TitNum,title);
  fNumerator = new TH3D(TitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  // set up denominator
  char TitDen[100] = "Den";
  strcat(TitDen,title);
  fDenominator = new TH3D(TitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  // set up uncorrected denominator
  char TitDenUncoul[100] = "DenNoCoul";
  strcat(TitDenUncoul,title);
  //  fUncorrectedDenominator = new TH3D(TitDenUncoul,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  // set up ratio
  char TitRat[100] = "Rat";
  strcat(TitRat,title);
  fRatio = new TH3D(TitRat,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  // set up ave qInv
  char TitQinv[100] = "Qinv";
  strcat(TitQinv,title);
  fQinvHisto = new TH3D(TitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);

  // to enable error bar calculation...
  fNumerator->Sumw2();
  fDenominator->Sumw2();
  //  fUncorrectedDenominator->Sumw2();
  fRatio->Sumw2();

  // Following histos are for the momentum resolution correction
  // they are filled only if a AliFemtoSmear object is plugged in
  // here comes the "idea" numerator and denominator and ratio...
  char TitNumID[100] = "IDNum";
  strcat(TitNumID,title);
  fIDNumHisto = new TH3D(TitNumID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  char TitDenID[100] = "IDDen";
  strcat(TitDenID,title);
  fIDDenHisto = new TH3D(TitDenID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  char TitRatID[100] = "IDRat";
  strcat(TitRatID,title);
  fIDRatHisto = new TH3D(TitRatID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);

  fIDNumHisto->Sumw2();
  fIDDenHisto->Sumw2();
  fIDRatHisto->Sumw2();

  //
  // here comes the "smeared" numerator and denominator...
  char TitNumSM[100] = "SMNum";
  strcat(TitNumSM,title);
  fSMNumHisto = new TH3D(TitNumSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  char TitDenSM[100] = "SMDen";
  strcat(TitDenSM,title);
  fSMDenHisto = new TH3D(TitDenSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  char TitRatSM[100] = "SMRat";
  strcat(TitRatSM,title);
  fSMRatHisto = new TH3D(TitRatSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  //
  fSMNumHisto->Sumw2();
  fSMDenHisto->Sumw2();
  fSMRatHisto->Sumw2();
  //
  // here comes the correction factor (which is just ratio of ideal ratio to smeared ratio)
  char TitCorrection[100] = "CorrectionFactor";
  strcat(TitCorrection,title);
  fCorrectionHisto = new TH3D(TitCorrection,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);  
  fCorrectionHisto->Sumw2();
  // here comes the fully corrected correlation function
  char TitCorrCF[100] = "CorrectedCF";
  strcat(TitCorrCF,title);
  fCorrCFHisto = new TH3D(TitCorrCF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  fCorrCFHisto->Sumw2();

  // user can (and should) override these defaults...
  fLambda = 0.6;
  fRout2 = 6.0*6.0;
  fRside2 = 6.0*6.0;
  fRlong2 = 7.0*7.0;

}

AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(const AliFemtoBPLCMS3DCorrFctn& aCorrFctn) :
  fIDNumHisto(0),
  fIDDenHisto(0),
  fIDRatHisto(0),
  fSMNumHisto(0),
  fSMDenHisto(0),
  fSMRatHisto(0),
  fCorrectionHisto(0),
  fCorrCFHisto(0),
  fNumerator(0),
  fDenominator(0),
  fRatio(0),
  fQinvHisto(0),
  fLambda(0),
  fRout2(0),
  fRside2(0),
  fRlong2(0),
  fPairCut(0), 
  fQinvNormLo(0),
  fQinvNormHi(0),
  fNumRealsNorm(0),
  fNumMixedNorm(0)
{
  fIDNumHisto = aCorrFctn.fIDNumHisto;
  fIDDenHisto = aCorrFctn.fIDDenHisto;
  fIDRatHisto = aCorrFctn.fIDRatHisto;
  fSMNumHisto = aCorrFctn.fSMNumHisto;
  fSMDenHisto = aCorrFctn.fSMDenHisto;
  fSMRatHisto = aCorrFctn.fSMRatHisto;
  fCorrectionHisto = aCorrFctn.fCorrectionHisto;
  fCorrCFHisto = aCorrFctn.fCorrCFHisto;
  fNumerator = aCorrFctn.fNumerator;
  fDenominator = aCorrFctn.fDenominator;
  fRatio = aCorrFctn.fRatio;
  fQinvHisto = aCorrFctn.fQinvHisto;
  fLambda = aCorrFctn.fLambda;
  fRout2 = aCorrFctn.fRout2;
  fRside2 = aCorrFctn.fRside2;
  fRlong2 = aCorrFctn.fRlong2;
  fPairCut = aCorrFctn.fPairCut; 
  fQinvNormLo = aCorrFctn.fQinvNormLo;
  fQinvNormHi = aCorrFctn.fQinvNormHi;
  fNumRealsNorm = aCorrFctn.fNumRealsNorm;
  fNumMixedNorm = aCorrFctn.fNumMixedNorm;
}
//____________________________
AliFemtoBPLCMS3DCorrFctn::~AliFemtoBPLCMS3DCorrFctn(){
  delete fNumerator;
  delete fDenominator;
  delete fRatio;
  delete fQinvHisto;
  delete fIDNumHisto;
  delete fIDDenHisto;
  delete fIDRatHisto;
  delete fSMNumHisto;
  delete fSMDenHisto;
  delete fSMRatHisto;
  delete fCorrectionHisto;
  delete fCorrCFHisto;
}
//_________________________
AliFemtoBPLCMS3DCorrFctn& AliFemtoBPLCMS3DCorrFctn::operator=(const AliFemtoBPLCMS3DCorrFctn& aCorrFctn)
{
  if (this == &aCorrFctn)
    return *this;
  if (fIDNumHisto) delete fIDNumHisto;
  fIDNumHisto = new TH3D(*aCorrFctn.fIDNumHisto);
  if (fIDDenHisto) delete fIDDenHisto;
  fIDDenHisto = new TH3D(*aCorrFctn.fIDDenHisto);
  if (fIDRatHisto) delete fIDRatHisto;
  fIDRatHisto = new TH3D(*aCorrFctn.fIDRatHisto);
  if (fSMNumHisto) delete fSMNumHisto;
  fSMNumHisto = new TH3D(*aCorrFctn.fSMNumHisto);
  if (fSMDenHisto) delete fSMDenHisto;
  fSMDenHisto = new TH3D(*aCorrFctn.fSMDenHisto);
  if (fSMRatHisto) delete fSMRatHisto;
  fSMRatHisto = new TH3D(*aCorrFctn.fSMRatHisto);

  if (fCorrectionHisto) delete fCorrectionHisto;
  fCorrectionHisto = new TH3D(*aCorrFctn.fCorrectionHisto);
  if (fCorrCFHisto) delete fCorrCFHisto;
  fCorrCFHisto = new TH3D(*aCorrFctn.fCorrCFHisto);
  if (fNumerator) delete fNumerator;
  fNumerator = new TH3D(*aCorrFctn.fNumerator);
  if (fDenominator) delete fDenominator;
  fDenominator = new TH3D(*aCorrFctn.fDenominator);
  if (fRatio) delete fRatio;
  fRatio = new TH3D(*aCorrFctn.fRatio);
  if (fQinvHisto) delete fQinvHisto;
  fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);

  fLambda = aCorrFctn.fLambda;
  fRout2 = aCorrFctn.fRout2;
  fRside2 = aCorrFctn.fRside2;
  fRlong2 = aCorrFctn.fRlong2;
  fPairCut = aCorrFctn.fPairCut; 
  fQinvNormLo = aCorrFctn.fQinvNormLo;
  fQinvNormHi = aCorrFctn.fQinvNormHi;
  fNumRealsNorm = aCorrFctn.fNumRealsNorm;
  fNumMixedNorm = aCorrFctn.fNumMixedNorm;
  
  return *this;
}

//_________________________
void AliFemtoBPLCMS3DCorrFctn::WriteOutHistos(){

  fNumerator->Write();
  fDenominator->Write();
  //  fUncorrectedDenominator->Write();
  fRatio->Write();
  fQinvHisto->Write();

  /*
    if (fSmearPair){
    fIDNumHisto->Write();
    fIDDenHisto->Write();
    fIDRatHisto->Write();
    //
    fSMNumHisto->Write();
    fSMDenHisto->Write();
    fSMRatHisto->Write();
    //
    fCorrectionHisto->Write();
    fCorrCFHisto->Write();
    }
  */
}

//_________________________
void AliFemtoBPLCMS3DCorrFctn::Finish(){
  // here is where we should normalize, fit, etc...
  double NumFact,DenFact;
  if ((fNumRealsNorm !=0) && (fNumMixedNorm !=0)){
    NumFact = double(fNumRealsNorm);
    DenFact = double(fNumMixedNorm);
  }
  // can happen that the fNumRealsNorm and fNumMixedNorm = 0 if you do non-standard
  //   things like making a new CorrFctn and just setting the Numerator and Denominator
  //   from OTHER CorrFctns which you read in (like when doing parallel processing) 
  else{
    cout << "Warning! - no normalization constants defined - I do the best I can..." << endl;
    int nbins = fNumerator->GetNbinsX();
    int half_way = nbins/2;
    NumFact = fNumerator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
    DenFact = fDenominator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
  }

  fRatio->Divide(fNumerator,fDenominator,DenFact,NumFact);
  //  fQinvHisto->Divide(fUncorrectedDenominator);
  fQinvHisto->Divide(fDenominator);

  /*
  // now do all the resolution correction stuff..
  if (fSmearPair){  // but only do it if we have been working with a SmearPair
  fIDRatHisto->Divide(fIDNumHisto,fIDDenHisto);
  fSMRatHisto->Divide(fSMNumHisto,fSMDenHisto);
  fCorrectionHisto->Divide(fIDRatHisto,fSMRatHisto);
  fCorrCFHisto->Multiply(fRatio,fCorrectionHisto);
  }
  */

}

//____________________________
AliFemtoString AliFemtoBPLCMS3DCorrFctn::Report(){
  string stemp = "LCMS Frame Bertsch-Pratt 3D Correlation Function Report:\n";
  char ctemp[100];
  sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
  stemp += ctemp;
  sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
  stemp += ctemp;
  sprintf(ctemp,"Number of entries in ratio:\t%E\n",fRatio->GetEntries());
  stemp += ctemp;
  sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",fQinvNormLo,fQinvNormHi);
  stemp += ctemp;
  sprintf(ctemp,"Number of pairs in Normalization region was:\n");
  stemp += ctemp;
  sprintf(ctemp,"In numerator:\t%lu\t In denominator:\t%lu\n",fNumRealsNorm,fNumMixedNorm);
  stemp += ctemp;
  /*  if (fCorrection)
      {
      float radius = fCorrection->GetRadius();
      sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius);
      }
      else
      {
      sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n");
      }
      stemp += ctemp;
  */

  if (fPairCut){
    sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
    stemp += ctemp;
    stemp += fPairCut->Report();
  }
  else{
    sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
    stemp += ctemp;
  }

  //  
  AliFemtoString returnThis = stemp;
  return returnThis;
}
//____________________________
void AliFemtoBPLCMS3DCorrFctn::AddRealPair(const AliFemtoPair* pair){

  if (fPairCut){
    if (!(fPairCut->Pass(pair))) return;
  }

  double Qinv = fabs(pair->qInv());   // note - qInv() will be negative for identical pairs...
  if ((Qinv < fQinvNormHi) && (Qinv > fQinvNormLo)) fNumRealsNorm++;
  double qOut = fabs(pair->qOutCMS());
  double qSide = fabs(pair->qSideCMS());
  double qLong = fabs(pair->qLongCMS());

  fNumerator->Fill(qOut,qSide,qLong);
}
//____________________________
void AliFemtoBPLCMS3DCorrFctn::AddMixedPair(const AliFemtoPair* pair){

  if (fPairCut){
    if (!(fPairCut->Pass(pair))) return;
  }

  //  double CoulombWeight = (fCorrection ? fCorrection->CoulombCorrect(pair) : 1.0);
  double CoulombWeight = 1.0;

  double Qinv = fabs(pair->qInv());   // note - qInv() will be negative for identical pairs...
  if ((Qinv < fQinvNormHi) && (Qinv > fQinvNormLo)) fNumMixedNorm++;
  double qOut = fabs(pair->qOutCMS());
  double qSide = fabs(pair->qSideCMS());
  double qLong = fabs(pair->qLongCMS());

  fDenominator->Fill(qOut,qSide,qLong,CoulombWeight);
  //  fUncorrectedDenominator->Fill(qOut,qSide,qLong,1.0);
  fQinvHisto->Fill(qOut,qSide,qLong,Qinv);

  /*
  // now for the momentum resolution stuff...
  if (fSmearPair){
      double CorrWeight =  1.0 + 
      fLambda*exp((-qOut*qOut*fRout2 -qSide*qSide*fRside2 -qLong*qLong*fRlong2)/0.038936366329);
    CorrWeight *= CoulombWeight;  // impt.

    fIDNumHisto->Fill(qOut,qSide,qLong,CorrWeight);
    fIDDenHisto->Fill(qOut,qSide,qLong,CoulombWeight);

    fSmearPair->SetUnsmearedPair(pair);
    double qOut_prime = fabs(fSmearPair->SmearedPair().qOutCMS());
    double qSide_prime = fabs(fSmearPair->SmearedPair().qSideCMS());
    double qLong_prime = fabs(fSmearPair->SmearedPair().qLongCMS());

    fSMNumHisto->Fill(qOut_prime,qSide_prime,qLong_prime,CorrWeight);

    double SmearedCoulombWeight = ( fCorrection ? 
				    fCorrection->CoulombCorrect(&(fSmearPair->SmearedPair())) : 
				    1.0);

    fSMDenHisto->Fill(qOut_prime,qSide_prime,qLong_prime,SmearedCoulombWeight);
  }
  */
}