[u/mrichter/AliRoot.git] / PWG2 / 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[100] = "Num3DTrue";
  strcat(tTitNumT,title);
  fNumerator3DTrue = new TH3D(tTitNumT,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  // set up fake numerator
  char tTitNumF[100] = "Num3DFake";
  strcat(tTitNumF,title);
  fNumerator3DFake = new TH3D(tTitNumF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  // set up denominator
  char tTitDen[100] = "Den3D";
  strcat(tTitDen,title);
  fDenominator3D = new TH3D(tTitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
  // set up ave qInv
  char tTitQinv[100] = "Qinv";
  strcat(tTitQinv,title);
  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];
  sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumeratorTrue->GetEntries());
  stemp += ctemp;
  sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
  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;
  */

  //  
  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;
	}
	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;
	}
	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);

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