[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / CorrFctn / AliFemtoBPLCMS3DCorrFctn.cxx

/***************************************************************************
 *
 * $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);
  }
  */
}
Commit	Line	Data
67427ff7	1	/***************************************************************************
	2	*
	3	* $Id$
	4	*
	5	* Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
	6	***************************************************************************
	7	*
	8	* Description: part of STAR HBT Framework: AliFemtoMaker package
	9	* 3D Bertsch-Pratt decomposition in the LCMS frame
	10	*
	11	***************************************************************************
	12	*
	13	* $Log$
0215f606	14	* Revision 1.1.1.1 2007/04/25 15:38:41 panos
	15	* Importing the HBT code dir
	16	*
67427ff7	17	* Revision 1.1.1.1 2007/03/07 10:14:49 mchojnacki
	18	* First version on CVS
	19	*
	20	* Revision 1.7 2003/01/31 19:20:54 magestro
	21	* Cleared up simple compiler warnings on i386_linux24
	22	*
	23	* Revision 1.6 2002/06/07 22:51:38 lisa
	24	* Widely used AliFemtoBPLCMS3DCorrFctn class now accumulates UNcorrected denominator and has a WriteOutHistos method
	25	*
	26	* Revision 1.5 2001/05/23 00:19:04 lisa
	27	* Add in Smearing classes and methods needed for momentum resolution studies and correction
	28	*
	29	* Revision 1.4 2000/10/26 19:48:49 rcwells
	30	* Added functionality for Coulomb correction of <qInv> in 3D correltions
	31	*
	32	* Revision 1.3 2000/09/14 18:36:53 lisa
	33	* Added Qinv and ExitSep pair cuts and AliFemtoBPLCMS3DCorrFctn_SIM CorrFctn
	34	*
	35	* Revision 1.2 2000/08/23 19:43:43 lisa
	36	* added alternate normalization algorithm to 3d CorrFctns in case normal one fails
	37	*
	38	* Revision 1.1 2000/08/17 20:48:39 lisa
	39	* Adding correlationfunction in LCMS frame
	40	*
	41	*
	42	**************************************************************************/
	43
	44	#include "CorrFctn/AliFemtoBPLCMS3DCorrFctn.h"
	45	//#include "Infrastructure/AliFemtoHisto.h"
	46	#include <cstdio>
	47
	48	#ifdef __ROOT__
	49	ClassImp(AliFemtoBPLCMS3DCorrFctn)
	50	#endif
	51
	52	//____________________________
0215f606	53	AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi)
	54	:
	55	fIDNumHisto(0),
	56	fIDDenHisto(0),
	57	fIDRatHisto(0),
	58	fSMNumHisto(0),
	59	fSMDenHisto(0),
	60	fSMRatHisto(0),
	61	fCorrectionHisto(0),
	62	fCorrCFHisto(0),
	63	fNumerator(0),
	64	fDenominator(0),
	65	fRatio(0),
	66	fQinvHisto(0),
	67	fLambda(0),
	68	fRout2(0),
	69	fRside2(0),
	70	fRlong2(0),
	71	fPairCut(0),
	72	fQinvNormLo(0),
	73	fQinvNormHi(0),
	74	fNumRealsNorm(0),
	75	fNumMixedNorm(0)
	76	{
67427ff7	77
	78	// set some stuff...
	79	fQinvNormLo = 0.15;
	80	fQinvNormHi = 0.18;
	81	fNumRealsNorm = 0;
	82	fNumMixedNorm = 0;
	83	// fCorrection = 0; // pointer to Coulomb Correction object
	84
	85	fPairCut = 0; // added Sept2000 - CorrFctn-specific PairCut
	86
	87	// fSmearPair = 0; // no resolution correction unless user sets SmearPair
	88
	89	// set up numerator
	90	char TitNum[100] = "Num";
	91	strcat(TitNum,title);
	92	fNumerator = new TH3D(TitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	93	// set up denominator
	94	char TitDen[100] = "Den";
	95	strcat(TitDen,title);
	96	fDenominator = new TH3D(TitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	97	// set up uncorrected denominator
	98	char TitDenUncoul[100] = "DenNoCoul";
	99	strcat(TitDenUncoul,title);
	100	// fUncorrectedDenominator = new TH3D(TitDenUncoul,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	101	// set up ratio
	102	char TitRat[100] = "Rat";
	103	strcat(TitRat,title);
	104	fRatio = new TH3D(TitRat,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	105	// set up ave qInv
	106	char TitQinv[100] = "Qinv";
	107	strcat(TitQinv,title);
	108	fQinvHisto = new TH3D(TitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	109
	110	// to enable error bar calculation...
	111	fNumerator->Sumw2();
	112	fDenominator->Sumw2();
	113	// fUncorrectedDenominator->Sumw2();
	114	fRatio->Sumw2();
	115
	116	// Following histos are for the momentum resolution correction
	117	// they are filled only if a AliFemtoSmear object is plugged in
	118	// here comes the "idea" numerator and denominator and ratio...
	119	char TitNumID[100] = "IDNum";
	120	strcat(TitNumID,title);
	121	fIDNumHisto = new TH3D(TitNumID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	122	char TitDenID[100] = "IDDen";
	123	strcat(TitDenID,title);
	124	fIDDenHisto = new TH3D(TitDenID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	125	char TitRatID[100] = "IDRat";
	126	strcat(TitRatID,title);
	127	fIDRatHisto = new TH3D(TitRatID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	128
	129	fIDNumHisto->Sumw2();
	130	fIDDenHisto->Sumw2();
	131	fIDRatHisto->Sumw2();
	132
	133	//
	134	// here comes the "smeared" numerator and denominator...
	135	char TitNumSM[100] = "SMNum";
	136	strcat(TitNumSM,title);
	137	fSMNumHisto = new TH3D(TitNumSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	138	char TitDenSM[100] = "SMDen";
	139	strcat(TitDenSM,title);
	140	fSMDenHisto = new TH3D(TitDenSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
141	char TitRatSM[100] = "SMRat";
142	strcat(TitRatSM,title);
143	fSMRatHisto = new TH3D(TitRatSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
144	//
145	fSMNumHisto->Sumw2();
146	fSMDenHisto->Sumw2();
147	fSMRatHisto->Sumw2();
148	//
149	// here comes the correction factor (which is just ratio of ideal ratio to smeared ratio)
150	char TitCorrection[100] = "CorrectionFactor";
151	strcat(TitCorrection,title);
152	fCorrectionHisto = new TH3D(TitCorrection,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
153	fCorrectionHisto->Sumw2();
154	// here comes the fully corrected correlation function
155	char TitCorrCF[100] = "CorrectedCF";
156	strcat(TitCorrCF,title);
157	fCorrCFHisto = new TH3D(TitCorrCF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
158	fCorrCFHisto->Sumw2();
159
160	// user can (and should) override these defaults...
161	fLambda = 0.6;
162	fRout2 = 6.0*6.0;
163	fRside2 = 6.0*6.0;
164	fRlong2 = 7.0*7.0;
165
166	}
167
0215f606	168	AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(const AliFemtoBPLCMS3DCorrFctn& aCorrFctn) :
	169	fIDNumHisto(0),
	170	fIDDenHisto(0),
	171	fIDRatHisto(0),
	172	fSMNumHisto(0),
	173	fSMDenHisto(0),
	174	fSMRatHisto(0),
	175	fCorrectionHisto(0),
	176	fCorrCFHisto(0),
	177	fNumerator(0),
	178	fDenominator(0),
	179	fRatio(0),
	180	fQinvHisto(0),
	181	fLambda(0),
	182	fRout2(0),
	183	fRside2(0),
	184	fRlong2(0),
	185	fPairCut(0),
	186	fQinvNormLo(0),
	187	fQinvNormHi(0),
	188	fNumRealsNorm(0),
	189	fNumMixedNorm(0)
	190	{
	191	fIDNumHisto = aCorrFctn.fIDNumHisto;
	192	fIDDenHisto = aCorrFctn.fIDDenHisto;
	193	fIDRatHisto = aCorrFctn.fIDRatHisto;
	194	fSMNumHisto = aCorrFctn.fSMNumHisto;
	195	fSMDenHisto = aCorrFctn.fSMDenHisto;
	196	fSMRatHisto = aCorrFctn.fSMRatHisto;
	197	fCorrectionHisto = aCorrFctn.fCorrectionHisto;
	198	fCorrCFHisto = aCorrFctn.fCorrCFHisto;
	199	fNumerator = aCorrFctn.fNumerator;
	200	fDenominator = aCorrFctn.fDenominator;
	201	fRatio = aCorrFctn.fRatio;
	202	fQinvHisto = aCorrFctn.fQinvHisto;
	203	fLambda = aCorrFctn.fLambda;
	204	fRout2 = aCorrFctn.fRout2;
	205	fRside2 = aCorrFctn.fRside2;
	206	fRlong2 = aCorrFctn.fRlong2;
	207	fPairCut = aCorrFctn.fPairCut;
	208	fQinvNormLo = aCorrFctn.fQinvNormLo;
	209	fQinvNormHi = aCorrFctn.fQinvNormHi;
	210	fNumRealsNorm = aCorrFctn.fNumRealsNorm;
	211	fNumMixedNorm = aCorrFctn.fNumMixedNorm;
	212	}
67427ff7	213	//____________________________
	214	AliFemtoBPLCMS3DCorrFctn::~AliFemtoBPLCMS3DCorrFctn(){
	215	delete fNumerator;
	216	delete fDenominator;
	217	delete fRatio;
	218	delete fQinvHisto;
	219	delete fIDNumHisto;
	220	delete fIDDenHisto;
	221	delete fIDRatHisto;
	222	delete fSMNumHisto;
	223	delete fSMDenHisto;
	224	delete fSMRatHisto;
	225	delete fCorrectionHisto;
	226	delete fCorrCFHisto;
	227	}
0215f606	228	//_________________________
	229	AliFemtoBPLCMS3DCorrFctn& AliFemtoBPLCMS3DCorrFctn::operator=(const AliFemtoBPLCMS3DCorrFctn& aCorrFctn)
	230	{
	231	if (this == &aCorrFctn)
	232	return *this;
	233	if (fIDNumHisto) delete fIDNumHisto;
	234	fIDNumHisto = new TH3D(*aCorrFctn.fIDNumHisto);
	235	if (fIDDenHisto) delete fIDDenHisto;
	236	fIDDenHisto = new TH3D(*aCorrFctn.fIDDenHisto);
	237	if (fIDRatHisto) delete fIDRatHisto;
	238	fIDRatHisto = new TH3D(*aCorrFctn.fIDRatHisto);
	239	if (fSMNumHisto) delete fSMNumHisto;
	240	fSMNumHisto = new TH3D(*aCorrFctn.fSMNumHisto);
	241	if (fSMDenHisto) delete fSMDenHisto;
	242	fSMDenHisto = new TH3D(*aCorrFctn.fSMDenHisto);
	243	if (fSMRatHisto) delete fSMRatHisto;
	244	fSMRatHisto = new TH3D(*aCorrFctn.fSMRatHisto);
	245
	246	if (fCorrectionHisto) delete fCorrectionHisto;
	247	fCorrectionHisto = new TH3D(*aCorrFctn.fCorrectionHisto);
	248	if (fCorrCFHisto) delete fCorrCFHisto;
	249	fCorrCFHisto = new TH3D(*aCorrFctn.fCorrCFHisto);
	250	if (fNumerator) delete fNumerator;
	251	fNumerator = new TH3D(*aCorrFctn.fNumerator);
	252	if (fDenominator) delete fDenominator;
	253	fDenominator = new TH3D(*aCorrFctn.fDenominator);
	254	if (fRatio) delete fRatio;
	255	fRatio = new TH3D(*aCorrFctn.fRatio);
	256	if (fQinvHisto) delete fQinvHisto;
	257	fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
	258
	259	fLambda = aCorrFctn.fLambda;
	260	fRout2 = aCorrFctn.fRout2;
	261	fRside2 = aCorrFctn.fRside2;
	262	fRlong2 = aCorrFctn.fRlong2;
	263	fPairCut = aCorrFctn.fPairCut;
	264	fQinvNormLo = aCorrFctn.fQinvNormLo;
	265	fQinvNormHi = aCorrFctn.fQinvNormHi;
	266	fNumRealsNorm = aCorrFctn.fNumRealsNorm;
	267	fNumMixedNorm = aCorrFctn.fNumMixedNorm;
	268
	269	return *this;
	270	}
67427ff7	271
	272	//_________________________
	273	void AliFemtoBPLCMS3DCorrFctn::WriteOutHistos(){
	274
	275	fNumerator->Write();
	276	fDenominator->Write();
	277	// fUncorrectedDenominator->Write();
	278	fRatio->Write();
	279	fQinvHisto->Write();
	280
	281	/*
	282	if (fSmearPair){
	283	fIDNumHisto->Write();
	284	fIDDenHisto->Write();
	285	fIDRatHisto->Write();
	286	//
	287	fSMNumHisto->Write();
	288	fSMDenHisto->Write();
	289	fSMRatHisto->Write();
	290	//
	291	fCorrectionHisto->Write();
	292	fCorrCFHisto->Write();
	293	}
	294	*/
	295	}
	296
	297	//_________________________
	298	void AliFemtoBPLCMS3DCorrFctn::Finish(){
	299	// here is where we should normalize, fit, etc...
	300	double NumFact,DenFact;
	301	if ((fNumRealsNorm !=0) && (fNumMixedNorm !=0)){
	302	NumFact = double(fNumRealsNorm);
	303	DenFact = double(fNumMixedNorm);
	304	}
	305	// can happen that the fNumRealsNorm and fNumMixedNorm = 0 if you do non-standard
	306	// things like making a new CorrFctn and just setting the Numerator and Denominator
	307	// from OTHER CorrFctns which you read in (like when doing parallel processing)
	308	else{
	309	cout << "Warning! - no normalization constants defined - I do the best I can..." << endl;
	310	int nbins = fNumerator->GetNbinsX();
	311	int half_way = nbins/2;
	312	NumFact = fNumerator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
	313	DenFact = fDenominator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
	314	}
	315
	316	fRatio->Divide(fNumerator,fDenominator,DenFact,NumFact);
	317	// fQinvHisto->Divide(fUncorrectedDenominator);
	318	fQinvHisto->Divide(fDenominator);
	319
	320	/*
	321	// now do all the resolution correction stuff..
	322	if (fSmearPair){ // but only do it if we have been working with a SmearPair
	323	fIDRatHisto->Divide(fIDNumHisto,fIDDenHisto);
	324	fSMRatHisto->Divide(fSMNumHisto,fSMDenHisto);
	325	fCorrectionHisto->Divide(fIDRatHisto,fSMRatHisto);
	326	fCorrCFHisto->Multiply(fRatio,fCorrectionHisto);
	327	}
	328	*/
	329
	330	}
	331
	332	//____________________________
	333	AliFemtoString AliFemtoBPLCMS3DCorrFctn::Report(){
	334	string stemp = "LCMS Frame Bertsch-Pratt 3D Correlation Function Report:\n";
335	char ctemp[100];
336	sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
337	stemp += ctemp;
338	sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
339	stemp += ctemp;
340	sprintf(ctemp,"Number of entries in ratio:\t%E\n",fRatio->GetEntries());
341	stemp += ctemp;
342	sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",fQinvNormLo,fQinvNormHi);
343	stemp += ctemp;
344	sprintf(ctemp,"Number of pairs in Normalization region was:\n");
345	stemp += ctemp;
346	sprintf(ctemp,"In numerator:\t%lu\t In denominator:\t%lu\n",fNumRealsNorm,fNumMixedNorm);
347	stemp += ctemp;
348	/* if (fCorrection)
349	{
350	float radius = fCorrection->GetRadius();
351	sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius);
352	}
353	else
354	{
355	sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n");
356	}
357	stemp += ctemp;
358	*/
359
360	if (fPairCut){
361	sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
362	stemp += ctemp;
363	stemp += fPairCut->Report();
364	}
365	else{
366	sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
367	stemp += ctemp;
368	}
369
370	//
371	AliFemtoString returnThis = stemp;
372	return returnThis;
373	}
374	//____________________________
375	void AliFemtoBPLCMS3DCorrFctn::AddRealPair(const AliFemtoPair* pair){
376
377	if (fPairCut){
378	if (!(fPairCut->Pass(pair))) return;
379	}
380
381	double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
382	if ((Qinv < fQinvNormHi) && (Qinv > fQinvNormLo)) fNumRealsNorm++;
383	double qOut = fabs(pair->qOutCMS());
384	double qSide = fabs(pair->qSideCMS());
385	double qLong = fabs(pair->qLongCMS());
386
387	fNumerator->Fill(qOut,qSide,qLong);
388	}
389	//____________________________
390	void AliFemtoBPLCMS3DCorrFctn::AddMixedPair(const AliFemtoPair* pair){
391
392	if (fPairCut){
393	if (!(fPairCut->Pass(pair))) return;
394	}
395
396	// double CoulombWeight = (fCorrection ? fCorrection->CoulombCorrect(pair) : 1.0);
397	double CoulombWeight = 1.0;
398
399	double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
400	if ((Qinv < fQinvNormHi) && (Qinv > fQinvNormLo)) fNumMixedNorm++;
401	double qOut = fabs(pair->qOutCMS());
402	double qSide = fabs(pair->qSideCMS());
403	double qLong = fabs(pair->qLongCMS());
404
405	fDenominator->Fill(qOut,qSide,qLong,CoulombWeight);
406	// fUncorrectedDenominator->Fill(qOut,qSide,qLong,1.0);
407	fQinvHisto->Fill(qOut,qSide,qLong,Qinv);
408
409	/*
410	// now for the momentum resolution stuff...
411	if (fSmearPair){
412	double CorrWeight = 1.0 +
413	fLambdaexp((-qOutqOutfRout2 -qSideqSidefRside2 -qLongqLong*fRlong2)/0.038936366329);
414	CorrWeight *= CoulombWeight; // impt.
415
416	fIDNumHisto->Fill(qOut,qSide,qLong,CorrWeight);
417	fIDDenHisto->Fill(qOut,qSide,qLong,CoulombWeight);
418
419	fSmearPair->SetUnsmearedPair(pair);
420	double qOut_prime = fabs(fSmearPair->SmearedPair().qOutCMS());
421	double qSide_prime = fabs(fSmearPair->SmearedPair().qSideCMS());
422	double qLong_prime = fabs(fSmearPair->SmearedPair().qLongCMS());
423
424	fSMNumHisto->Fill(qOut_prime,qSide_prime,qLong_prime,CorrWeight);
425
426	double SmearedCoulombWeight = ( fCorrection ?
427	fCorrection->CoulombCorrect(&(fSmearPair->SmearedPair())) :
428	1.0);
429
430	fSMDenHisto->Fill(qOut_prime,qSide_prime,qLong_prime,SmearedCoulombWeight);
431	}
432	*/
433	}
434
435