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

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
// AliFemtoGammaMonitor - A correlation function that analyzes            //
// two particle correlations with respect to the azimuthal angle (phi)        //
// and pseudorapidity (eta) difference                                        //
//                                                                            //
// Authors: Adam Kisiel Adam.Kisiel@cern.ch                                   //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

#include "AliFemtoCorrFctnGammaMonitor.h"
//#include "AliFemtoHisto.hh"
#include <cstdio>
#include <TMath.h>

#ifdef __ROOT__ 
ClassImp(AliFemtoCorrFctnGammaMonitor)
#endif

//____________________________
AliFemtoCorrFctnGammaMonitor::AliFemtoCorrFctnGammaMonitor(char* title, const int& aMinvBins=20, const int& aDThetaBins=20):
  AliFemtoCorrFctn(),
  fNumPMinvDTheta(0),
  fDenPMinvDTheta(0),
  fNumNMinvDTheta(0),
  fDenNMinvDTheta(0)
{
  // set up numerator
  char tTitNumD[101] = "NumPMinvTheta";
  strncat(tTitNumD,title, 100);
  fNumPMinvDTheta = new TH2D(tTitNumD,title,aMinvBins,0.0,0.2,aDThetaBins,0.0,0.2);
  // set up denominator
  char tTitDenD[101] = "DenPMinvTheta";
  strncat(tTitDenD,title, 100);
  fDenPMinvDTheta = new TH2D(tTitDenD,title,aMinvBins,0.0,0.2,aDThetaBins,0.0,0.2);

  // set up numerator
  char tTitNumR[101] = "NumNMinvTheta";
  strncat(tTitNumR,title, 100);
  fNumNMinvDTheta = new TH2D(tTitNumR,title,aMinvBins,0.0,0.2,aDThetaBins,0.0,0.2);
  // set up denominator
  char tTitDenR[101] = "DenNMinvTheta";
  strncat(tTitDenR,title, 100);
  fDenNMinvDTheta = new TH2D(tTitDenR,title,aMinvBins,0.0,0.2,aDThetaBins,0.0,0.2);

  // to enable error bar calculation...
  fNumPMinvDTheta->Sumw2();
  fDenPMinvDTheta->Sumw2();
  fNumNMinvDTheta->Sumw2();
  fDenNMinvDTheta->Sumw2();
}

//____________________________
AliFemtoCorrFctnGammaMonitor::AliFemtoCorrFctnGammaMonitor(const AliFemtoCorrFctnGammaMonitor& aCorrFctn) :
  AliFemtoCorrFctn(),
  fNumPMinvDTheta(0),
  fDenPMinvDTheta(0),
  fNumNMinvDTheta(0),
  fDenNMinvDTheta(0)
{
  // copy constructor
  if (aCorrFctn.fNumPMinvDTheta)
    fNumPMinvDTheta = new TH2D(*aCorrFctn.fNumPMinvDTheta);
  else
    fNumPMinvDTheta = 0;
  if (aCorrFctn.fDenPMinvDTheta)
    fDenPMinvDTheta = new TH2D(*aCorrFctn.fDenPMinvDTheta);
  else
    fDenPMinvDTheta = 0;

  if (aCorrFctn.fNumNMinvDTheta)
    fNumNMinvDTheta = new TH2D(*aCorrFctn.fNumNMinvDTheta);
  else
    fNumNMinvDTheta = 0;
  if (aCorrFctn.fDenNMinvDTheta)
    fDenNMinvDTheta = new TH2D(*aCorrFctn.fDenNMinvDTheta);
  else
    fDenNMinvDTheta = 0;

}
//____________________________
AliFemtoCorrFctnGammaMonitor::~AliFemtoCorrFctnGammaMonitor(){
  // destructor
  delete fNumPMinvDTheta;
  delete fDenPMinvDTheta;
  delete fNumNMinvDTheta;
  delete fDenNMinvDTheta;
}
//_________________________
AliFemtoCorrFctnGammaMonitor& AliFemtoCorrFctnGammaMonitor::operator=(const AliFemtoCorrFctnGammaMonitor& aCorrFctn)
{
  // assignment operator
  if (this == &aCorrFctn)
    return *this;

  if (aCorrFctn.fNumPMinvDTheta)
    fNumPMinvDTheta = new TH2D(*aCorrFctn.fNumPMinvDTheta);
  else
    fNumPMinvDTheta = 0;
  if (aCorrFctn.fDenPMinvDTheta)
    fDenPMinvDTheta = new TH2D(*aCorrFctn.fDenPMinvDTheta);
  else
    fDenPMinvDTheta = 0;

  if (aCorrFctn.fNumNMinvDTheta)
    fNumNMinvDTheta = new TH2D(*aCorrFctn.fNumNMinvDTheta);
  else
    fNumNMinvDTheta = 0;
  if (aCorrFctn.fDenNMinvDTheta)
    fDenNMinvDTheta = new TH2D(*aCorrFctn.fDenNMinvDTheta);
  else
    fDenNMinvDTheta = 0;


  return *this;
}
//_________________________
void AliFemtoCorrFctnGammaMonitor::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.
  //  mShareNumerator->Draw();
  //mShareDenominator->Draw();
  //mRatio->Draw();

}

//____________________________
AliFemtoString AliFemtoCorrFctnGammaMonitor::Report(){
  // create report
  string stemp = "Gamma Monitor Function Report:\n";
  char ctemp[100];
  snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fNumPMinvDTheta->GetEntries());
  stemp += ctemp;
  snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDenPMinvDTheta->GetEntries());
  stemp += ctemp;
  //  stemp += mCoulombWeight->Report();
  AliFemtoString returnThis = stemp;
  return returnThis;
}
//____________________________
void AliFemtoCorrFctnGammaMonitor::AddRealPair( AliFemtoPair* pair){
  // add real (effect) pair
  double me = 0.000511;

  double theta1 = pair->Track1()->Track()->P().Theta();
  double theta2 = pair->Track2()->Track()->P().Theta();
  double dtheta = TMath::Abs(theta1 - theta2);

  double e1 = TMath::Sqrt(me*me + pair->Track1()->Track()->P().Mag2());
  double e2 = TMath::Sqrt(me*me + pair->Track2()->Track()->P().Mag2());

  double minv = 2*me*me + 2*(e1*e2 - 
			     pair->Track1()->Track()->P().x()*pair->Track2()->Track()->P().x() -
			     pair->Track1()->Track()->P().y()*pair->Track2()->Track()->P().y() -
			     pair->Track1()->Track()->P().z()*pair->Track2()->Track()->P().z());

  if (pair->KSide()>0.0) 
    fNumPMinvDTheta->Fill(minv, dtheta);
  else
    fNumNMinvDTheta->Fill(minv, dtheta);
}
//____________________________
void AliFemtoCorrFctnGammaMonitor::AddMixedPair( AliFemtoPair* pair){
  // add mixed (background) pair
  double me = 0.000511;

  double theta1 = pair->Track1()->Track()->P().Theta();
  double theta2 = pair->Track2()->Track()->P().Theta();
  double dtheta = TMath::Abs(theta1 - theta2);

  double e1 = TMath::Sqrt(me*me + pair->Track1()->Track()->P().Mag2());
  double e2 = TMath::Sqrt(me*me + pair->Track2()->Track()->P().Mag2());

  double minv = 2*me*me + 2*(e1*e2 - 
			     pair->Track1()->Track()->P().x()*pair->Track2()->Track()->P().x() -
			     pair->Track1()->Track()->P().y()*pair->Track2()->Track()->P().y() -
			     pair->Track1()->Track()->P().z()*pair->Track2()->Track()->P().z());

  if (pair->KSide()>0.0) 
    fDenPMinvDTheta->Fill(minv, dtheta);
  else
    fDenNMinvDTheta->Fill(minv, dtheta);
}


void AliFemtoCorrFctnGammaMonitor::WriteHistos()
{
  // Write out result histograms
  fNumPMinvDTheta->Write();
  fDenPMinvDTheta->Write();
  fNumNMinvDTheta->Write();
  fDenNMinvDTheta->Write();
}

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

  tOutputList->Add(fNumPMinvDTheta);
  tOutputList->Add(fDenPMinvDTheta);
  tOutputList->Add(fNumNMinvDTheta);
  tOutputList->Add(fDenNMinvDTheta);

  return tOutputList;
}
Commit	Line	Data
76ce4b5b	1	////////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoGammaMonitor - A correlation function that analyzes //
	4	// two particle correlations with respect to the azimuthal angle (phi) //
	5	// and pseudorapidity (eta) difference //
	6	// //
	7	// Authors: Adam Kisiel Adam.Kisiel@cern.ch //
	8	// //
	9	////////////////////////////////////////////////////////////////////////////////
	10
	11	#include "AliFemtoCorrFctnGammaMonitor.h"
	12	//#include "AliFemtoHisto.hh"
	13	#include <cstdio>
	14	#include <TMath.h>
	15
	16	#ifdef __ROOT__
	17	ClassImp(AliFemtoCorrFctnGammaMonitor)
	18	#endif
	19
	20	//____________________________
	21	AliFemtoCorrFctnGammaMonitor::AliFemtoCorrFctnGammaMonitor(char* title, const int& aMinvBins=20, const int& aDThetaBins=20):
	22	AliFemtoCorrFctn(),
	23	fNumPMinvDTheta(0),
	24	fDenPMinvDTheta(0),
	25	fNumNMinvDTheta(0),
	26	fDenNMinvDTheta(0)
	27	{
	28	// set up numerator
	29	char tTitNumD[101] = "NumPMinvTheta";
	30	strncat(tTitNumD,title, 100);
	31	fNumPMinvDTheta = new TH2D(tTitNumD,title,aMinvBins,0.0,0.2,aDThetaBins,0.0,0.2);
	32	// set up denominator
	33	char tTitDenD[101] = "DenPMinvTheta";
	34	strncat(tTitDenD,title, 100);
	35	fDenPMinvDTheta = new TH2D(tTitDenD,title,aMinvBins,0.0,0.2,aDThetaBins,0.0,0.2);
	36
	37	// set up numerator
	38	char tTitNumR[101] = "NumNMinvTheta";
	39	strncat(tTitNumR,title, 100);
	40	fNumNMinvDTheta = new TH2D(tTitNumR,title,aMinvBins,0.0,0.2,aDThetaBins,0.0,0.2);
	41	// set up denominator
	42	char tTitDenR[101] = "DenNMinvTheta";
	43	strncat(tTitDenR,title, 100);
	44	fDenNMinvDTheta = new TH2D(tTitDenR,title,aMinvBins,0.0,0.2,aDThetaBins,0.0,0.2);
	45
	46	// to enable error bar calculation...
	47	fNumPMinvDTheta->Sumw2();
	48	fDenPMinvDTheta->Sumw2();
	49	fNumNMinvDTheta->Sumw2();
	50	fDenNMinvDTheta->Sumw2();
	51	}
	52
	53	//____________________________
	54	AliFemtoCorrFctnGammaMonitor::AliFemtoCorrFctnGammaMonitor(const AliFemtoCorrFctnGammaMonitor& aCorrFctn) :
	55	AliFemtoCorrFctn(),
	56	fNumPMinvDTheta(0),
	57	fDenPMinvDTheta(0),
	58	fNumNMinvDTheta(0),
	59	fDenNMinvDTheta(0)
	60	{
	61	// copy constructor
	62	if (aCorrFctn.fNumPMinvDTheta)
	63	fNumPMinvDTheta = new TH2D(*aCorrFctn.fNumPMinvDTheta);
	64	else
65	fNumPMinvDTheta = 0;
66	if (aCorrFctn.fDenPMinvDTheta)
67	fDenPMinvDTheta = new TH2D(*aCorrFctn.fDenPMinvDTheta);
68	else
69	fDenPMinvDTheta = 0;
70
71	if (aCorrFctn.fNumNMinvDTheta)
72	fNumNMinvDTheta = new TH2D(*aCorrFctn.fNumNMinvDTheta);
73	else
74	fNumNMinvDTheta = 0;
75	if (aCorrFctn.fDenNMinvDTheta)
76	fDenNMinvDTheta = new TH2D(*aCorrFctn.fDenNMinvDTheta);
77	else
78	fDenNMinvDTheta = 0;
79
80	}
81	//____________________________
82	AliFemtoCorrFctnGammaMonitor::~AliFemtoCorrFctnGammaMonitor(){
83	// destructor
84	delete fNumPMinvDTheta;
85	delete fDenPMinvDTheta;
86	delete fNumNMinvDTheta;
87	delete fDenNMinvDTheta;
88	}
89	//_________________________
90	AliFemtoCorrFctnGammaMonitor& AliFemtoCorrFctnGammaMonitor::operator=(const AliFemtoCorrFctnGammaMonitor& aCorrFctn)
91	{
92	// assignment operator
93	if (this == &aCorrFctn)
94	return *this;
95
96	if (aCorrFctn.fNumPMinvDTheta)
97	fNumPMinvDTheta = new TH2D(*aCorrFctn.fNumPMinvDTheta);
98	else
99	fNumPMinvDTheta = 0;
100	if (aCorrFctn.fDenPMinvDTheta)
101	fDenPMinvDTheta = new TH2D(*aCorrFctn.fDenPMinvDTheta);
102	else
103	fDenPMinvDTheta = 0;
104
105	if (aCorrFctn.fNumNMinvDTheta)
106	fNumNMinvDTheta = new TH2D(*aCorrFctn.fNumNMinvDTheta);
107	else
108	fNumNMinvDTheta = 0;
109	if (aCorrFctn.fDenNMinvDTheta)
110	fDenNMinvDTheta = new TH2D(*aCorrFctn.fDenNMinvDTheta);
111	else
112	fDenNMinvDTheta = 0;
113
114
115	return *this;
116	}
117	//_________________________
118	void AliFemtoCorrFctnGammaMonitor::Finish(){
119	// here is where we should normalize, fit, etc...
120	// we should NOT Draw() the histos (as I had done it below),
121	// since we want to insulate ourselves from root at this level
122	// of the code. Do it instead at root command line with browser.
123	// mShareNumerator->Draw();
124	//mShareDenominator->Draw();
125	//mRatio->Draw();
126
127	}
128
129	//____________________________
130	AliFemtoString AliFemtoCorrFctnGammaMonitor::Report(){
131	// create report
132	string stemp = "Gamma Monitor Function Report:\n";
133	char ctemp[100];
134	snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fNumPMinvDTheta->GetEntries());
135	stemp += ctemp;
136	snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDenPMinvDTheta->GetEntries());
137	stemp += ctemp;
138	// stemp += mCoulombWeight->Report();
139	AliFemtoString returnThis = stemp;
140	return returnThis;
141	}
142	//____________________________
143	void AliFemtoCorrFctnGammaMonitor::AddRealPair( AliFemtoPair* pair){
144	// add real (effect) pair
145	double me = 0.000511;
146
147	double theta1 = pair->Track1()->Track()->P().Theta();
148	double theta2 = pair->Track2()->Track()->P().Theta();
149	double dtheta = TMath::Abs(theta1 - theta2);
150
151	double e1 = TMath::Sqrt(me*me + pair->Track1()->Track()->P().Mag2());
152	double e2 = TMath::Sqrt(me*me + pair->Track2()->Track()->P().Mag2());
153
154	double minv = 2meme + 2(e1e2 -
155	pair->Track1()->Track()->P().x()*pair->Track2()->Track()->P().x() -
156	pair->Track1()->Track()->P().y()*pair->Track2()->Track()->P().y() -
157	pair->Track1()->Track()->P().z()*pair->Track2()->Track()->P().z());
158
159	if (pair->KSide()>0.0)
160	fNumPMinvDTheta->Fill(minv, dtheta);
161	else
162	fNumNMinvDTheta->Fill(minv, dtheta);
163	}
164	//____________________________
165	void AliFemtoCorrFctnGammaMonitor::AddMixedPair( AliFemtoPair* pair){
166	// add mixed (background) pair
167	double me = 0.000511;
168
169	double theta1 = pair->Track1()->Track()->P().Theta();
170	double theta2 = pair->Track2()->Track()->P().Theta();
171	double dtheta = TMath::Abs(theta1 - theta2);
172
173	double e1 = TMath::Sqrt(me*me + pair->Track1()->Track()->P().Mag2());
174	double e2 = TMath::Sqrt(me*me + pair->Track2()->Track()->P().Mag2());
175
176	double minv = 2meme + 2(e1e2 -
177	pair->Track1()->Track()->P().x()*pair->Track2()->Track()->P().x() -
178	pair->Track1()->Track()->P().y()*pair->Track2()->Track()->P().y() -
179	pair->Track1()->Track()->P().z()*pair->Track2()->Track()->P().z());
180
181	if (pair->KSide()>0.0)
182	fDenPMinvDTheta->Fill(minv, dtheta);
183	else
184	fDenNMinvDTheta->Fill(minv, dtheta);
185	}
186
187
188	void AliFemtoCorrFctnGammaMonitor::WriteHistos()
189	{
190	// Write out result histograms
191	fNumPMinvDTheta->Write();
192	fDenPMinvDTheta->Write();
193	fNumNMinvDTheta->Write();
194	fDenNMinvDTheta->Write();
195	}
196
197	TList* AliFemtoCorrFctnGammaMonitor::GetOutputList()
198	{
199	// Prepare the list of objects to be written to the output
200	TList *tOutputList = new TList();
201
202	tOutputList->Add(fNumPMinvDTheta);
203	tOutputList->Add(fDenPMinvDTheta);
204	tOutputList->Add(fNumNMinvDTheta);
205	tOutputList->Add(fDenNMinvDTheta);
206
207	return tOutputList;
208	}