[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / CorrFctn / AliFemtoShareQualityCorrFctn.cxx

////////////////////////////////////////////////////////////////////////////////
///                                                                          ///
/// AliFemtoShareQualityCorrFctn - A correlation function that saves the     ///
/// amount of sharing and splitting hits per pair as a function of qinv      ///
/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
///                                                                          ///
////////////////////////////////////////////////////////////////////////////////

#include "AliFemtoShareQualityCorrFctn.h"
//#include "Infrastructure/AliFemtoHisto.hh"
#include <cstdio>

#ifdef __ROOT__ 
ClassImp(AliFemtoShareQualityCorrFctn)
#endif

//____________________________
AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
  fShareNumerator(0),
  fShareDenominator(0),
  fQualityNumerator(0),
  fQualityDenominator(0)
{
  // set up numerator
  //  title = "Num Qinv (MeV/c)";
  char TitNum[100] = "NumShare";
  strcat(TitNum,title);
  fShareNumerator = new TH2D(TitNum,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
  // set up denominator
  //title = "Den Qinv (MeV/c)";
  char TitDen[100] = "DenShare";
  strcat(TitDen,title);
  fShareDenominator = new TH2D(TitDen,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);

  char Tit2Num[100] = "NumQuality";
  strcat(Tit2Num,title);
  fQualityNumerator = new TH2D(Tit2Num,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
  // set up denominator
  //title = "Den Qinv (MeV/c)";
  char Tit2Den[100] = "DenQuality";
  strcat(Tit2Den,title);
  fQualityDenominator = new TH2D(Tit2Den,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
  // set up ratio
  //title = "Ratio Qinv (MeV/c)";
  // this next bit is unfortunately needed so that we can have many histos of same "title"
  // it is neccessary if we typedef TH2D to TH1d (which we do)
  //mShareNumerator->SetDirectory(0);
  //mShareDenominator->SetDirectory(0);
  //mRatio->SetDirectory(0);

  // to enable error bar calculation...
  fShareNumerator->Sumw2();
  fShareDenominator->Sumw2();

  fQualityNumerator->Sumw2();
  fQualityDenominator->Sumw2();
}

//____________________________
AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(const AliFemtoShareQualityCorrFctn& aCorrFctn) :
  fShareNumerator(0),
  fShareDenominator(0),
  fQualityNumerator(0),
  fQualityDenominator(0)
{
  if (aCorrFctn.fShareNumerator)
    fShareNumerator = new TH2D(*aCorrFctn.fShareNumerator);
  if (aCorrFctn.fShareDenominator)
    fShareDenominator = new TH2D(*aCorrFctn.fShareDenominator);
  if (aCorrFctn.fQualityNumerator)
    fQualityNumerator = new TH2D(*aCorrFctn.fQualityNumerator);
  if (aCorrFctn.fQualityDenominator)
    fQualityDenominator = new TH2D(*aCorrFctn.fQualityDenominator);
}
//____________________________
AliFemtoShareQualityCorrFctn::~AliFemtoShareQualityCorrFctn(){
  delete fShareNumerator;
  delete fShareDenominator;
  delete fQualityNumerator;
  delete fQualityDenominator;
}
//_________________________
AliFemtoShareQualityCorrFctn& AliFemtoShareQualityCorrFctn::operator=(const AliFemtoShareQualityCorrFctn& aCorrFctn)
{
  if (this == &aCorrFctn)
    return *this;

  if (aCorrFctn.fShareNumerator)
    fShareNumerator = new TH2D(*aCorrFctn.fShareNumerator);
  else
    fShareNumerator = 0;
  if (aCorrFctn.fShareDenominator)
    fShareDenominator = new TH2D(*aCorrFctn.fShareDenominator);
  else
    fShareDenominator = 0;
  if (aCorrFctn.fQualityNumerator)
    fQualityNumerator = new TH2D(*aCorrFctn.fQualityNumerator);
  else
    fQualityNumerator = 0;
  if (aCorrFctn.fQualityDenominator)
    fQualityDenominator = new TH2D(*aCorrFctn.fQualityDenominator);
  else
    fQualityDenominator = 0;

  return *this;
}
//_________________________
void AliFemtoShareQualityCorrFctn::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 AliFemtoShareQualityCorrFctn::Report(){
  string stemp = "Qinv Correlation Function Report:\n";
  char ctemp[100];
  sprintf(ctemp,"Number of entries in numerator:\t%E\n",fShareNumerator->GetEntries());
  stemp += ctemp;
  sprintf(ctemp,"Number of entries in denominator:\t%E\n",fShareDenominator->GetEntries());
  stemp += ctemp;
  //  stemp += mCoulombWeight->Report();
  AliFemtoString returnThis = stemp;
  return returnThis;
}
//____________________________
void AliFemtoShareQualityCorrFctn::AddRealPair(const AliFemtoPair* pair){
  double Qinv = fabs(pair->qInv());   // note - qInv() will be negative for identical pairs...
  Int_t nh = 0;
  Int_t an = 0;
  Int_t ns = 0;
  
  for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
    // If both have clusters in the same row
    if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) && 
	pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
      // Do they share it ?
      if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap) ||
	  pair->track2()->Track()->TPCsharing().TestBitNumber(imap))
	{
	  if (Qinv < 0.01) {
	    cout << "Shared cluster in row " << imap << endl; 
	  }
	  an++;
	  nh+=2;
	  ns+=2;
	}
      
      // Different hits on the same padrow
      else {
	an--;
	nh+=2;
      }
    }
    else if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) ||
	     pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
      // One track has a hit, the other does not
      an++;
      nh++;
    }
  }
  if (Qinv < 0.01) {
    cout << "Qinv of the pair is " << Qinv << endl;
    cout << "Clusters: " << endl;
    for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
      cout << imap ;
      if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
      else cout << " 0 " ;
      if (pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
      else cout << " 0 " ;
      cout << "     ";
      if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
      else cout << " X ";
      if (pair->track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
      else cout << " X ";
      cout << endl;
    }
  }

  Float_t hsmval = 0.0;
  Float_t hsfval = 0.0;

  if (nh >0) {
    hsmval = an*1.0/nh;
    hsfval = ns*1.0/nh;
  }

  if (Qinv < 0.01) {
    cout << "Quality  Sharity " << hsmval << " " << hsfval << " " << pair->track1()->Track() << " " << pair->track2()->Track() << endl;
  }

  fShareNumerator->Fill(Qinv, hsfval);
  fQualityNumerator->Fill(Qinv, hsmval);
  //  cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << Qinv <<
  //" " << pair->track1().FourMomentum() << " " << pair->track2().FourMomentum() << endl;
}
//____________________________
void AliFemtoShareQualityCorrFctn::AddMixedPair(const AliFemtoPair* pair){
  double weight = 1.0;
  double Qinv = fabs(pair->qInv());   // note - qInv() will be negative for identical pairs...
  Int_t nh = 0;
  Int_t an = 0;
  Int_t ns = 0;
  
  for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
    // If both have clusters in the same row
    if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) && 
	pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
      // Do they share it ?
      if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap) ||
	  pair->track2()->Track()->TPCsharing().TestBitNumber(imap))
	{
	  //	  cout << "A shared cluster !!!" << endl;
	  //	cout << "imap idx1 idx2 " 
	  //	     << imap << " "
	  //	     << tP1idx[imap] << " " << tP2idx[imap] << endl;
	  an++;
	  nh+=2;
	  ns+=2;
	}
      
      // Different hits on the same padrow
      else {
	an--;
	nh+=2;
      }
    }
    else if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) ||
	     pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
      // One track has a hit, the other does not
      an++;
      nh++;
    }
  }
  
  Float_t hsmval = 0.0;
  Float_t hsfval = 0.0;

  if (nh >0) {
    hsmval = an*1.0/nh;
    hsfval = ns*1.0/nh;
  }

  fShareDenominator->Fill(Qinv,hsfval,weight);
  fQualityDenominator->Fill(Qinv,hsmval,weight);
}


void AliFemtoShareQualityCorrFctn::WriteHistos()
{
  fShareNumerator->Write();
  fShareDenominator->Write();
  fQualityNumerator->Write();
  fQualityDenominator->Write();
  
}
Commit	Line	Data
0215f606	1	////////////////////////////////////////////////////////////////////////////////
	2	/// ///
	3	/// AliFemtoShareQualityCorrFctn - A correlation function that saves the ///
	4	/// amount of sharing and splitting hits per pair as a function of qinv ///
	5	/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu ///
	6	/// ///
	7	////////////////////////////////////////////////////////////////////////////////
67427ff7	8
b2f60a91	9	#include "AliFemtoShareQualityCorrFctn.h"
67427ff7	10	//#include "Infrastructure/AliFemtoHisto.hh"
	11	#include <cstdio>
	12
	13	#ifdef __ROOT__
	14	ClassImp(AliFemtoShareQualityCorrFctn)
	15	#endif
	16
	17	//____________________________
0215f606	18	AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
	19	fShareNumerator(0),
	20	fShareDenominator(0),
	21	fQualityNumerator(0),
	22	fQualityDenominator(0)
	23	{
67427ff7	24	// set up numerator
	25	// title = "Num Qinv (MeV/c)";
	26	char TitNum[100] = "NumShare";
	27	strcat(TitNum,title);
	28	fShareNumerator = new TH2D(TitNum,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
	29	// set up denominator
	30	//title = "Den Qinv (MeV/c)";
	31	char TitDen[100] = "DenShare";
	32	strcat(TitDen,title);
	33	fShareDenominator = new TH2D(TitDen,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
	34
	35	char Tit2Num[100] = "NumQuality";
	36	strcat(Tit2Num,title);
	37	fQualityNumerator = new TH2D(Tit2Num,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
	38	// set up denominator
	39	//title = "Den Qinv (MeV/c)";
	40	char Tit2Den[100] = "DenQuality";
	41	strcat(Tit2Den,title);
	42	fQualityDenominator = new TH2D(Tit2Den,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
	43	// set up ratio
	44	//title = "Ratio Qinv (MeV/c)";
	45	// this next bit is unfortunately needed so that we can have many histos of same "title"
	46	// it is neccessary if we typedef TH2D to TH1d (which we do)
	47	//mShareNumerator->SetDirectory(0);
	48	//mShareDenominator->SetDirectory(0);
	49	//mRatio->SetDirectory(0);
	50
	51	// to enable error bar calculation...
	52	fShareNumerator->Sumw2();
	53	fShareDenominator->Sumw2();
	54
	55	fQualityNumerator->Sumw2();
	56	fQualityDenominator->Sumw2();
	57	}
	58
0215f606	59	//____________________________
	60	AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(const AliFemtoShareQualityCorrFctn& aCorrFctn) :
	61	fShareNumerator(0),
	62	fShareDenominator(0),
	63	fQualityNumerator(0),
	64	fQualityDenominator(0)
	65	{
	66	if (aCorrFctn.fShareNumerator)
	67	fShareNumerator = new TH2D(*aCorrFctn.fShareNumerator);
	68	if (aCorrFctn.fShareDenominator)
	69	fShareDenominator = new TH2D(*aCorrFctn.fShareDenominator);
	70	if (aCorrFctn.fQualityNumerator)
	71	fQualityNumerator = new TH2D(*aCorrFctn.fQualityNumerator);
	72	if (aCorrFctn.fQualityDenominator)
	73	fQualityDenominator = new TH2D(*aCorrFctn.fQualityDenominator);
	74	}
67427ff7	75	//____________________________
	76	AliFemtoShareQualityCorrFctn::~AliFemtoShareQualityCorrFctn(){
	77	delete fShareNumerator;
	78	delete fShareDenominator;
	79	delete fQualityNumerator;
	80	delete fQualityDenominator;
	81	}
	82	//_________________________
0215f606	83	AliFemtoShareQualityCorrFctn& AliFemtoShareQualityCorrFctn::operator=(const AliFemtoShareQualityCorrFctn& aCorrFctn)
	84	{
	85	if (this == &aCorrFctn)
	86	return *this;
	87
	88	if (aCorrFctn.fShareNumerator)
	89	fShareNumerator = new TH2D(*aCorrFctn.fShareNumerator);
	90	else
	91	fShareNumerator = 0;
	92	if (aCorrFctn.fShareDenominator)
	93	fShareDenominator = new TH2D(*aCorrFctn.fShareDenominator);
	94	else
	95	fShareDenominator = 0;
	96	if (aCorrFctn.fQualityNumerator)
	97	fQualityNumerator = new TH2D(*aCorrFctn.fQualityNumerator);
	98	else
	99	fQualityNumerator = 0;
	100	if (aCorrFctn.fQualityDenominator)
	101	fQualityDenominator = new TH2D(*aCorrFctn.fQualityDenominator);
	102	else
	103	fQualityDenominator = 0;
	104
	105	return *this;
	106	}
	107	//_________________________
67427ff7	108	void AliFemtoShareQualityCorrFctn::Finish(){
	109	// here is where we should normalize, fit, etc...
	110	// we should NOT Draw() the histos (as I had done it below),
	111	// since we want to insulate ourselves from root at this level
	112	// of the code. Do it instead at root command line with browser.
	113	// mShareNumerator->Draw();
	114	//mShareDenominator->Draw();
	115	//mRatio->Draw();
	116
	117	}
	118
	119	//____________________________
	120	AliFemtoString AliFemtoShareQualityCorrFctn::Report(){
	121	string stemp = "Qinv Correlation Function Report:\n";
	122	char ctemp[100];
	123	sprintf(ctemp,"Number of entries in numerator:\t%E\n",fShareNumerator->GetEntries());
	124	stemp += ctemp;
	125	sprintf(ctemp,"Number of entries in denominator:\t%E\n",fShareDenominator->GetEntries());
	126	stemp += ctemp;
	127	// stemp += mCoulombWeight->Report();
	128	AliFemtoString returnThis = stemp;
	129	return returnThis;
	130	}
	131	//____________________________
	132	void AliFemtoShareQualityCorrFctn::AddRealPair(const AliFemtoPair* pair){
	133	double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
	134	Int_t nh = 0;
	135	Int_t an = 0;
	136	Int_t ns = 0;
	137
	138	for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
	139	// If both have clusters in the same row
	140	if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) &&
	141	pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
	142	// Do they share it ?
	143	if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap) \|\|
	144	pair->track2()->Track()->TPCsharing().TestBitNumber(imap))
	145	{
	146	if (Qinv < 0.01) {
	147	cout << "Shared cluster in row " << imap << endl;
	148	}
	149	an++;
	150	nh+=2;
	151	ns+=2;
	152	}
	153
	154	// Different hits on the same padrow
	155	else {
	156	an--;
	157	nh+=2;
	158	}
	159	}
	160	else if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) \|\|
	161	pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
	162	// One track has a hit, the other does not
	163	an++;
	164	nh++;
	165	}
	166	}
	167	if (Qinv < 0.01) {
	168	cout << "Qinv of the pair is " << Qinv << endl;
	169	cout << "Clusters: " << endl;
	170	for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
	171	cout << imap ;
172	if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
173	else cout << " 0 " ;
174	if (pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
175	else cout << " 0 " ;
176	cout << " ";
177	if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
178	else cout << " X ";
179	if (pair->track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
180	else cout << " X ";
181	cout << endl;
182	}
183	}
184
185	Float_t hsmval = 0.0;
186	Float_t hsfval = 0.0;
187
188	if (nh >0) {
189	hsmval = an*1.0/nh;
190	hsfval = ns*1.0/nh;
191	}
192
193	if (Qinv < 0.01) {
194	cout << "Quality Sharity " << hsmval << " " << hsfval << " " << pair->track1()->Track() << " " << pair->track2()->Track() << endl;
195	}
196
197	fShareNumerator->Fill(Qinv, hsfval);
198	fQualityNumerator->Fill(Qinv, hsmval);
199	// cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << Qinv <<
200	//" " << pair->track1().FourMomentum() << " " << pair->track2().FourMomentum() << endl;
201	}
202	//____________________________
203	void AliFemtoShareQualityCorrFctn::AddMixedPair(const AliFemtoPair* pair){
204	double weight = 1.0;
205	double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
206	Int_t nh = 0;
207	Int_t an = 0;
208	Int_t ns = 0;
209
210	for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
211	// If both have clusters in the same row
212	if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) &&
213	pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
214	// Do they share it ?
215	if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap) \|\|
216	pair->track2()->Track()->TPCsharing().TestBitNumber(imap))
217	{
218	// cout << "A shared cluster !!!" << endl;
219	// cout << "imap idx1 idx2 "
220	// << imap << " "
221	// << tP1idx[imap] << " " << tP2idx[imap] << endl;
222	an++;
223	nh+=2;
224	ns+=2;
225	}
226
227	// Different hits on the same padrow
228	else {
229	an--;
230	nh+=2;
231	}
232	}
233	else if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) \|\|
234	pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
235	// One track has a hit, the other does not
236	an++;
237	nh++;
238	}
239	}
240
241	Float_t hsmval = 0.0;
242	Float_t hsfval = 0.0;
243
244	if (nh >0) {
245	hsmval = an*1.0/nh;
246	hsfval = ns*1.0/nh;
247	}
248
249	fShareDenominator->Fill(Qinv,hsfval,weight);
250	fQualityDenominator->Fill(Qinv,hsmval,weight);
251	}
252
253
254	void AliFemtoShareQualityCorrFctn::WriteHistos()
255	{
256	fShareNumerator->Write();
257	fShareDenominator->Write();
258	fQualityNumerator->Write();
259	fQualityDenominator->Write();
260
261	}