[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemtoUser / AliFemtoModelBPLCMSCorrFctnKK.cxx

////////////////////////////////////////////////////////////////////////////////
///                                                                          ///
/// AliFemtoModelBPLCMSCorrFctnKK - the class for correlation function which   ///
/// uses the model framework and weight generation and calculated the 3D     ///
/// correlation function in the Bertsh-Pratt LCMS system                     ///
/// Authors: Adam Kisiel, kisiel@mps.ohio-state.edu                          ///
///                                                                          ///
////////////////////////////////////////////////////////////////////////////////
#include "AliFemtoModelBPLCMSCorrFctnKK.h"
#include "AliFemtoPair.h"
#include "AliFemtoModelManager.h"
#include "AliFemtoKTPairCut.h"
#include "AliFemtoAnalysisReactionPlane.h"
#include <cstdio>

#ifdef __ROOT__
ClassImp(AliFemtoModelBPLCMSCorrFctnKK)
#endif

//____________________________
AliFemtoModelBPLCMSCorrFctnKK::AliFemtoModelBPLCMSCorrFctnKK(char* title, const int& nbins, const float& QLo, const float& QHi)
  :
  AliFemtoModelCorrFctn(title, nbins, QLo, QHi),
  fNumerator3DTrue(0),
  fNumerator3DFake(0),
  fDenominator3D(0),
  fQinvHisto(0),
  fPairCut(0),
  fUseRPSelection(0),
  fNumerator3DTrueIdeal(0),
  fNumerator3DFakeIdeal(0),
  fDenominator3DIdeal(0)
{
  // set up true numerator
  char tTitNumT[101] = "Num3DTrue";
  strncat(tTitNumT,title, 101);
  fNumerator3DTrue = new TH3D(tTitNumT,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  // set up fake numerator
  char tTitNumF[101] = "Num3DFake";
  strncat(tTitNumF,title, 100);
  fNumerator3DFake = new TH3D(tTitNumF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  // set up denominator
  char tTitDen[101] = "Den3D";
  strncat(tTitDen,title, 100);
  fDenominator3D = new TH3D(tTitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  // set up ave qInv
  char tTitQinv[101] = "Qinv";
  strncat(tTitQinv,title, 100);
  fQinvHisto = new TH3D(tTitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);

 // set up true Ideal numerator
  char tTitNumTI[101] = "Num3DTrueIdeal";
  strncat(tTitNumTI,title, 100);
  fNumerator3DTrueIdeal = new TH3D(tTitNumTI,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);

 // set up fake Ideal numerator
  char tTitNumFI[101] = "Num3DFakeIdeal";
  strncat(tTitNumFI,title, 100);
  fNumerator3DFakeIdeal = new TH3D(tTitNumFI,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);

// set up denominator
  char tTitDenI[101] = "Den3DIdeal";
  strncat(tTitDenI,title, 100);
  fDenominator3DIdeal = new TH3D(tTitDenI,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);


  // to enable error bar calculation...
  fNumerator3DTrue->Sumw2();
  fNumerator3DFake->Sumw2();
  fDenominator3D->Sumw2();
  fNumerator3DTrueIdeal->Sumw2();
  fNumerator3DFakeIdeal->Sumw2();
  fDenominator3DIdeal->Sumw2();
}

AliFemtoModelBPLCMSCorrFctnKK::AliFemtoModelBPLCMSCorrFctnKK(const AliFemtoModelBPLCMSCorrFctnKK& aCorrFctn) :
  AliFemtoModelCorrFctn(aCorrFctn),
  fNumerator3DTrue(0),
  fNumerator3DFake(0),
  fDenominator3D(0),
  fQinvHisto(0),
  fPairCut(0),
  fUseRPSelection(0),
  fNumerator3DTrueIdeal(0),
  fNumerator3DFakeIdeal(0),
  fDenominator3DIdeal(0)
{
  // Copy constructor
  fNumerator3DTrue = new TH3D(*aCorrFctn.fNumerator3DTrue);
  fNumerator3DFake = new TH3D(*aCorrFctn.fNumerator3DFake);
  fDenominator3D   = new TH3D(*aCorrFctn.fDenominator3D);
  fNumerator3DTrueIdeal = new TH3D(*aCorrFctn.fNumerator3DTrueIdeal);
  fNumerator3DFakeIdeal = new TH3D(*aCorrFctn.fNumerator3DFakeIdeal);
  fDenominator3DIdeal   = new TH3D(*aCorrFctn.fDenominator3DIdeal);
  fQinvHisto       = new TH3D(*aCorrFctn.fQinvHisto);
  fPairCut         = aCorrFctn.fPairCut->Clone();
}
//____________________________
AliFemtoModelBPLCMSCorrFctnKK::~AliFemtoModelBPLCMSCorrFctnKK()
{
  // destructor
  if (fNumeratorTrue) delete fNumeratorTrue;
  if (fNumeratorFake) delete fNumeratorFake;
  if (fDenominator) delete fDenominator;
  delete fNumerator3DTrue;
  delete fNumerator3DFake;
  delete fDenominator3D;
  delete fNumerator3DTrueIdeal;
  delete fNumerator3DFakeIdeal;
  delete fDenominator3DIdeal;
  delete fQinvHisto;
  if (fPairCut) delete fPairCut;
}
//_________________________
AliFemtoModelBPLCMSCorrFctnKK& AliFemtoModelBPLCMSCorrFctnKK::operator=(const AliFemtoModelBPLCMSCorrFctnKK& aCorrFctn)
{
  // Assignment operator
  if (this == &aCorrFctn)
    return *this;
  if (fNumerator3DTrue) delete fNumerator3DTrue;
  fNumerator3DTrue = new TH3D(*aCorrFctn.fNumerator3DTrue);
  if (fNumerator3DFake) delete fNumerator3DFake;
  fNumerator3DFake = new TH3D(*aCorrFctn.fNumerator3DFake);
  if (fDenominator3D) delete fDenominator3D;
  fDenominator3D = new TH3D(*aCorrFctn.fDenominator3D);

 if (fNumerator3DTrueIdeal) delete fNumerator3DTrueIdeal;
  fNumerator3DTrueIdeal = new TH3D(*aCorrFctn.fNumerator3DTrueIdeal);
 if (fNumerator3DFakeIdeal) delete fNumerator3DFakeIdeal;
  fNumerator3DFakeIdeal = new TH3D(*aCorrFctn.fNumerator3DFakeIdeal);
if (fDenominator3DIdeal) delete fDenominator3DIdeal;
  fDenominator3DIdeal = new TH3D(*aCorrFctn.fDenominator3DIdeal);
  if (fQinvHisto) delete fQinvHisto;
  fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
  fPairCut = aCorrFctn.fPairCut->Clone();

  return *this;
}

//_________________________
void AliFemtoModelBPLCMSCorrFctnKK::Write(){
  // Write out data histograms
  AliFemtoModelCorrFctn::Write();
  fNumerator3DTrue->Write();
  fNumerator3DFake->Write();
  fDenominator3D->Write();
  fNumerator3DTrueIdeal->Write();
  fNumerator3DFakeIdeal->Write();
  fDenominator3DIdeal->Write();
  fQinvHisto->Write();
}
//________________________
TList* AliFemtoModelBPLCMSCorrFctnKK::GetOutputList()
{
  // Prepare the list of objects to be written to the output
  TList *tOutputList = AliFemtoModelCorrFctn::GetOutputList();

  tOutputList->Add(fNumerator3DTrue);
  tOutputList->Add(fNumerator3DFake);
  tOutputList->Add(fDenominator3D);
  tOutputList->Add(fNumerator3DTrueIdeal);
  tOutputList->Add(fNumerator3DFakeIdeal);
  tOutputList->Add(fDenominator3DIdeal);
  tOutputList->Add(fQinvHisto);

  return tOutputList;
}

//_________________________
void AliFemtoModelBPLCMSCorrFctnKK::Finish(){
  fQinvHisto->Divide(fDenominator);
}

//____________________________
AliFemtoString AliFemtoModelBPLCMSCorrFctnKK::Report(){
  // Prepare a report from the execution
  string stemp = "LCMS Frame Bertsch-Pratt 3D Model Correlation Function Report:\n";
  char ctemp[100];
  snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fNumeratorTrue->GetEntries());
  stemp += ctemp;
  snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
  stemp += ctemp;

  /*  if (fCorrection)
      {
      float radius = fCorrection->GetRadius();
      snprintf(ctemp , 100, "Coulomb correction used radius of\t%E\n",radius);
      }
      else
      {
      snprintf(ctemp , 100, "No Coulomb Correction applied to this CorrFctn\n");
      }
      stemp += ctemp;
  */

  //
  AliFemtoString returnThis = stemp;
  return returnThis;
}
//____________________________
void AliFemtoModelBPLCMSCorrFctnKK::AddRealPair( AliFemtoPair* pair)
{
  // Store a real pair in numerator
  if (fPairCut){
    if (fUseRPSelection) {
      AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
      if (!ktc) {
	cout << "RP aware cut requested, but not connected to the CF" << endl;
	if (!(fPairCut->Pass(pair))) return;
      }
      else {
	AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
	if (!arp) {
	  cout << "RP aware cut requested, but not connected to the CF" << endl;
	  if (!(fPairCut->Pass(pair))) return;
	}
	else if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
      }
    }
    else
      if (!(fPairCut->Pass(pair))) return;
  }
//   if (fPairCut){
//     if (!(fPairCut->Pass(pair))) return;
//   }

  Double_t weight = fManager->GetWeight(pair);

  double qOut = (pair->QOutCMS());
  double qSide = (pair->QSideCMS());
  double qLong = (pair->QLongCMS());

  double qOutTrue = GetQoutTrue(pair);
  double qSideTrue = GetQsideTrue(pair);
  double qLongTrue =  GetQlongTrue(pair);


  fNumerator3DTrue->Fill(qOut, qSide, qLong, weight);
  fNumeratorTrue->Fill(pair->QInv(), weight);

  fNumerator3DTrueIdeal->Fill(qOutTrue, qSideTrue, qLongTrue, weight);

}
//____________________________
void AliFemtoModelBPLCMSCorrFctnKK::AddMixedPair( AliFemtoPair* pair){
  // store mixed pair in denominator
  if (fPairCut){
    if (fUseRPSelection) {
      AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
      if (!ktc) {
	cout << "RP aware cut requested, but not connected to the CF" << endl;
	if (!(fPairCut->Pass(pair))) return;
      }
      else {
	AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
	if (!arp) {
	  cout << "RP aware cut requested, but not connected to the CF" << endl;
	  if (!(fPairCut->Pass(pair))) return;
	}
	else if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
      }
    }
    else
      if (!(fPairCut->Pass(pair))) return;
  }
//   if (fPairCut){
//     if (!(fPairCut->Pass(pair))) return;
//   }

  Double_t weight = fManager->GetWeight(pair);

  double qOut = (pair->QOutCMS());
  double qSide = (pair->QSideCMS());
  double qLong = (pair->QLongCMS());

  double qOutTrue = GetQoutTrue(pair);
  double qSideTrue = GetQsideTrue(pair);
  double qLongTrue =  GetQlongTrue(pair);

  fNumerator3DFake->Fill(qOut, qSide, qLong, weight);
  fDenominator3D->Fill(qOut, qSide, qLong, 1.0);
  fNumeratorFake->Fill(pair->QInv(), weight);
  fDenominator->Fill(pair->QInv(), 1.0);

  fNumerator3DFakeIdeal->Fill(qOutTrue, qSideTrue, qLongTrue, weight);
  fDenominator3DIdeal->Fill(qOutTrue, qSideTrue, qLongTrue, 1.0);

}
//_______________________
AliFemtoModelCorrFctn* AliFemtoModelBPLCMSCorrFctnKK::Clone()
{
  // Clone the correlation function
  AliFemtoModelBPLCMSCorrFctnKK *tCopy = new AliFemtoModelBPLCMSCorrFctnKK(*this);

  return tCopy;
}

void AliFemtoModelBPLCMSCorrFctnKK::SetSpecificPairCut(AliFemtoPairCut* aCut)
{
  fPairCut = aCut;
}

void AliFemtoModelBPLCMSCorrFctnKK::SetUseRPSelection(unsigned short aRPSel)
{
  fUseRPSelection = aRPSel;
}


//_______________________
Double_t AliFemtoModelBPLCMSCorrFctnKK::GetQinvTrue(AliFemtoPair* aPair)
{
  AliFemtoTrack *inf1 = (AliFemtoTrack *) aPair->Track1()->Track();
  AliFemtoTrack *inf2 = (AliFemtoTrack *) aPair->Track2()->Track();

  AliFemtoLorentzVector fm1;
  AliFemtoThreeVector* temp = inf1->GetTrueMomentum();
  fm1.SetVect(*temp);
  double ener = TMath::Sqrt(temp->Mag2()+(aPair->Track1()->Track()->GetMass())*(aPair->Track1()->Track()->GetMass()));
  fm1.SetE(ener);

  AliFemtoLorentzVector fm2;
  AliFemtoThreeVector* temp2 = inf2->GetTrueMomentum();
  fm2.SetVect(*temp2);
  ener = TMath::Sqrt(temp2->Mag2()+(aPair->Track2()->Track()->GetMass())*(aPair->Track2()->Track()->GetMass()));
  fm2.SetE(ener);

  AliFemtoLorentzVector tQinvTrueVec = (fm1-fm2);
  Double_t tQinvTrue = -1.* tQinvTrueVec.m();

  return tQinvTrue;
}

//_______________________
Double_t AliFemtoModelBPLCMSCorrFctnKK::GetQoutTrue(AliFemtoPair* aPair)
{

 // relative momentum out component in lab frame
  AliFemtoTrack *inf1 = (AliFemtoTrack *) aPair->Track1()->Track();
  AliFemtoTrack *inf2 = (AliFemtoTrack *) aPair->Track2()->Track();

  // AliFemtoLorentzVector fm1;
  AliFemtoThreeVector* tmp1 = inf1->GetTrueMomentum();
  /*
 fm1.SetVect(*tmp1);
  double ener = TMath::Sqrt(tmp1->Mag2()+(aPair->Track1()->Track()->GetMass())*(aPair->Track1()->Track()->GetMass()));
  fm1.SetE(ener);
  */
  //  AliFemtoLorentzVector fm2;
  AliFemtoThreeVector* tmp2 = inf2->GetTrueMomentum();
  /*
  fm2.SetVect(*tmp2);
  ener = TMath::Sqrt(tmp2->Mag2()+(aPair->Track2()->Track()->GetMass())*(aPair->Track2()->Track()->GetMass()));
  fm2.SetE(ener);
  */
    double dx = tmp1->x() - tmp2->x();
    double xt = tmp1->x() + tmp2->x();

    double dy = tmp1->y() - tmp2->y();
    double yt = tmp1->y() + tmp2->y();

    double k1 = (::sqrt(xt*xt+yt*yt));
    double k2 = (dx*xt+dy*yt);
    double tmp;

    if(k1!=0) tmp= k2/k1;
    else tmp=0;

    return (tmp);

}


//_________________
Double_t AliFemtoModelBPLCMSCorrFctnKK::GetQsideTrue(AliFemtoPair* aPair)
{

  AliFemtoTrack *inf1 = (AliFemtoTrack *) aPair->Track1()->Track();
  AliFemtoTrack *inf2 = (AliFemtoTrack *) aPair->Track2()->Track();
  AliFemtoThreeVector* tmp1 = inf1->GetTrueMomentum();
  AliFemtoThreeVector* tmp2 = inf2->GetTrueMomentum();

  // relative momentum side component in lab frame

    double x1 = tmp1->x();  double y1 = tmp1->y();
    double x2 = tmp2->x();  double y2 = tmp2->y();

    double xt = x1+x2;  double yt = y1+y2;
    double k1 = ::sqrt(xt*xt+yt*yt);

    double tmp;
    if(k1!=0) tmp= 2.0*(x2*y1-x1*y2)/k1;
    else tmp=0;

    return (tmp);
}

//_________________________
double AliFemtoModelBPLCMSCorrFctnKK::GetQlongTrue(AliFemtoPair* aPair)
{
  // relative momentum component in lab frame

  AliFemtoTrack *inf1 = (AliFemtoTrack *) aPair->Track1()->Track();
  AliFemtoTrack *inf2 = (AliFemtoTrack *) aPair->Track2()->Track();
  AliFemtoThreeVector* temp1 = inf1->GetTrueMomentum();
  AliFemtoThreeVector* temp2 = inf2->GetTrueMomentum();

  AliFemtoLorentzVector tmp1;
  tmp1.SetVect(*temp1);
 double ener = TMath::Sqrt(temp1->Mag2()+(aPair->Track1()->Track()->GetMass())*(aPair->Track1()->Track()->GetMass()));
  tmp1.SetE(ener);

AliFemtoLorentzVector tmp2;
  tmp2.SetVect(*temp2);
 double ener2 = TMath::Sqrt(temp2->Mag2()+(aPair->Track2()->Track()->GetMass())*(aPair->Track2()->Track()->GetMass()));
  tmp2.SetE(ener2);

    double dz = tmp1.z() - tmp2.z();
    double zz = tmp1.z() + tmp2.z();

    double dt = tmp1.t() - tmp2.t();
    double tt = tmp1.t() + tmp2.t();

    double beta = zz/tt;
    double gamma = 1.0/TMath::Sqrt((1.-beta)*(1.+beta));

    double temp = gamma*(dz - beta*dt);
    return (temp);
}
Commit	Line	Data
c47b1f4b	1	////////////////////////////////////////////////////////////////////////////////
	2	/// ///
	3	/// AliFemtoModelBPLCMSCorrFctnKK - the class for correlation function which ///
	4	/// uses the model framework and weight generation and calculated the 3D ///
	5	/// correlation function in the Bertsh-Pratt LCMS system ///
	6	/// Authors: Adam Kisiel, kisiel@mps.ohio-state.edu ///
	7	/// ///
	8	////////////////////////////////////////////////////////////////////////////////
	9	#include "AliFemtoModelBPLCMSCorrFctnKK.h"
	10	#include "AliFemtoPair.h"
	11	#include "AliFemtoModelManager.h"
	12	#include "AliFemtoKTPairCut.h"
	13	#include "AliFemtoAnalysisReactionPlane.h"
	14	#include <cstdio>
	15
348cfdf7	16	#ifdef __ROOT__
c47b1f4b	17	ClassImp(AliFemtoModelBPLCMSCorrFctnKK)
	18	#endif
	19
	20	//____________________________
	21	AliFemtoModelBPLCMSCorrFctnKK::AliFemtoModelBPLCMSCorrFctnKK(char* title, const int& nbins, const float& QLo, const float& QHi)
	22	:
	23	AliFemtoModelCorrFctn(title, nbins, QLo, QHi),
	24	fNumerator3DTrue(0),
	25	fNumerator3DFake(0),
	26	fDenominator3D(0),
	27	fQinvHisto(0),
	28	fPairCut(0),
	29	fUseRPSelection(0),
	30	fNumerator3DTrueIdeal(0),
	31	fNumerator3DFakeIdeal(0),
	32	fDenominator3DIdeal(0)
	33	{
	34	// set up true numerator
	35	char tTitNumT[101] = "Num3DTrue";
348cfdf7	36	strncat(tTitNumT,title, 101);
c47b1f4b	37	fNumerator3DTrue = new TH3D(tTitNumT,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	38	// set up fake numerator
	39	char tTitNumF[101] = "Num3DFake";
	40	strncat(tTitNumF,title, 100);
	41	fNumerator3DFake = new TH3D(tTitNumF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	42	// set up denominator
	43	char tTitDen[101] = "Den3D";
	44	strncat(tTitDen,title, 100);
	45	fDenominator3D = new TH3D(tTitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	46	// set up ave qInv
	47	char tTitQinv[101] = "Qinv";
	48	strncat(tTitQinv,title, 100);
	49	fQinvHisto = new TH3D(tTitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	50
	51	// set up true Ideal numerator
	52	char tTitNumTI[101] = "Num3DTrueIdeal";
	53	strncat(tTitNumTI,title, 100);
	54	fNumerator3DTrueIdeal = new TH3D(tTitNumTI,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	55
	56	// set up fake Ideal numerator
	57	char tTitNumFI[101] = "Num3DFakeIdeal";
	58	strncat(tTitNumFI,title, 100);
	59	fNumerator3DFakeIdeal = new TH3D(tTitNumFI,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	60
	61	// set up denominator
	62	char tTitDenI[101] = "Den3DIdeal";
	63	strncat(tTitDenI,title, 100);
	64	fDenominator3DIdeal = new TH3D(tTitDenI,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	65
	66
	67	// to enable error bar calculation...
	68	fNumerator3DTrue->Sumw2();
	69	fNumerator3DFake->Sumw2();
	70	fDenominator3D->Sumw2();
	71	fNumerator3DTrueIdeal->Sumw2();
	72	fNumerator3DFakeIdeal->Sumw2();
	73	fDenominator3DIdeal->Sumw2();
	74	}
	75
	76	AliFemtoModelBPLCMSCorrFctnKK::AliFemtoModelBPLCMSCorrFctnKK(const AliFemtoModelBPLCMSCorrFctnKK& aCorrFctn) :
	77	AliFemtoModelCorrFctn(aCorrFctn),
	78	fNumerator3DTrue(0),
	79	fNumerator3DFake(0),
	80	fDenominator3D(0),
	81	fQinvHisto(0),
	82	fPairCut(0),
	83	fUseRPSelection(0),
	84	fNumerator3DTrueIdeal(0),
	85	fNumerator3DFakeIdeal(0),
	86	fDenominator3DIdeal(0)
	87	{
	88	// Copy constructor
	89	fNumerator3DTrue = new TH3D(*aCorrFctn.fNumerator3DTrue);
	90	fNumerator3DFake = new TH3D(*aCorrFctn.fNumerator3DFake);
	91	fDenominator3D = new TH3D(*aCorrFctn.fDenominator3D);
	92	fNumerator3DTrueIdeal = new TH3D(*aCorrFctn.fNumerator3DTrueIdeal);
	93	fNumerator3DFakeIdeal = new TH3D(*aCorrFctn.fNumerator3DFakeIdeal);
	94	fDenominator3DIdeal = new TH3D(*aCorrFctn.fDenominator3DIdeal);
	95	fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
	96	fPairCut = aCorrFctn.fPairCut->Clone();
	97	}
	98	//____________________________
	99	AliFemtoModelBPLCMSCorrFctnKK::~AliFemtoModelBPLCMSCorrFctnKK()
	100	{
101	// destructor
102	if (fNumeratorTrue) delete fNumeratorTrue;
103	if (fNumeratorFake) delete fNumeratorFake;
104	if (fDenominator) delete fDenominator;
105	delete fNumerator3DTrue;
106	delete fNumerator3DFake;
107	delete fDenominator3D;
108	delete fNumerator3DTrueIdeal;
109	delete fNumerator3DFakeIdeal;
110	delete fDenominator3DIdeal;
111	delete fQinvHisto;
112	if (fPairCut) delete fPairCut;
113	}
114	//_________________________
115	AliFemtoModelBPLCMSCorrFctnKK& AliFemtoModelBPLCMSCorrFctnKK::operator=(const AliFemtoModelBPLCMSCorrFctnKK& aCorrFctn)
116	{
117	// Assignment operator
118	if (this == &aCorrFctn)
119	return *this;
120	if (fNumerator3DTrue) delete fNumerator3DTrue;
121	fNumerator3DTrue = new TH3D(*aCorrFctn.fNumerator3DTrue);
122	if (fNumerator3DFake) delete fNumerator3DFake;
123	fNumerator3DFake = new TH3D(*aCorrFctn.fNumerator3DFake);
124	if (fDenominator3D) delete fDenominator3D;
125	fDenominator3D = new TH3D(*aCorrFctn.fDenominator3D);
126
127	if (fNumerator3DTrueIdeal) delete fNumerator3DTrueIdeal;
128	fNumerator3DTrueIdeal = new TH3D(*aCorrFctn.fNumerator3DTrueIdeal);
129	if (fNumerator3DFakeIdeal) delete fNumerator3DFakeIdeal;
130	fNumerator3DFakeIdeal = new TH3D(*aCorrFctn.fNumerator3DFakeIdeal);
131	if (fDenominator3DIdeal) delete fDenominator3DIdeal;
132	fDenominator3DIdeal = new TH3D(*aCorrFctn.fDenominator3DIdeal);
133	if (fQinvHisto) delete fQinvHisto;
134	fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
135	fPairCut = aCorrFctn.fPairCut->Clone();
136
137	return *this;
138	}
139
140	//_________________________
141	void AliFemtoModelBPLCMSCorrFctnKK::Write(){
142	// Write out data histograms
143	AliFemtoModelCorrFctn::Write();
144	fNumerator3DTrue->Write();
145	fNumerator3DFake->Write();
146	fDenominator3D->Write();
147	fNumerator3DTrueIdeal->Write();
148	fNumerator3DFakeIdeal->Write();
149	fDenominator3DIdeal->Write();
150	fQinvHisto->Write();
151	}
152	//________________________
153	TList* AliFemtoModelBPLCMSCorrFctnKK::GetOutputList()
154	{
155	// Prepare the list of objects to be written to the output
156	TList *tOutputList = AliFemtoModelCorrFctn::GetOutputList();
157
348cfdf7	158	tOutputList->Add(fNumerator3DTrue);
	159	tOutputList->Add(fNumerator3DFake);
	160	tOutputList->Add(fDenominator3D);
	161	tOutputList->Add(fNumerator3DTrueIdeal);
	162	tOutputList->Add(fNumerator3DFakeIdeal);
	163	tOutputList->Add(fDenominator3DIdeal);
	164	tOutputList->Add(fQinvHisto);
c47b1f4b	165
	166	return tOutputList;
	167	}
	168
	169	//_________________________
	170	void AliFemtoModelBPLCMSCorrFctnKK::Finish(){
	171	fQinvHisto->Divide(fDenominator);
	172	}
	173
	174	//____________________________
	175	AliFemtoString AliFemtoModelBPLCMSCorrFctnKK::Report(){
	176	// Prepare a report from the execution
	177	string stemp = "LCMS Frame Bertsch-Pratt 3D Model Correlation Function Report:\n";
	178	char ctemp[100];
	179	snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fNumeratorTrue->GetEntries());
	180	stemp += ctemp;
	181	snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
	182	stemp += ctemp;
	183
	184	/* if (fCorrection)
	185	{
	186	float radius = fCorrection->GetRadius();
	187	snprintf(ctemp , 100, "Coulomb correction used radius of\t%E\n",radius);
	188	}
	189	else
	190	{
	191	snprintf(ctemp , 100, "No Coulomb Correction applied to this CorrFctn\n");
	192	}
	193	stemp += ctemp;
	194	*/
	195
348cfdf7	196	//
c47b1f4b	197	AliFemtoString returnThis = stemp;
	198	return returnThis;
	199	}
	200	//____________________________
	201	void AliFemtoModelBPLCMSCorrFctnKK::AddRealPair( AliFemtoPair* pair)
	202	{
	203	// Store a real pair in numerator
	204	if (fPairCut){
	205	if (fUseRPSelection) {
	206	AliFemtoKTPairCut ktc = dynamic_cast<AliFemtoKTPairCut >(fPairCut);
348cfdf7	207	if (!ktc) {
c47b1f4b	208	cout << "RP aware cut requested, but not connected to the CF" << endl;
	209	if (!(fPairCut->Pass(pair))) return;
	210	}
	211	else {
	212	AliFemtoAnalysisReactionPlane arp = dynamic_cast<AliFemtoAnalysisReactionPlane > (HbtAnalysis());
	213	if (!arp) {
	214	cout << "RP aware cut requested, but not connected to the CF" << endl;
	215	if (!(fPairCut->Pass(pair))) return;
	216	}
	217	else if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
	218	}
	219	}
	220	else
	221	if (!(fPairCut->Pass(pair))) return;
	222	}
	223	// if (fPairCut){
	224	// if (!(fPairCut->Pass(pair))) return;
	225	// }
348cfdf7	226
c47b1f4b	227	Double_t weight = fManager->GetWeight(pair);
	228
	229	double qOut = (pair->QOutCMS());
	230	double qSide = (pair->QSideCMS());
	231	double qLong = (pair->QLongCMS());
348cfdf7	232
c47b1f4b	233	double qOutTrue = GetQoutTrue(pair);
	234	double qSideTrue = GetQsideTrue(pair);
	235	double qLongTrue = GetQlongTrue(pair);
	236
	237
	238	fNumerator3DTrue->Fill(qOut, qSide, qLong, weight);
	239	fNumeratorTrue->Fill(pair->QInv(), weight);
	240
	241	fNumerator3DTrueIdeal->Fill(qOutTrue, qSideTrue, qLongTrue, weight);
	242
	243	}
	244	//____________________________
	245	void AliFemtoModelBPLCMSCorrFctnKK::AddMixedPair( AliFemtoPair* pair){
	246	// store mixed pair in denominator
	247	if (fPairCut){
	248	if (fUseRPSelection) {
	249	AliFemtoKTPairCut ktc = dynamic_cast<AliFemtoKTPairCut >(fPairCut);
348cfdf7	250	if (!ktc) {
c47b1f4b	251	cout << "RP aware cut requested, but not connected to the CF" << endl;
	252	if (!(fPairCut->Pass(pair))) return;
	253	}
	254	else {
	255	AliFemtoAnalysisReactionPlane arp = dynamic_cast<AliFemtoAnalysisReactionPlane > (HbtAnalysis());
	256	if (!arp) {
	257	cout << "RP aware cut requested, but not connected to the CF" << endl;
	258	if (!(fPairCut->Pass(pair))) return;
	259	}
	260	else if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
	261	}
	262	}
	263	else
	264	if (!(fPairCut->Pass(pair))) return;
	265	}
	266	// if (fPairCut){
	267	// if (!(fPairCut->Pass(pair))) return;
	268	// }
	269
	270	Double_t weight = fManager->GetWeight(pair);
	271
	272	double qOut = (pair->QOutCMS());
	273	double qSide = (pair->QSideCMS());
	274	double qLong = (pair->QLongCMS());
	275
	276	double qOutTrue = GetQoutTrue(pair);
	277	double qSideTrue = GetQsideTrue(pair);
	278	double qLongTrue = GetQlongTrue(pair);
	279
	280	fNumerator3DFake->Fill(qOut, qSide, qLong, weight);
	281	fDenominator3D->Fill(qOut, qSide, qLong, 1.0);
	282	fNumeratorFake->Fill(pair->QInv(), weight);
	283	fDenominator->Fill(pair->QInv(), 1.0);
	284
	285	fNumerator3DFakeIdeal->Fill(qOutTrue, qSideTrue, qLongTrue, weight);
	286	fDenominator3DIdeal->Fill(qOutTrue, qSideTrue, qLongTrue, 1.0);
	287
	288	}
	289	//_______________________
	290	AliFemtoModelCorrFctn* AliFemtoModelBPLCMSCorrFctnKK::Clone()
	291	{
	292	// Clone the correlation function
	293	AliFemtoModelBPLCMSCorrFctnKK tCopy = new AliFemtoModelBPLCMSCorrFctnKK(this);
348cfdf7	294
c47b1f4b	295	return tCopy;
	296	}
	297
	298	void AliFemtoModelBPLCMSCorrFctnKK::SetSpecificPairCut(AliFemtoPairCut* aCut)
	299	{
	300	fPairCut = aCut;
	301	}
	302
	303	void AliFemtoModelBPLCMSCorrFctnKK::SetUseRPSelection(unsigned short aRPSel)
	304	{
	305	fUseRPSelection = aRPSel;
	306	}
	307
	308
	309	//_______________________
	310	Double_t AliFemtoModelBPLCMSCorrFctnKK::GetQinvTrue(AliFemtoPair* aPair)
	311	{
	312	AliFemtoTrack inf1 = (AliFemtoTrack ) aPair->Track1()->Track();
	313	AliFemtoTrack inf2 = (AliFemtoTrack ) aPair->Track2()->Track();
	314
	315	AliFemtoLorentzVector fm1;
	316	AliFemtoThreeVector* temp = inf1->GetTrueMomentum();
	317	fm1.SetVect(*temp);
	318	double ener = TMath::Sqrt(temp->Mag2()+(aPair->Track1()->Track()->GetMass())*(aPair->Track1()->Track()->GetMass()));
	319	fm1.SetE(ener);
	320
	321	AliFemtoLorentzVector fm2;
	322	AliFemtoThreeVector* temp2 = inf2->GetTrueMomentum();
	323	fm2.SetVect(*temp2);
	324	ener = TMath::Sqrt(temp2->Mag2()+(aPair->Track2()->Track()->GetMass())*(aPair->Track2()->Track()->GetMass()));
	325	fm2.SetE(ener);
	326
	327	AliFemtoLorentzVector tQinvTrueVec = (fm1-fm2);
	328	Double_t tQinvTrue = -1.* tQinvTrueVec.m();
	329
	330	return tQinvTrue;
	331	}
	332
	333	//_______________________
	334	Double_t AliFemtoModelBPLCMSCorrFctnKK::GetQoutTrue(AliFemtoPair* aPair)
	335	{
	336
	337	// relative momentum out component in lab frame
	338	AliFemtoTrack inf1 = (AliFemtoTrack ) aPair->Track1()->Track();
	339	AliFemtoTrack inf2 = (AliFemtoTrack ) aPair->Track2()->Track();
	340
	341	// AliFemtoLorentzVector fm1;
	342	AliFemtoThreeVector* tmp1 = inf1->GetTrueMomentum();
348cfdf7	343	/*
c47b1f4b	344	fm1.SetVect(*tmp1);
	345	double ener = TMath::Sqrt(tmp1->Mag2()+(aPair->Track1()->Track()->GetMass())*(aPair->Track1()->Track()->GetMass()));
	346	fm1.SetE(ener);
	347	*/
	348	// AliFemtoLorentzVector fm2;
	349	AliFemtoThreeVector* tmp2 = inf2->GetTrueMomentum();
	350	/*
	351	fm2.SetVect(*tmp2);
	352	ener = TMath::Sqrt(tmp2->Mag2()+(aPair->Track2()->Track()->GetMass())*(aPair->Track2()->Track()->GetMass()));
	353	fm2.SetE(ener);
	354	*/
	355	double dx = tmp1->x() - tmp2->x();
	356	double xt = tmp1->x() + tmp2->x();
348cfdf7	357
c47b1f4b	358	double dy = tmp1->y() - tmp2->y();
	359	double yt = tmp1->y() + tmp2->y();
	360
	361	double k1 = (::sqrt(xtxt+ytyt));
	362	double k2 = (dxxt+dyyt);
	363	double tmp;
	364
	365	if(k1!=0) tmp= k2/k1;
	366	else tmp=0;
	367
	368	return (tmp);
	369
	370	}
	371
	372
	373
	374	//_________________
	375	Double_t AliFemtoModelBPLCMSCorrFctnKK::GetQsideTrue(AliFemtoPair* aPair)
	376	{
	377
	378	AliFemtoTrack inf1 = (AliFemtoTrack ) aPair->Track1()->Track();
	379	AliFemtoTrack inf2 = (AliFemtoTrack ) aPair->Track2()->Track();
	380	AliFemtoThreeVector* tmp1 = inf1->GetTrueMomentum();
	381	AliFemtoThreeVector* tmp2 = inf2->GetTrueMomentum();
	382
	383	// relative momentum side component in lab frame
348cfdf7	384
c47b1f4b	385	double x1 = tmp1->x(); double y1 = tmp1->y();
	386	double x2 = tmp2->x(); double y2 = tmp2->y();
	387
	388	double xt = x1+x2; double yt = y1+y2;
	389	double k1 = ::sqrt(xtxt+ytyt);
348cfdf7	390
c47b1f4b	391	double tmp;
	392	if(k1!=0) tmp= 2.0(x2y1-x1*y2)/k1;
	393	else tmp=0;
	394
	395	return (tmp);
	396	}
	397
	398	//_________________________
	399	double AliFemtoModelBPLCMSCorrFctnKK::GetQlongTrue(AliFemtoPair* aPair)
	400	{
	401	// relative momentum component in lab frame
	402
	403	AliFemtoTrack inf1 = (AliFemtoTrack ) aPair->Track1()->Track();
	404	AliFemtoTrack inf2 = (AliFemtoTrack ) aPair->Track2()->Track();
	405	AliFemtoThreeVector* temp1 = inf1->GetTrueMomentum();
	406	AliFemtoThreeVector* temp2 = inf2->GetTrueMomentum();
	407
	408	AliFemtoLorentzVector tmp1;
	409	tmp1.SetVect(*temp1);
	410	double ener = TMath::Sqrt(temp1->Mag2()+(aPair->Track1()->Track()->GetMass())*(aPair->Track1()->Track()->GetMass()));
	411	tmp1.SetE(ener);
	412
	413	AliFemtoLorentzVector tmp2;
	414	tmp2.SetVect(*temp2);
	415	double ener2 = TMath::Sqrt(temp2->Mag2()+(aPair->Track2()->Track()->GetMass())*(aPair->Track2()->Track()->GetMass()));
	416	tmp2.SetE(ener2);
	417
	418	double dz = tmp1.z() - tmp2.z();
	419	double zz = tmp1.z() + tmp2.z();
	420
	421	double dt = tmp1.t() - tmp2.t();
	422	double tt = tmp1.t() + tmp2.t();
	423
	424	double beta = zz/tt;
	425	double gamma = 1.0/TMath::Sqrt((1.-beta)*(1.+beta));
	426
	427	double temp = gamma(dz - betadt);
	428	return (temp);
	429	}