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

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
// AliFemtoCorrFctnNonIdDR - correlation function for non-identical particles //
// uses k* as a function variable. Stores the correlation function separately //
// for positive and negative signs of k* projections into out, side and long  //
// directions, enabling the calculations of double ratios                     //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

#include "AliFemtoCorrFctnNonIdDR.h"
//#include "AliFemtoHisto.h"
#include <cstdio>

#ifdef __ROOT__ 
ClassImp(AliFemtoCorrFctnNonIdDR)
#endif

//____________________________
AliFemtoCorrFctnNonIdDR::AliFemtoCorrFctnNonIdDR(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
  fNumOutP(0), 
  fNumOutN(0),  
  fNumSideP(0), 
  fNumSideN(0), 
  fNumLongP(0), 
  fNumLongN(0), 
  fDenOutP(0),  
  fDenOutN(0),  
  fDenSideP(0), 
  fDenSideN(0), 
  fDenLongP(0), 
  fDenLongN(0)
{
  // Default constructor
  // set up numerators
  char bufname[200];
  snprintf(bufname, 200, "NumOutP%s", title);
  fNumOutP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  snprintf(bufname, 200, "NumOutN%s", title);
  fNumOutN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  snprintf(bufname, 200, "NumSideP%s", title);
  fNumSideP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  snprintf(bufname, 200, "NumSideN%s", title);
  fNumSideN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  snprintf(bufname, 200, "NumLongP%s", title);
  fNumLongP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  snprintf(bufname, 200, "NumLongN%s", title);
  fNumLongN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);

  // set up denominators
  snprintf(bufname, 200, "DenOutP%s", title);
  fDenOutP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  snprintf(bufname, 200, "DenOutN%s", title);
  fDenOutN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  snprintf(bufname, 200, "DenSideP%s", title);
  fDenSideP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  snprintf(bufname, 200, "DenSideN%s", title);
  fDenSideN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  snprintf(bufname, 200, "DenLongP%s", title);
  fDenLongP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
  snprintf(bufname, 200, "DenLongN%s", title);
  fDenLongN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);

  // to enable error bar calculation...
  fNumOutP->Sumw2(); 
  fNumOutN->Sumw2();  
  fNumSideP->Sumw2(); 
  fNumSideN->Sumw2(); 
  fNumLongP->Sumw2(); 
  fNumLongN->Sumw2(); 
  fDenOutP->Sumw2();  
  fDenOutN->Sumw2();  
  fDenSideP->Sumw2(); 
  fDenSideN->Sumw2(); 
  fDenLongP->Sumw2(); 
  fDenLongN->Sumw2();
}

//____________________________
AliFemtoCorrFctnNonIdDR::AliFemtoCorrFctnNonIdDR(const AliFemtoCorrFctnNonIdDR& aCorrFctn) :
  fNumOutP(0), 
  fNumOutN(0),  
  fNumSideP(0), 
  fNumSideN(0), 
  fNumLongP(0), 
  fNumLongN(0), 
  fDenOutP(0),  
  fDenOutN(0),  
  fDenSideP(0), 
  fDenSideN(0), 
  fDenLongP(0), 
  fDenLongN(0)
{
  // copy constructor
  if (aCorrFctn.fNumOutP)
    fNumOutP = new TH1D(*aCorrFctn.fNumOutP);
  if (aCorrFctn.fNumOutN)
    fNumOutN = new TH1D(*aCorrFctn.fNumOutN);
  if (aCorrFctn.fNumSideP)
    fNumSideP = new TH1D(*aCorrFctn.fNumSideP);
  if (aCorrFctn.fNumSideN)
    fNumSideN = new TH1D(*aCorrFctn.fNumSideN);
  if (aCorrFctn.fNumLongP)
    fNumLongP = new TH1D(*aCorrFctn.fNumLongP);
  if (aCorrFctn.fNumLongN)
    fNumLongN = new TH1D(*aCorrFctn.fNumLongN);

  if (aCorrFctn.fDenOutP)
    fDenOutP = new TH1D(*aCorrFctn.fDenOutP);
  if (aCorrFctn.fDenOutN)
    fDenOutN = new TH1D(*aCorrFctn.fDenOutN);
  if (aCorrFctn.fDenSideP)
    fDenSideP = new TH1D(*aCorrFctn.fDenSideP);
  if (aCorrFctn.fDenSideN)
    fDenSideN = new TH1D(*aCorrFctn.fDenSideN);
  if (aCorrFctn.fDenLongP)
    fDenLongP = new TH1D(*aCorrFctn.fDenLongP);
  if (aCorrFctn.fDenLongN)
    fDenLongN = new TH1D(*aCorrFctn.fDenLongN);
}
//____________________________
AliFemtoCorrFctnNonIdDR::~AliFemtoCorrFctnNonIdDR(){
  delete fNumOutP; 
  delete fNumOutN;  
  delete fNumSideP; 
  delete fNumSideN; 
  delete fNumLongP; 
  delete fNumLongN; 
  delete fDenOutP;  
  delete fDenOutN;  
  delete fDenSideP; 
  delete fDenSideN; 
  delete fDenLongP; 
  delete fDenLongN;
}
//_________________________
AliFemtoCorrFctnNonIdDR& AliFemtoCorrFctnNonIdDR::operator=(const AliFemtoCorrFctnNonIdDR& aCorrFctn)
{
  // assignment operator
  if (this == &aCorrFctn)
    return *this;

  if (aCorrFctn.fNumOutP)
    fNumOutP = new TH1D(*aCorrFctn.fNumOutP);
  if (aCorrFctn.fNumOutN)
    fNumOutN = new TH1D(*aCorrFctn.fNumOutN);
  if (aCorrFctn.fNumSideP)
    fNumSideP = new TH1D(*aCorrFctn.fNumSideP);
  if (aCorrFctn.fNumSideN)
    fNumSideN = new TH1D(*aCorrFctn.fNumSideN);
  if (aCorrFctn.fNumLongP)
    fNumLongP = new TH1D(*aCorrFctn.fNumLongP);
  if (aCorrFctn.fNumLongN)
    fNumLongN = new TH1D(*aCorrFctn.fNumLongN);

  if (aCorrFctn.fDenOutP)
    fDenOutP = new TH1D(*aCorrFctn.fDenOutP);
  if (aCorrFctn.fDenOutN)
    fDenOutN = new TH1D(*aCorrFctn.fDenOutN);
  if (aCorrFctn.fDenSideP)
    fDenSideP = new TH1D(*aCorrFctn.fDenSideP);
  if (aCorrFctn.fDenSideN)
    fDenSideN = new TH1D(*aCorrFctn.fDenSideN);
  if (aCorrFctn.fDenLongP)
    fDenLongP = new TH1D(*aCorrFctn.fDenLongP);
  if (aCorrFctn.fDenLongN)
    fDenLongN = new TH1D(*aCorrFctn.fDenLongN);

  return *this;
}

//_________________________
void AliFemtoCorrFctnNonIdDR::Finish(){
  // here is where we should normalize, fit, etc...
  // we should NOT Draw() the histos (as I had done it below),
  // since we want to insulate ourselves from root at this level
  // of the code.  Do it instead at root command line with browser.
  //  fNumerator->Draw();
  //fDenominator->Draw();
  //fRatio->Draw();
  //  fRatio->Divide(fNumerator,fDenominator,1.0,1.0);

}

//____________________________
AliFemtoString AliFemtoCorrFctnNonIdDR::Report(){
  // construct report
  string stemp = "Non-identical particles Correlation Function Report:\n";
  char ctemp[100];
  sprintf(ctemp,"Number of entries in numerators:\t%E\n",fNumOutP->GetEntries()+fNumOutN->GetEntries());
  stemp += ctemp;
  sprintf(ctemp,"Number of entries in denominators:\t%E\n",fDenOutP->GetEntries()+fDenOutN->GetEntries());
  stemp += ctemp;
  //  stemp += mCoulombWeight->Report();
  AliFemtoString returnThis = stemp;
  return returnThis;
}
//____________________________
void AliFemtoCorrFctnNonIdDR::AddRealPair(AliFemtoPair* pair){
  // add true pair
  double tKStar = pair->KStar();
  if (pair->KOut()>0.0)
    fNumOutP->Fill(tKStar);
  else
    fNumOutN->Fill(tKStar);

  if (pair->KSide()>0.0)
    fNumSideP->Fill(tKStar);
  else
    fNumSideN->Fill(tKStar);

  if (pair->KLong()>0.0)
    fNumLongP->Fill(tKStar);
  else
    fNumLongN->Fill(tKStar);

}
//____________________________
void AliFemtoCorrFctnNonIdDR::AddMixedPair(AliFemtoPair* pair){
  // add mixed (background) pair
  double tKStar = pair->KStar();
  if (pair->KOut()>0.0)
    fDenOutP->Fill(tKStar);
  else
    fDenOutN->Fill(tKStar);

  if (pair->KSide()>0.0)
    fDenSideP->Fill(tKStar);
  else
    fDenSideN->Fill(tKStar);

  if (pair->KLong()>0.0)
    fDenLongP->Fill(tKStar);
  else
    fDenLongN->Fill(tKStar);
}
//____________________________
void AliFemtoCorrFctnNonIdDR::Write(){
  fNumOutP->Write(); 
  fNumOutN->Write();  
  fNumSideP->Write(); 
  fNumSideN->Write(); 
  fNumLongP->Write(); 
  fNumLongN->Write(); 
  fDenOutP->Write();  
  fDenOutN->Write();  
  fDenSideP->Write(); 
  fDenSideN->Write(); 
  fDenLongP->Write(); 
  fDenLongN->Write();
}

TList* AliFemtoCorrFctnNonIdDR::GetOutputList()
{
  // Prepare the list of objects to be written to the output
  TList *tOutputList = new TList();

  tOutputList->Add(fNumOutP); 
  tOutputList->Add(fNumOutN);  
  tOutputList->Add(fNumSideP); 
  tOutputList->Add(fNumSideN); 
  tOutputList->Add(fNumLongP); 
  tOutputList->Add(fNumLongN); 
  tOutputList->Add(fDenOutP);  
  tOutputList->Add(fDenOutN);  
  tOutputList->Add(fDenSideP); 
  tOutputList->Add(fDenSideN); 
  tOutputList->Add(fDenLongP); 
  tOutputList->Add(fDenLongN);

  return tOutputList;
}
Commit	Line	Data
0b3bd1ac	1	////////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoCorrFctnNonIdDR - correlation function for non-identical particles //
	4	// uses k* as a function variable. Stores the correlation function separately //
	5	// for positive and negative signs of k* projections into out, side and long //
	6	// directions, enabling the calculations of double ratios //
	7	// //
	8	////////////////////////////////////////////////////////////////////////////////
	9
	10	#include "AliFemtoCorrFctnNonIdDR.h"
	11	//#include "AliFemtoHisto.h"
	12	#include <cstdio>
	13
	14	#ifdef __ROOT__
	15	ClassImp(AliFemtoCorrFctnNonIdDR)
	16	#endif
	17
	18	//____________________________
	19	AliFemtoCorrFctnNonIdDR::AliFemtoCorrFctnNonIdDR(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
	20	fNumOutP(0),
	21	fNumOutN(0),
	22	fNumSideP(0),
	23	fNumSideN(0),
	24	fNumLongP(0),
	25	fNumLongN(0),
	26	fDenOutP(0),
	27	fDenOutN(0),
	28	fDenSideP(0),
	29	fDenSideN(0),
	30	fDenLongP(0),
	31	fDenLongN(0)
	32	{
	33	// Default constructor
	34	// set up numerators
	35	char bufname[200];
	36	snprintf(bufname, 200, "NumOutP%s", title);
	37	fNumOutP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
	38	snprintf(bufname, 200, "NumOutN%s", title);
	39	fNumOutN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
	40	snprintf(bufname, 200, "NumSideP%s", title);
	41	fNumSideP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
	42	snprintf(bufname, 200, "NumSideN%s", title);
	43	fNumSideN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
	44	snprintf(bufname, 200, "NumLongP%s", title);
	45	fNumLongP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
	46	snprintf(bufname, 200, "NumLongN%s", title);
	47	fNumLongN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
	48
	49	// set up denominators
	50	snprintf(bufname, 200, "DenOutP%s", title);
	51	fDenOutP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
	52	snprintf(bufname, 200, "DenOutN%s", title);
	53	fDenOutN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
	54	snprintf(bufname, 200, "DenSideP%s", title);
	55	fDenSideP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
	56	snprintf(bufname, 200, "DenSideN%s", title);
	57	fDenSideN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
	58	snprintf(bufname, 200, "DenLongP%s", title);
	59	fDenLongP = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
	60	snprintf(bufname, 200, "DenLongN%s", title);
	61	fDenLongN = new TH1D(bufname,title,nbins,QinvLo,QinvHi);
	62
	63	// to enable error bar calculation...
	64	fNumOutP->Sumw2();
65	fNumOutN->Sumw2();
66	fNumSideP->Sumw2();
67	fNumSideN->Sumw2();
68	fNumLongP->Sumw2();
69	fNumLongN->Sumw2();
70	fDenOutP->Sumw2();
71	fDenOutN->Sumw2();
72	fDenSideP->Sumw2();
73	fDenSideN->Sumw2();
74	fDenLongP->Sumw2();
75	fDenLongN->Sumw2();
76	}
77
78	//____________________________
79	AliFemtoCorrFctnNonIdDR::AliFemtoCorrFctnNonIdDR(const AliFemtoCorrFctnNonIdDR& aCorrFctn) :
80	fNumOutP(0),
81	fNumOutN(0),
82	fNumSideP(0),
83	fNumSideN(0),
84	fNumLongP(0),
85	fNumLongN(0),
86	fDenOutP(0),
87	fDenOutN(0),
88	fDenSideP(0),
89	fDenSideN(0),
90	fDenLongP(0),
91	fDenLongN(0)
92	{
93	// copy constructor
94	if (aCorrFctn.fNumOutP)
95	fNumOutP = new TH1D(*aCorrFctn.fNumOutP);
96	if (aCorrFctn.fNumOutN)
97	fNumOutN = new TH1D(*aCorrFctn.fNumOutN);
98	if (aCorrFctn.fNumSideP)
99	fNumSideP = new TH1D(*aCorrFctn.fNumSideP);
100	if (aCorrFctn.fNumSideN)
101	fNumSideN = new TH1D(*aCorrFctn.fNumSideN);
102	if (aCorrFctn.fNumLongP)
103	fNumLongP = new TH1D(*aCorrFctn.fNumLongP);
104	if (aCorrFctn.fNumLongN)
105	fNumLongN = new TH1D(*aCorrFctn.fNumLongN);
106
107	if (aCorrFctn.fDenOutP)
108	fDenOutP = new TH1D(*aCorrFctn.fDenOutP);
109	if (aCorrFctn.fDenOutN)
110	fDenOutN = new TH1D(*aCorrFctn.fDenOutN);
111	if (aCorrFctn.fDenSideP)
112	fDenSideP = new TH1D(*aCorrFctn.fDenSideP);
113	if (aCorrFctn.fDenSideN)
114	fDenSideN = new TH1D(*aCorrFctn.fDenSideN);
115	if (aCorrFctn.fDenLongP)
116	fDenLongP = new TH1D(*aCorrFctn.fDenLongP);
117	if (aCorrFctn.fDenLongN)
118	fDenLongN = new TH1D(*aCorrFctn.fDenLongN);
119	}
120	//____________________________
121	AliFemtoCorrFctnNonIdDR::~AliFemtoCorrFctnNonIdDR(){
122	delete fNumOutP;
123	delete fNumOutN;
124	delete fNumSideP;
125	delete fNumSideN;
126	delete fNumLongP;
127	delete fNumLongN;
128	delete fDenOutP;
129	delete fDenOutN;
130	delete fDenSideP;
131	delete fDenSideN;
132	delete fDenLongP;
133	delete fDenLongN;
134	}
135	//_________________________
136	AliFemtoCorrFctnNonIdDR& AliFemtoCorrFctnNonIdDR::operator=(const AliFemtoCorrFctnNonIdDR& aCorrFctn)
137	{
138	// assignment operator
139	if (this == &aCorrFctn)
140	return *this;
141
142	if (aCorrFctn.fNumOutP)
143	fNumOutP = new TH1D(*aCorrFctn.fNumOutP);
144	if (aCorrFctn.fNumOutN)
145	fNumOutN = new TH1D(*aCorrFctn.fNumOutN);
146	if (aCorrFctn.fNumSideP)
147	fNumSideP = new TH1D(*aCorrFctn.fNumSideP);
148	if (aCorrFctn.fNumSideN)
149	fNumSideN = new TH1D(*aCorrFctn.fNumSideN);
150	if (aCorrFctn.fNumLongP)
151	fNumLongP = new TH1D(*aCorrFctn.fNumLongP);
152	if (aCorrFctn.fNumLongN)
153	fNumLongN = new TH1D(*aCorrFctn.fNumLongN);
154
155	if (aCorrFctn.fDenOutP)
156	fDenOutP = new TH1D(*aCorrFctn.fDenOutP);
157	if (aCorrFctn.fDenOutN)
158	fDenOutN = new TH1D(*aCorrFctn.fDenOutN);
159	if (aCorrFctn.fDenSideP)
160	fDenSideP = new TH1D(*aCorrFctn.fDenSideP);
161	if (aCorrFctn.fDenSideN)
162	fDenSideN = new TH1D(*aCorrFctn.fDenSideN);
163	if (aCorrFctn.fDenLongP)
164	fDenLongP = new TH1D(*aCorrFctn.fDenLongP);
165	if (aCorrFctn.fDenLongN)
166	fDenLongN = new TH1D(*aCorrFctn.fDenLongN);
167
168	return *this;
169	}
170
171	//_________________________
172	void AliFemtoCorrFctnNonIdDR::Finish(){
173	// here is where we should normalize, fit, etc...
174	// we should NOT Draw() the histos (as I had done it below),
175	// since we want to insulate ourselves from root at this level
176	// of the code. Do it instead at root command line with browser.
177	// fNumerator->Draw();
178	//fDenominator->Draw();
179	//fRatio->Draw();
180	// fRatio->Divide(fNumerator,fDenominator,1.0,1.0);
181
182	}
183
184	//____________________________
185	AliFemtoString AliFemtoCorrFctnNonIdDR::Report(){
186	// construct report
187	string stemp = "Non-identical particles Correlation Function Report:\n";
188	char ctemp[100];
189	sprintf(ctemp,"Number of entries in numerators:\t%E\n",fNumOutP->GetEntries()+fNumOutN->GetEntries());
190	stemp += ctemp;
191	sprintf(ctemp,"Number of entries in denominators:\t%E\n",fDenOutP->GetEntries()+fDenOutN->GetEntries());
192	stemp += ctemp;
193	// stemp += mCoulombWeight->Report();
194	AliFemtoString returnThis = stemp;
195	return returnThis;
196	}
197	//____________________________
198	void AliFemtoCorrFctnNonIdDR::AddRealPair(AliFemtoPair* pair){
199	// add true pair
200	double tKStar = pair->KStar();
201	if (pair->KOut()>0.0)
202	fNumOutP->Fill(tKStar);
203	else
204	fNumOutN->Fill(tKStar);
205
206	if (pair->KSide()>0.0)
207	fNumSideP->Fill(tKStar);
208	else
209	fNumSideN->Fill(tKStar);
210
211	if (pair->KLong()>0.0)
212	fNumLongP->Fill(tKStar);
213	else
214	fNumLongN->Fill(tKStar);
215
216	}
217	//____________________________
218	void AliFemtoCorrFctnNonIdDR::AddMixedPair(AliFemtoPair* pair){
219	// add mixed (background) pair
220	double tKStar = pair->KStar();
221	if (pair->KOut()>0.0)
222	fDenOutP->Fill(tKStar);
223	else
224	fDenOutN->Fill(tKStar);
225
226	if (pair->KSide()>0.0)
227	fDenSideP->Fill(tKStar);
228	else
229	fDenSideN->Fill(tKStar);
230
231	if (pair->KLong()>0.0)
232	fDenLongP->Fill(tKStar);
233	else
234	fDenLongN->Fill(tKStar);
235	}
236	//____________________________
237	void AliFemtoCorrFctnNonIdDR::Write(){
238	fNumOutP->Write();
239	fNumOutN->Write();
240	fNumSideP->Write();
241	fNumSideN->Write();
242	fNumLongP->Write();
243	fNumLongN->Write();
244	fDenOutP->Write();
245	fDenOutN->Write();
246	fDenSideP->Write();
247	fDenSideN->Write();
248	fDenLongP->Write();
249	fDenLongN->Write();
250	}
251
252	TList* AliFemtoCorrFctnNonIdDR::GetOutputList()
253	{
254	// Prepare the list of objects to be written to the output
255	TList *tOutputList = new TList();
256
257	tOutputList->Add(fNumOutP);
258	tOutputList->Add(fNumOutN);
259	tOutputList->Add(fNumSideP);
260	tOutputList->Add(fNumSideN);
261	tOutputList->Add(fNumLongP);
262	tOutputList->Add(fNumLongN);
263	tOutputList->Add(fDenOutP);
264	tOutputList->Add(fDenOutN);
265	tOutputList->Add(fDenSideP);
266	tOutputList->Add(fDenSideN);
267	tOutputList->Add(fDenLongP);
268	tOutputList->Add(fDenLongN);
269
270	return tOutputList;
271	}