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