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

////////////////////////////////////////////////////////////////////////////////
///                                                                          ///
/// AliFemtoModelBPLCMSCorrFctn - 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 "AliFemtoModelBPLCMSCorrFctn.h"
#include "AliFemtoPair.h"
#include "AliFemtoModelManager.h"
#include "AliFemtoKTPairCut.h"
#include "AliFemtoAnalysisReactionPlane.h"
#include <cstdio>

#ifdef __ROOT__
ClassImp(AliFemtoModelBPLCMSCorrFctn)
#endif

//____________________________
AliFemtoModelBPLCMSCorrFctn::AliFemtoModelBPLCMSCorrFctn(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)
{
  // set up true numerator
  char tTitNumT[101] = "Num3DTrue";
  strncat(tTitNumT,title, 100);
  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);

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

AliFemtoModelBPLCMSCorrFctn::AliFemtoModelBPLCMSCorrFctn(const AliFemtoModelBPLCMSCorrFctn& aCorrFctn) :
  AliFemtoModelCorrFctn(aCorrFctn),
  fNumerator3DTrue(0),
  fNumerator3DFake(0),
  fDenominator3D(0),
  fQinvHisto(0),
  fPairCut(0),
  fUseRPSelection(0)
{
  // Copy constructor
  fNumerator3DTrue = new TH3D(*aCorrFctn.fNumerator3DTrue);
  fNumerator3DFake = new TH3D(*aCorrFctn.fNumerator3DFake);
  fDenominator3D   = new TH3D(*aCorrFctn.fDenominator3D);
  fQinvHisto       = new TH3D(*aCorrFctn.fQinvHisto);
  fPairCut         = aCorrFctn.fPairCut->Clone();
}
//____________________________
AliFemtoModelBPLCMSCorrFctn::~AliFemtoModelBPLCMSCorrFctn()
{
  // destructor
  if (fNumeratorTrue) delete fNumeratorTrue;
  if (fNumeratorFake) delete fNumeratorFake;
  if (fDenominator) delete fDenominator;
  delete fNumerator3DTrue;
  delete fNumerator3DFake;
  delete fDenominator3D;
  delete fQinvHisto;
  if (fPairCut) delete fPairCut;
}
//_________________________
AliFemtoModelBPLCMSCorrFctn& AliFemtoModelBPLCMSCorrFctn::operator=(const AliFemtoModelBPLCMSCorrFctn& 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 (fQinvHisto) delete fQinvHisto;
  fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
  fPairCut = aCorrFctn.fPairCut->Clone();

  return *this;
}

//_________________________
void AliFemtoModelBPLCMSCorrFctn::Write(){
  // Write out data histograms
  AliFemtoModelCorrFctn::Write();
  fNumerator3DTrue->Write();
  fNumerator3DFake->Write();
  fDenominator3D->Write();
  fQinvHisto->Write();
}
//________________________
TList* AliFemtoModelBPLCMSCorrFctn::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(fQinvHisto);

  return tOutputList;
}

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

//____________________________
AliFemtoString AliFemtoModelBPLCMSCorrFctn::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 AliFemtoModelBPLCMSCorrFctn::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());

  fNumerator3DTrue->Fill(qOut, qSide, qLong, weight);
  fNumeratorTrue->Fill(pair->QInv(), weight);
}
//____________________________
void AliFemtoModelBPLCMSCorrFctn::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());

  fNumerator3DFake->Fill(qOut, qSide, qLong, weight);
  fDenominator3D->Fill(qOut, qSide, qLong, 1.0);
  fNumeratorFake->Fill(pair->QInv(), weight);
  fDenominator->Fill(pair->QInv(), 1.0);
  fQinvHisto->Fill(qOut, qSide, qLong, pair->QInv() );
}
//_______________________
AliFemtoModelCorrFctn* AliFemtoModelBPLCMSCorrFctn::Clone()
{
  // Clone the correlation function
  AliFemtoModelBPLCMSCorrFctn *tCopy = new AliFemtoModelBPLCMSCorrFctn(*this);

  return tCopy;
}

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

void AliFemtoModelBPLCMSCorrFctn::SetUseRPSelection(unsigned short aRPSel)
{
  fUseRPSelection = aRPSel;
}
Commit	Line	Data
76ce4b5b	1	////////////////////////////////////////////////////////////////////////////////
	2	/// ///
	3	/// AliFemtoModelBPLCMSCorrFctn - 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 "AliFemtoModelBPLCMSCorrFctn.h"
	10	#include "AliFemtoPair.h"
	11	#include "AliFemtoModelManager.h"
	12	#include "AliFemtoKTPairCut.h"
	13	#include "AliFemtoAnalysisReactionPlane.h"
	14	#include <cstdio>
	15
6f24ad05	16	#ifdef __ROOT__
76ce4b5b	17	ClassImp(AliFemtoModelBPLCMSCorrFctn)
	18	#endif
	19
	20	//____________________________
	21	AliFemtoModelBPLCMSCorrFctn::AliFemtoModelBPLCMSCorrFctn(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	{
	31	// set up true numerator
	32	char tTitNumT[101] = "Num3DTrue";
	33	strncat(tTitNumT,title, 100);
	34	fNumerator3DTrue = new TH3D(tTitNumT,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	35	// set up fake numerator
	36	char tTitNumF[101] = "Num3DFake";
	37	strncat(tTitNumF,title, 100);
	38	fNumerator3DFake = new TH3D(tTitNumF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	39	// set up denominator
	40	char tTitDen[101] = "Den3D";
	41	strncat(tTitDen,title, 100);
	42	fDenominator3D = new TH3D(tTitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	43	// set up ave qInv
	44	char tTitQinv[101] = "Qinv";
	45	strncat(tTitQinv,title, 100);
	46	fQinvHisto = new TH3D(tTitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
	47
	48	// to enable error bar calculation...
	49	fNumerator3DTrue->Sumw2();
	50	fNumerator3DFake->Sumw2();
	51	fDenominator3D->Sumw2();
	52	}
	53
	54	AliFemtoModelBPLCMSCorrFctn::AliFemtoModelBPLCMSCorrFctn(const AliFemtoModelBPLCMSCorrFctn& aCorrFctn) :
	55	AliFemtoModelCorrFctn(aCorrFctn),
	56	fNumerator3DTrue(0),
	57	fNumerator3DFake(0),
	58	fDenominator3D(0),
	59	fQinvHisto(0),
	60	fPairCut(0),
	61	fUseRPSelection(0)
	62	{
	63	// Copy constructor
	64	fNumerator3DTrue = new TH3D(*aCorrFctn.fNumerator3DTrue);
	65	fNumerator3DFake = new TH3D(*aCorrFctn.fNumerator3DFake);
	66	fDenominator3D = new TH3D(*aCorrFctn.fDenominator3D);
	67	fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
	68	fPairCut = aCorrFctn.fPairCut->Clone();
	69	}
	70	//____________________________
	71	AliFemtoModelBPLCMSCorrFctn::~AliFemtoModelBPLCMSCorrFctn()
	72	{
	73	// destructor
	74	if (fNumeratorTrue) delete fNumeratorTrue;
	75	if (fNumeratorFake) delete fNumeratorFake;
	76	if (fDenominator) delete fDenominator;
	77	delete fNumerator3DTrue;
	78	delete fNumerator3DFake;
	79	delete fDenominator3D;
	80	delete fQinvHisto;
81	if (fPairCut) delete fPairCut;
82	}
83	//_________________________
84	AliFemtoModelBPLCMSCorrFctn& AliFemtoModelBPLCMSCorrFctn::operator=(const AliFemtoModelBPLCMSCorrFctn& aCorrFctn)
85	{
86	// Assignment operator
87	if (this == &aCorrFctn)
88	return *this;
89	if (fNumerator3DTrue) delete fNumerator3DTrue;
90	fNumerator3DTrue = new TH3D(*aCorrFctn.fNumerator3DTrue);
91	if (fNumerator3DFake) delete fNumerator3DFake;
92	fNumerator3DFake = new TH3D(*aCorrFctn.fNumerator3DFake);
93	if (fDenominator3D) delete fDenominator3D;
94	fDenominator3D = new TH3D(*aCorrFctn.fDenominator3D);
95	if (fQinvHisto) delete fQinvHisto;
96	fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
97	fPairCut = aCorrFctn.fPairCut->Clone();
98
99	return *this;
100	}
101
102	//_________________________
103	void AliFemtoModelBPLCMSCorrFctn::Write(){
104	// Write out data histograms
105	AliFemtoModelCorrFctn::Write();
106	fNumerator3DTrue->Write();
107	fNumerator3DFake->Write();
108	fDenominator3D->Write();
109	fQinvHisto->Write();
110	}
111	//________________________
112	TList* AliFemtoModelBPLCMSCorrFctn::GetOutputList()
113	{
114	// Prepare the list of objects to be written to the output
115	TList *tOutputList = AliFemtoModelCorrFctn::GetOutputList();
116
6f24ad05	117	tOutputList->Add(fNumerator3DTrue);
	118	tOutputList->Add(fNumerator3DFake);
	119	tOutputList->Add(fDenominator3D);
	120	tOutputList->Add(fQinvHisto);
76ce4b5b	121
	122	return tOutputList;
	123	}
	124
	125	//_________________________
	126	void AliFemtoModelBPLCMSCorrFctn::Finish(){
109b4bb5	127	//fQinvHisto->Divide(fDenominator);
76ce4b5b	128	}
	129
	130	//____________________________
	131	AliFemtoString AliFemtoModelBPLCMSCorrFctn::Report(){
	132	// Prepare a report from the execution
	133	string stemp = "LCMS Frame Bertsch-Pratt 3D Model Correlation Function Report:\n";
	134	char ctemp[100];
	135	snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fNumeratorTrue->GetEntries());
	136	stemp += ctemp;
	137	snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
	138	stemp += ctemp;
	139	/* if (fCorrection)
	140	{
	141	float radius = fCorrection->GetRadius();
	142	snprintf(ctemp , 100, "Coulomb correction used radius of\t%E\n",radius);
	143	}
	144	else
	145	{
	146	snprintf(ctemp , 100, "No Coulomb Correction applied to this CorrFctn\n");
	147	}
	148	stemp += ctemp;
	149	*/
	150
6f24ad05	151	//
76ce4b5b	152	AliFemtoString returnThis = stemp;
	153	return returnThis;
	154	}
	155	//____________________________
	156	void AliFemtoModelBPLCMSCorrFctn::AddRealPair( AliFemtoPair* pair)
	157	{
	158	// Store a real pair in numerator
	159	if (fPairCut){
	160	if (fUseRPSelection) {
	161	AliFemtoKTPairCut ktc = dynamic_cast<AliFemtoKTPairCut >(fPairCut);
6f24ad05	162	if (!ktc) {
76ce4b5b	163	cout << "RP aware cut requested, but not connected to the CF" << endl;
	164	if (!(fPairCut->Pass(pair))) return;
	165	}
	166	else {
	167	AliFemtoAnalysisReactionPlane arp = dynamic_cast<AliFemtoAnalysisReactionPlane > (HbtAnalysis());
	168	if (!arp) {
	169	cout << "RP aware cut requested, but not connected to the CF" << endl;
	170	if (!(fPairCut->Pass(pair))) return;
	171	}
	172	else if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
	173	}
	174	}
	175	else
	176	if (!(fPairCut->Pass(pair))) return;
	177	}
	178	// if (fPairCut){
	179	// if (!(fPairCut->Pass(pair))) return;
	180	// }
6f24ad05	181
76ce4b5b	182	Double_t weight = fManager->GetWeight(pair);
	183
	184	double qOut = (pair->QOutCMS());
	185	double qSide = (pair->QSideCMS());
	186	double qLong = (pair->QLongCMS());
	187
	188	fNumerator3DTrue->Fill(qOut, qSide, qLong, weight);
	189	fNumeratorTrue->Fill(pair->QInv(), weight);
	190	}
	191	//____________________________
	192	void AliFemtoModelBPLCMSCorrFctn::AddMixedPair( AliFemtoPair* pair){
	193	// store mixed pair in denominator
	194	if (fPairCut){
	195	if (fUseRPSelection) {
	196	AliFemtoKTPairCut ktc = dynamic_cast<AliFemtoKTPairCut >(fPairCut);
6f24ad05	197	if (!ktc) {
109b4bb5	198	cout << "RP aware cut requested, but not connected to the CF" << endl;
109b4bb5	199	if (!(fPairCut->Pass(pair))) return;
76ce4b5b	200	}
76ce4b5b	201	else {
109b4bb5	202	AliFemtoAnalysisReactionPlane arp = dynamic_cast<AliFemtoAnalysisReactionPlane > (HbtAnalysis());
	203	if (!arp) {
	204	cout << "RP aware cut requested, but not connected to the CF" << endl;
	205	if (!(fPairCut->Pass(pair))) return;
	206	}
	207	else if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
76ce4b5b	208	}
	209	}
	210	else
	211	if (!(fPairCut->Pass(pair))) return;
	212	}
	213	// if (fPairCut){
	214	// if (!(fPairCut->Pass(pair))) return;
	215	// }
	216
	217	Double_t weight = fManager->GetWeight(pair);
	218
	219	double qOut = (pair->QOutCMS());
	220	double qSide = (pair->QSideCMS());
	221	double qLong = (pair->QLongCMS());
	222
	223	fNumerator3DFake->Fill(qOut, qSide, qLong, weight);
	224	fDenominator3D->Fill(qOut, qSide, qLong, 1.0);
	225	fNumeratorFake->Fill(pair->QInv(), weight);
	226	fDenominator->Fill(pair->QInv(), 1.0);
109b4bb5	227	fQinvHisto->Fill(qOut, qSide, qLong, pair->QInv() );
76ce4b5b	228	}
	229	//_______________________
	230	AliFemtoModelCorrFctn* AliFemtoModelBPLCMSCorrFctn::Clone()
	231	{
	232	// Clone the correlation function
	233	AliFemtoModelBPLCMSCorrFctn tCopy = new AliFemtoModelBPLCMSCorrFctn(this);
6f24ad05	234
76ce4b5b	235	return tCopy;
	236	}
	237
	238	void AliFemtoModelBPLCMSCorrFctn::SetSpecificPairCut(AliFemtoPairCut* aCut)
	239	{
	240	fPairCut = aCut;
	241	}
	242
	243	void AliFemtoModelBPLCMSCorrFctn::SetUseRPSelection(unsigned short aRPSel)
	244	{
	245	fUseRPSelection = aRPSel;
	246	}