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

/***************************************************************************
 *
 * $Id$
 *
 * Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
 ***************************************************************************
 *
 * Description: part of STAR HBT Framework: AliFemtoMaker package
 *   a simple Q-invariant correlation function           
 *
 ***************************************************************************
 *
 * $Log$
 * Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
 * First version on CVS
 *
 * Revision 1.4  2000/01/25 17:34:45  laue
 * I. In order to run the stand alone version of the AliFemtoMaker the following
 * changes have been done:
 * a) all ClassDefs and ClassImps have been put into #ifdef __ROOT__ statements
 * b) unnecessary includes of StMaker.h have been removed
 * c) the subdirectory AliFemtoMaker/doc/Make has been created including everything
 * needed for the stand alone version
 *
 * II. To reduce the amount of compiler warning
 * a) some variables have been type casted
 * b) some destructors have been declared as virtual
 *
 * Revision 1.3  1999/07/29 02:47:09  lisa
 * 1) add OpeningAngle correlation function 2) add AliFemtoMcEventReader 3) make histos in CorrFctns do errors correctly
 *
 * Revision 1.2  1999/07/06 22:33:20  lisa
 * Adjusted all to work in pro and new - dev itself is broken
 *
 * Revision 1.1.1.1  1999/06/29 16:02:57  lisa
 * Installation of AliFemtoMaker
 *
 **************************************************************************/

#include "CorrFctn/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){
  // 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(){
  delete fShareNumerator;
  delete fShareDenominator;
  delete fQualityNumerator;
  delete fQualityDenominator;
}
//_________________________
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
67427ff7	1	/***************************************************************************
	2	*
	3	* $Id$
	4	*
	5	* Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
	6	***************************************************************************
	7	*
	8	* Description: part of STAR HBT Framework: AliFemtoMaker package
	9	* a simple Q-invariant correlation function
	10	*
	11	***************************************************************************
	12	*
	13	* $Log$
	14	* Revision 1.1.1.1 2007/03/07 10:14:49 mchojnacki
	15	* First version on CVS
	16	*
	17	* Revision 1.4 2000/01/25 17:34:45 laue
	18	* I. In order to run the stand alone version of the AliFemtoMaker the following
	19	* changes have been done:
	20	* a) all ClassDefs and ClassImps have been put into #ifdef __ROOT__ statements
	21	* b) unnecessary includes of StMaker.h have been removed
	22	* c) the subdirectory AliFemtoMaker/doc/Make has been created including everything
	23	* needed for the stand alone version
	24	*
	25	* II. To reduce the amount of compiler warning
	26	* a) some variables have been type casted
	27	* b) some destructors have been declared as virtual
	28	*
	29	* Revision 1.3 1999/07/29 02:47:09 lisa
	30	* 1) add OpeningAngle correlation function 2) add AliFemtoMcEventReader 3) make histos in CorrFctns do errors correctly
	31	*
	32	* Revision 1.2 1999/07/06 22:33:20 lisa
	33	* Adjusted all to work in pro and new - dev itself is broken
	34	*
	35	* Revision 1.1.1.1 1999/06/29 16:02:57 lisa
	36	* Installation of AliFemtoMaker
	37	*
	38	**************************************************************************/
	39
	40	#include "CorrFctn/AliFemtoShareQualityCorrFctn.h"
	41	//#include "Infrastructure/AliFemtoHisto.hh"
	42	#include <cstdio>
	43
	44	#ifdef __ROOT__
	45	ClassImp(AliFemtoShareQualityCorrFctn)
	46	#endif
	47
	48	//____________________________
	49	AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi){
	50	// set up numerator
	51	// title = "Num Qinv (MeV/c)";
	52	char TitNum[100] = "NumShare";
	53	strcat(TitNum,title);
	54	fShareNumerator = new TH2D(TitNum,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
	55	// set up denominator
	56	//title = "Den Qinv (MeV/c)";
	57	char TitDen[100] = "DenShare";
	58	strcat(TitDen,title);
	59	fShareDenominator = new TH2D(TitDen,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
	60
	61	char Tit2Num[100] = "NumQuality";
	62	strcat(Tit2Num,title);
	63	fQualityNumerator = new TH2D(Tit2Num,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
	64	// set up denominator
65	//title = "Den Qinv (MeV/c)";
66	char Tit2Den[100] = "DenQuality";
67	strcat(Tit2Den,title);
68	fQualityDenominator = new TH2D(Tit2Den,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
69	// set up ratio
70	//title = "Ratio Qinv (MeV/c)";
71	// this next bit is unfortunately needed so that we can have many histos of same "title"
72	// it is neccessary if we typedef TH2D to TH1d (which we do)
73	//mShareNumerator->SetDirectory(0);
74	//mShareDenominator->SetDirectory(0);
75	//mRatio->SetDirectory(0);
76
77	// to enable error bar calculation...
78	fShareNumerator->Sumw2();
79	fShareDenominator->Sumw2();
80
81	fQualityNumerator->Sumw2();
82	fQualityDenominator->Sumw2();
83	}
84
85	//____________________________
86	AliFemtoShareQualityCorrFctn::~AliFemtoShareQualityCorrFctn(){
87	delete fShareNumerator;
88	delete fShareDenominator;
89	delete fQualityNumerator;
90	delete fQualityDenominator;
91	}
92	//_________________________
93	void AliFemtoShareQualityCorrFctn::Finish(){
94	// here is where we should normalize, fit, etc...
95	// we should NOT Draw() the histos (as I had done it below),
96	// since we want to insulate ourselves from root at this level
97	// of the code. Do it instead at root command line with browser.
98	// mShareNumerator->Draw();
99	//mShareDenominator->Draw();
100	//mRatio->Draw();
101
102	}
103
104	//____________________________
105	AliFemtoString AliFemtoShareQualityCorrFctn::Report(){
106	string stemp = "Qinv Correlation Function Report:\n";
107	char ctemp[100];
108	sprintf(ctemp,"Number of entries in numerator:\t%E\n",fShareNumerator->GetEntries());
109	stemp += ctemp;
110	sprintf(ctemp,"Number of entries in denominator:\t%E\n",fShareDenominator->GetEntries());
111	stemp += ctemp;
112	// stemp += mCoulombWeight->Report();
113	AliFemtoString returnThis = stemp;
114	return returnThis;
115	}
116	//____________________________
117	void AliFemtoShareQualityCorrFctn::AddRealPair(const AliFemtoPair* pair){
118	double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
119	Int_t nh = 0;
120	Int_t an = 0;
121	Int_t ns = 0;
122
123	for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
124	// If both have clusters in the same row
125	if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) &&
126	pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
127	// Do they share it ?
128	if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap) \|\|
129	pair->track2()->Track()->TPCsharing().TestBitNumber(imap))
130	{
131	if (Qinv < 0.01) {
132	cout << "Shared cluster in row " << imap << endl;
133	}
134	an++;
135	nh+=2;
136	ns+=2;
137	}
138
139	// Different hits on the same padrow
140	else {
141	an--;
142	nh+=2;
143	}
144	}
145	else if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) \|\|
146	pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
147	// One track has a hit, the other does not
148	an++;
149	nh++;
150	}
151	}
152	if (Qinv < 0.01) {
153	cout << "Qinv of the pair is " << Qinv << endl;
154	cout << "Clusters: " << endl;
155	for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
156	cout << imap ;
157	if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
158	else cout << " 0 " ;
159	if (pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
160	else cout << " 0 " ;
161	cout << " ";
162	if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
163	else cout << " X ";
164	if (pair->track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
165	else cout << " X ";
166	cout << endl;
167	}
168	}
169
170	Float_t hsmval = 0.0;
171	Float_t hsfval = 0.0;
172
173	if (nh >0) {
174	hsmval = an*1.0/nh;
175	hsfval = ns*1.0/nh;
176	}
177
178	if (Qinv < 0.01) {
179	cout << "Quality Sharity " << hsmval << " " << hsfval << " " << pair->track1()->Track() << " " << pair->track2()->Track() << endl;
180	}
181
182	fShareNumerator->Fill(Qinv, hsfval);
183	fQualityNumerator->Fill(Qinv, hsmval);
184	// cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << Qinv <<
185	//" " << pair->track1().FourMomentum() << " " << pair->track2().FourMomentum() << endl;
186	}
187	//____________________________
188	void AliFemtoShareQualityCorrFctn::AddMixedPair(const AliFemtoPair* pair){
189	double weight = 1.0;
190	double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
191	Int_t nh = 0;
192	Int_t an = 0;
193	Int_t ns = 0;
194
195	for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
196	// If both have clusters in the same row
197	if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) &&
198	pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
199	// Do they share it ?
200	if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap) \|\|
201	pair->track2()->Track()->TPCsharing().TestBitNumber(imap))
202	{
203	// cout << "A shared cluster !!!" << endl;
204	// cout << "imap idx1 idx2 "
205	// << imap << " "
206	// << tP1idx[imap] << " " << tP2idx[imap] << endl;
207	an++;
208	nh+=2;
209	ns+=2;
210	}
211
212	// Different hits on the same padrow
213	else {
214	an--;
215	nh+=2;
216	}
217	}
218	else if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) \|\|
219	pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
220	// One track has a hit, the other does not
221	an++;
222	nh++;
223	}
224	}
225
226	Float_t hsmval = 0.0;
227	Float_t hsfval = 0.0;
228
229	if (nh >0) {
230	hsmval = an*1.0/nh;
231	hsfval = ns*1.0/nh;
232	}
233
234	fShareDenominator->Fill(Qinv,hsfval,weight);
235	fQualityDenominator->Fill(Qinv,hsmval,weight);
236	}
237
238
239	void AliFemtoShareQualityCorrFctn::WriteHistos()
240	{
241	fShareNumerator->Write();
242	fShareDenominator->Write();
243	fQualityNumerator->Write();
244	fQualityDenominator->Write();
245
246	}