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

////////////////////////////////////////////////////////////////////////////////
///                                                                          ///
/// AliFemtoTPCInnerCorrFctn - A correlation function that saves the         ///
/// distance at the entrance to the TPC between two tracks as a function     ///
/// of qinv                                                                  ///
/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
///                                                                          ///
////////////////////////////////////////////////////////////////////////////////

#include "AliFemtoTPCInnerCorrFctn.h"
//#include "AliFemtoHisto.hh"
#include <cstdio>

#ifdef __ROOT__ 
ClassImp(AliFemtoTPCInnerCorrFctn)
#endif

//____________________________
AliFemtoTPCInnerCorrFctn::AliFemtoTPCInnerCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
  fDTPCNumerator(0),
  fDTPCDenominator(0),
  fRadDNumerator(0),
  fRadDDenominator(0),
  fRadius(100)
{
  // set up numerator
  //  title = "Num Qinv (MeV/c)";
  char tTitNum[101] = "NumDTPC";
  strncat(tTitNum,title, 100);
  fDTPCNumerator = new TH2D(tTitNum,title,nbins,QinvLo,QinvHi,100,0.0,20.0);
  // set up denominator
  //title = "Den Qinv (MeV/c)";
  char tTitDen[101] = "DenDTPC";
  strncat(tTitDen,title, 100);
  fDTPCDenominator = new TH2D(tTitDen,title,nbins,QinvLo,QinvHi,100,0.0,20.0);

  char tTitNumR[101] = "NumRadD";
  strncat(tTitNumR,title, 100);
  fRadDNumerator = new TH2D(tTitNumR,title,50,-0.1,0.1,50,-0.1,0.1);
  // set up denominator
  //title = "Den Qinv (MeV/c)";
  char tTitDenR[101] = "DenRadD";
  strncat(tTitDenR,title, 100);
  fRadDDenominator = new TH2D(tTitDenR,title,50,-0.1,0.1,50,-0.1,0.1);

  // to enable error bar calculation...
  fDTPCNumerator->Sumw2();
  fDTPCDenominator->Sumw2();
  fRadDNumerator->Sumw2();
  fRadDDenominator->Sumw2();
}

//____________________________
AliFemtoTPCInnerCorrFctn::AliFemtoTPCInnerCorrFctn(const AliFemtoTPCInnerCorrFctn& aCorrFctn) :
  AliFemtoCorrFctn(),
  fDTPCNumerator(0),
  fDTPCDenominator(0),
  fRadDNumerator(0),
  fRadDDenominator(0),
  fRadius(100)
{
  // copy constructor
  if (aCorrFctn.fDTPCNumerator)
    fDTPCNumerator = new TH2D(*aCorrFctn.fDTPCNumerator);
  if (aCorrFctn.fDTPCDenominator)
    fDTPCDenominator = new TH2D(*aCorrFctn.fDTPCDenominator);
  if (aCorrFctn.fRadDNumerator)
    fRadDNumerator = new TH2D(*aCorrFctn.fRadDNumerator);
  if (aCorrFctn.fRadDDenominator)
    fRadDDenominator = new TH2D(*aCorrFctn.fRadDDenominator);
}
//____________________________
AliFemtoTPCInnerCorrFctn::~AliFemtoTPCInnerCorrFctn(){
  // destructor
  delete fDTPCNumerator;
  delete fDTPCDenominator;
  delete fRadDNumerator;
  delete fRadDDenominator;
}
//_________________________
AliFemtoTPCInnerCorrFctn& AliFemtoTPCInnerCorrFctn::operator=(const AliFemtoTPCInnerCorrFctn& aCorrFctn)
{
  // assignment operator
  if (this == &aCorrFctn)
    return *this;

  if (aCorrFctn.fDTPCNumerator)
    fDTPCNumerator = new TH2D(*aCorrFctn.fDTPCNumerator);
  else
    fDTPCNumerator = 0;
  if (aCorrFctn.fDTPCDenominator)
    fDTPCDenominator = new TH2D(*aCorrFctn.fDTPCDenominator);
  else
    fDTPCDenominator = 0;

  if (aCorrFctn.fRadDNumerator)
    fRadDNumerator = new TH2D(*aCorrFctn.fRadDNumerator);
  else
    fRadDNumerator = 0;
  if (aCorrFctn.fRadDDenominator)
    fRadDDenominator = new TH2D(*aCorrFctn.fRadDDenominator);
  else
    fRadDDenominator = 0;

  fRadius = aCorrFctn.fRadius;

  return *this;
}
//_________________________
void AliFemtoTPCInnerCorrFctn::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 AliFemtoTPCInnerCorrFctn::Report(){
  // create report
  string stemp = "Entrace TPC distance Correlation Function Report:\n";
  char ctemp[100];
  snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fDTPCNumerator->GetEntries());
  stemp += ctemp;
  snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDTPCDenominator->GetEntries());
  stemp += ctemp;
  //  stemp += mCoulombWeight->Report();
  AliFemtoString returnThis = stemp;
  return returnThis;
}
//____________________________
void AliFemtoTPCInnerCorrFctn::AddRealPair( AliFemtoPair* pair){
  // add real (effect) pair
  if (fPairCut)
    if (!fPairCut->Pass(pair)) return;

  double pih = TMath::Pi();
  double pit = TMath::Pi()*2;

  double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
  double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
  double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
  double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
  double dist = sqrt(distx*distx + disty*disty + distz*distz);

  fDTPCNumerator->Fill(tQinv, dist);

  if (tQinv < 0.1) {
    double phi1 = pair->Track1()->Track()->P().Phi();
    double phi2 = pair->Track2()->Track()->P().Phi();
    double chg1 = pair->Track1()->Track()->Charge();
    double chg2 = pair->Track2()->Track()->Charge();
    double ptv1 = pair->Track1()->Track()->Pt();
    double ptv2 = pair->Track2()->Track()->Pt();
    double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
    double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
    double arg1 = -0.3 * 0.5 * chg1 * fRadius/(2*ptv1);
    double arg2 = -0.3 * 0.5 * chg2 * fRadius/(2*ptv2);
    double phid = phi2 - phi1 + TMath::ASin(arg2) - TMath::ASin(arg1);

    while (phid>pih) phid -= pit;
    while (phid<-pih) phid += pit;
    //    dist = phi2 - phi1 + TMath::ASin(-0.3 * 0.5 * chg2 * fRadius/(2*ptv2)) - TMath::ASin(-0.3 * 0.5 * chg1 * fRadius/(2*ptv1));
    double etad = eta2 - eta1;
    fRadDNumerator->Fill(phid, etad);
  }
}
//____________________________
void AliFemtoTPCInnerCorrFctn::AddMixedPair( AliFemtoPair* pair){
  // add mixed (background) pair
  if (fPairCut)
    if (!fPairCut->Pass(pair)) return;

  double pih = TMath::Pi();
  double pit = TMath::Pi()*2;

  double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
  double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
  double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
  double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
  double dist = sqrt(distx*distx + disty*disty + distz*distz);

  fDTPCDenominator->Fill(tQinv,dist);

  if (tQinv < 0.1) {
    double phi1 = pair->Track1()->Track()->P().Phi();
    double phi2 = pair->Track2()->Track()->P().Phi();
    double chg1 = pair->Track1()->Track()->Charge();
    double chg2 = pair->Track2()->Track()->Charge();
    double ptv1 = pair->Track1()->Track()->Pt();
    double ptv2 = pair->Track2()->Track()->Pt();
    double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
    double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
    double arg1 = -0.3 * 0.5 * chg1 * fRadius/(2*ptv1);
    double arg2 = -0.3 * 0.5 * chg2 * fRadius/(2*ptv2);
    double phid = phi2 - phi1 + TMath::ASin(arg2) - TMath::ASin(arg1);

    while (phid>pih) phid -= pit;
    while (phid<-pih) phid += pit;
    //    dist = phi2 - phi1 + TMath::ASin(-0.3 * 0.5 * chg2 * fRadius/(2*ptv2)) - TMath::ASin(-0.3 * 0.5 * chg1 * fRadius/(2*ptv1));
    double etad = eta2 - eta1;
    fRadDDenominator->Fill(phid, etad);
  }
}


void AliFemtoTPCInnerCorrFctn::WriteHistos()
{
  // Write out result histograms
  fDTPCNumerator->Write();
  fDTPCDenominator->Write();
  fRadDNumerator->Write();
  fRadDDenominator->Write();
}
//______________________________
TList* AliFemtoTPCInnerCorrFctn::GetOutputList()
{
  // Prepare the list of objects to be written to the output
  TList *tOutputList = new TList();

  tOutputList->Add(fDTPCNumerator); 
  tOutputList->Add(fDTPCDenominator);  
  tOutputList->Add(fRadDNumerator); 
  tOutputList->Add(fRadDDenominator);  

  return tOutputList;
}

void AliFemtoTPCInnerCorrFctn::SetRadius(double rad)
{
  fRadius = rad;
}
Commit	Line	Data
76ce4b5b	1	////////////////////////////////////////////////////////////////////////////////
	2	/// ///
	3	/// AliFemtoTPCInnerCorrFctn - A correlation function that saves the ///
	4	/// distance at the entrance to the TPC between two tracks as a function ///
	5	/// of qinv ///
	6	/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu ///
	7	/// ///
	8	////////////////////////////////////////////////////////////////////////////////
	9
	10	#include "AliFemtoTPCInnerCorrFctn.h"
	11	//#include "AliFemtoHisto.hh"
	12	#include <cstdio>
	13
	14	#ifdef __ROOT__
	15	ClassImp(AliFemtoTPCInnerCorrFctn)
	16	#endif
	17
	18	//____________________________
	19	AliFemtoTPCInnerCorrFctn::AliFemtoTPCInnerCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
	20	fDTPCNumerator(0),
	21	fDTPCDenominator(0),
	22	fRadDNumerator(0),
	23	fRadDDenominator(0),
	24	fRadius(100)
	25	{
	26	// set up numerator
	27	// title = "Num Qinv (MeV/c)";
	28	char tTitNum[101] = "NumDTPC";
	29	strncat(tTitNum,title, 100);
	30	fDTPCNumerator = new TH2D(tTitNum,title,nbins,QinvLo,QinvHi,100,0.0,20.0);
	31	// set up denominator
	32	//title = "Den Qinv (MeV/c)";
	33	char tTitDen[101] = "DenDTPC";
	34	strncat(tTitDen,title, 100);
	35	fDTPCDenominator = new TH2D(tTitDen,title,nbins,QinvLo,QinvHi,100,0.0,20.0);
	36
	37	char tTitNumR[101] = "NumRadD";
	38	strncat(tTitNumR,title, 100);
	39	fRadDNumerator = new TH2D(tTitNumR,title,50,-0.1,0.1,50,-0.1,0.1);
	40	// set up denominator
	41	//title = "Den Qinv (MeV/c)";
	42	char tTitDenR[101] = "DenRadD";
	43	strncat(tTitDenR,title, 100);
	44	fRadDDenominator = new TH2D(tTitDenR,title,50,-0.1,0.1,50,-0.1,0.1);
	45
	46	// to enable error bar calculation...
	47	fDTPCNumerator->Sumw2();
	48	fDTPCDenominator->Sumw2();
	49	fRadDNumerator->Sumw2();
	50	fRadDDenominator->Sumw2();
	51	}
	52
	53	//____________________________
	54	AliFemtoTPCInnerCorrFctn::AliFemtoTPCInnerCorrFctn(const AliFemtoTPCInnerCorrFctn& aCorrFctn) :
	55	AliFemtoCorrFctn(),
	56	fDTPCNumerator(0),
	57	fDTPCDenominator(0),
	58	fRadDNumerator(0),
	59	fRadDDenominator(0),
	60	fRadius(100)
	61	{
	62	// copy constructor
	63	if (aCorrFctn.fDTPCNumerator)
	64	fDTPCNumerator = new TH2D(*aCorrFctn.fDTPCNumerator);
65	if (aCorrFctn.fDTPCDenominator)
66	fDTPCDenominator = new TH2D(*aCorrFctn.fDTPCDenominator);
67	if (aCorrFctn.fRadDNumerator)
68	fRadDNumerator = new TH2D(*aCorrFctn.fRadDNumerator);
69	if (aCorrFctn.fRadDDenominator)
70	fRadDDenominator = new TH2D(*aCorrFctn.fRadDDenominator);
71	}
72	//____________________________
73	AliFemtoTPCInnerCorrFctn::~AliFemtoTPCInnerCorrFctn(){
74	// destructor
75	delete fDTPCNumerator;
76	delete fDTPCDenominator;
77	delete fRadDNumerator;
78	delete fRadDDenominator;
79	}
80	//_________________________
81	AliFemtoTPCInnerCorrFctn& AliFemtoTPCInnerCorrFctn::operator=(const AliFemtoTPCInnerCorrFctn& aCorrFctn)
82	{
83	// assignment operator
84	if (this == &aCorrFctn)
85	return *this;
86
87	if (aCorrFctn.fDTPCNumerator)
88	fDTPCNumerator = new TH2D(*aCorrFctn.fDTPCNumerator);
89	else
90	fDTPCNumerator = 0;
91	if (aCorrFctn.fDTPCDenominator)
92	fDTPCDenominator = new TH2D(*aCorrFctn.fDTPCDenominator);
93	else
94	fDTPCDenominator = 0;
95
96	if (aCorrFctn.fRadDNumerator)
97	fRadDNumerator = new TH2D(*aCorrFctn.fRadDNumerator);
98	else
99	fRadDNumerator = 0;
100	if (aCorrFctn.fRadDDenominator)
101	fRadDDenominator = new TH2D(*aCorrFctn.fRadDDenominator);
102	else
103	fRadDDenominator = 0;
104
105	fRadius = aCorrFctn.fRadius;
106
107	return *this;
108	}
109	//_________________________
110	void AliFemtoTPCInnerCorrFctn::Finish(){
111	// here is where we should normalize, fit, etc...
112	// we should NOT Draw() the histos (as I had done it below),
113	// since we want to insulate ourselves from root at this level
114	// of the code. Do it instead at root command line with browser.
115	// mShareNumerator->Draw();
116	//mShareDenominator->Draw();
117	//mRatio->Draw();
118
119	}
120
121	//____________________________
122	AliFemtoString AliFemtoTPCInnerCorrFctn::Report(){
123	// create report
124	string stemp = "Entrace TPC distance Correlation Function Report:\n";
125	char ctemp[100];
126	snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fDTPCNumerator->GetEntries());
127	stemp += ctemp;
128	snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDTPCDenominator->GetEntries());
129	stemp += ctemp;
130	// stemp += mCoulombWeight->Report();
131	AliFemtoString returnThis = stemp;
132	return returnThis;
133	}
134	//____________________________
135	void AliFemtoTPCInnerCorrFctn::AddRealPair( AliFemtoPair* pair){
136	// add real (effect) pair
137	if (fPairCut)
138	if (!fPairCut->Pass(pair)) return;
139
140	double pih = TMath::Pi();
141	double pit = TMath::Pi()*2;
142
143	double tQinv = fabs(pair->QInv()); // note - qInv() will be negative for identical pairs...
144	double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
145	double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
146	double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
147	double dist = sqrt(distxdistx + distydisty + distz*distz);
148
149	fDTPCNumerator->Fill(tQinv, dist);
150
151	if (tQinv < 0.1) {
152	double phi1 = pair->Track1()->Track()->P().Phi();
153	double phi2 = pair->Track2()->Track()->P().Phi();
154	double chg1 = pair->Track1()->Track()->Charge();
155	double chg2 = pair->Track2()->Track()->Charge();
156	double ptv1 = pair->Track1()->Track()->Pt();
157	double ptv2 = pair->Track2()->Track()->Pt();
158	double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
159	double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
160	double arg1 = -0.3 * 0.5 * chg1 * fRadius/(2*ptv1);
161	double arg2 = -0.3 * 0.5 * chg2 * fRadius/(2*ptv2);
162	double phid = phi2 - phi1 + TMath::ASin(arg2) - TMath::ASin(arg1);
163
164	while (phid>pih) phid -= pit;
165	while (phid<-pih) phid += pit;
166	// dist = phi2 - phi1 + TMath::ASin(-0.3 * 0.5 * chg2 * fRadius/(2ptv2)) - TMath::ASin(-0.3 0.5 * chg1 * fRadius/(2*ptv1));
167	double etad = eta2 - eta1;
168	fRadDNumerator->Fill(phid, etad);
169	}
170	}
171	//____________________________
172	void AliFemtoTPCInnerCorrFctn::AddMixedPair( AliFemtoPair* pair){
173	// add mixed (background) pair
174	if (fPairCut)
175	if (!fPairCut->Pass(pair)) return;
176
177	double pih = TMath::Pi();
178	double pit = TMath::Pi()*2;
179
180	double tQinv = fabs(pair->QInv()); // note - qInv() will be negative for identical pairs...
181	double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
182	double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
183	double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
184	double dist = sqrt(distxdistx + distydisty + distz*distz);
185
186	fDTPCDenominator->Fill(tQinv,dist);
187
188	if (tQinv < 0.1) {
189	double phi1 = pair->Track1()->Track()->P().Phi();
190	double phi2 = pair->Track2()->Track()->P().Phi();
191	double chg1 = pair->Track1()->Track()->Charge();
192	double chg2 = pair->Track2()->Track()->Charge();
193	double ptv1 = pair->Track1()->Track()->Pt();
194	double ptv2 = pair->Track2()->Track()->Pt();
195	double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
196	double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
197	double arg1 = -0.3 * 0.5 * chg1 * fRadius/(2*ptv1);
198	double arg2 = -0.3 * 0.5 * chg2 * fRadius/(2*ptv2);
199	double phid = phi2 - phi1 + TMath::ASin(arg2) - TMath::ASin(arg1);
200
201	while (phid>pih) phid -= pit;
202	while (phid<-pih) phid += pit;
203	// dist = phi2 - phi1 + TMath::ASin(-0.3 * 0.5 * chg2 * fRadius/(2ptv2)) - TMath::ASin(-0.3 0.5 * chg1 * fRadius/(2*ptv1));
204	double etad = eta2 - eta1;
205	fRadDDenominator->Fill(phid, etad);
206	}
207	}
208
209
210	void AliFemtoTPCInnerCorrFctn::WriteHistos()
211	{
212	// Write out result histograms
213	fDTPCNumerator->Write();
214	fDTPCDenominator->Write();
215	fRadDNumerator->Write();
216	fRadDDenominator->Write();
217	}
218	//______________________________
219	TList* AliFemtoTPCInnerCorrFctn::GetOutputList()
220	{
221	// Prepare the list of objects to be written to the output
222	TList *tOutputList = new TList();
223
224	tOutputList->Add(fDTPCNumerator);
225	tOutputList->Add(fDTPCDenominator);
226	tOutputList->Add(fRadDNumerator);
227	tOutputList->Add(fRadDDenominator);
228
229	return tOutputList;
230	}
231
232	void AliFemtoTPCInnerCorrFctn::SetRadius(double rad)
233	{
234	fRadius = rad;
235	}