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

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// AliFemtoShareQualityQAPairCut - a pair cut which checks for some pair     //
// qualities that attempt to identify slit/doubly reconstructed tracks     //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
/***************************************************************************
 *
 * $Id: AliFemtoShareQualityQAPairCut.cxx 24360 2008-03-10 09:48:27Z akisiel $
 *
 * Author: Adam Kisiel, Ohio State, kisiel@mps.ohio-state.edu
 ***************************************************************************
 *
 * Description: part of STAR HBT Framework: AliFemtoMaker package
 *   a cut to remove "shared" and "split" pairs
 *
 ***************************************************************************
 *
 *
 **************************************************************************/

#include "AliFemtoShareQualityQAPairCut.h"
#include <string>
#include <cstdio>

#ifdef __ROOT__
ClassImp(AliFemtoShareQualityQAPairCut)
#endif

//__________________
AliFemtoShareQualityQAPairCut::AliFemtoShareQualityQAPairCut():
  fNPairsPassed(0),
  fNPairsFailed(0),
  fShareQualityMax(1.0),
  fShareQualitymin(-0.5),
  fShareFractionMax(1.0),
  fShareFractionmin(0.0),
  fRemoveSameLabel(0),
  fShareQualityQASwitch(0),
  fShareFractionQASwitch(0)
{ 
  fShareQualityQASwitch  = false;
  fShareQualityQAExclusionZone[0] = -0.5;
  fShareQualityQAExclusionZone[1] = 1.0;
  fShareFractionQASwitch = false; 
  fShareFractionQAExclusionZone[0] = 0.0;
  fShareFractionQAExclusionZone[1] = 1.0;
}

//__________________
AliFemtoShareQualityQAPairCut::~AliFemtoShareQualityQAPairCut(){
  /* no-op */
}
//__________________
AliFemtoShareQualityQAPairCut& AliFemtoShareQualityQAPairCut::operator=(const AliFemtoShareQualityQAPairCut& cut)
{
  if (this != &cut) {
    AliFemtoPairCut::operator=(cut);
    fNPairsPassed = 0;
    fNPairsFailed = 0;
    fShareQualityMax = 1.0;
    fShareQualitymin = -0.5;
    fShareFractionMax = 1.0;
    fShareFractionmin = 0.0;
    fRemoveSameLabel = 0;
    fShareQualityQASwitch = 0;
    fShareFractionQASwitch = 0;
    fShareQualityQASwitch  = cut.fShareQualityQASwitch;
    fShareQualityQAExclusionZone[0]  = cut.fShareQualityQAExclusionZone[0];
    fShareQualityQAExclusionZone[1]  = cut.fShareQualityQAExclusionZone[1];
    fShareFractionQASwitch = cut.fShareFractionQASwitch; 
    fShareFractionQAExclusionZone[0]  = cut.fShareFractionQAExclusionZone[0];
    fShareFractionQAExclusionZone[1]  = cut.fShareFractionQAExclusionZone[1];
  }

  return *this;
  
}
//__________________
bool AliFemtoShareQualityQAPairCut::Pass(const AliFemtoPair* pair){
  // Check for pairs that are possibly shared/double reconstruction
  bool pass;
  
  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;
  }
  //  if (hsmval > -0.4) {
//   cout << "Pair quality: " << hsmval << " " << an << " " << nh << " " 
//        << (pair->Track1()->Track()) << " " 
//        << (pair->Track2()->Track()) << endl;
//   cout << "Bits: " << pair->Track1()->Track()->TPCclusters().GetNbits() << endl;
    //  }
//   if (hsfval > 0.0) {
//     cout << "Pair sharity: " << hsfval << " " << ns << " " << nh << "    " << hsmval << " " << an << " " << nh << endl;
//   }

  // Determine if pair pass/fail cuts:
  if (fShareQualityQASwitch) {
    pass = (((hsmval >= fShareQualitymin)  && (hsmval < fShareQualityQAExclusionZone[0]))  || 
	    ((hsmval >= fShareQualityQAExclusionZone[1])  && (hsmval < fShareQualityMax))) &&
      (hsfval >= fShareFractionmin) && (hsfval < fShareFractionMax);
  }
  else if (fShareFractionQASwitch) {
    pass = (((hsfval >= fShareFractionmin)  && (hsfval < fShareFractionQAExclusionZone[0]))  || 
	    ((hsfval >= fShareFractionQAExclusionZone[1])  && (hsfval < fShareFractionMax))) &&
      (hsmval >= fShareQualitymin) && (hsmval < fShareQualityMax);
  }
  else {
    pass = (hsmval >= fShareQualitymin)  && (hsmval < fShareQualityMax) && 
      (hsfval >= fShareFractionmin) && (hsfval < fShareFractionMax);
  }

  if (fRemoveSameLabel) {
    if (abs(pair->Track1()->Track()->Label()) == abs(pair->Track2()->Track()->Label())) {
      cout << "Found a pair with same label " << pair->Track1()->Track()->Label() << endl;
      cout << "Quality Sharity Passed " << hsmval << " " << hsfval << " " << pair->QInv() << " " << pass << endl;
      pass = kFALSE;
    }
  }

  pass ? fNPairsPassed++ : fNPairsFailed++;
  return pass;
}
//__________________
AliFemtoString AliFemtoShareQualityQAPairCut::Report(){
  // Prepare the report from the execution
  string stemp = "AliFemtoShareQuality Pair Cut - remove shared and split pairs\n";  char ctemp[100];
  snprintf(ctemp , 100, "Number of pairs which passed:\t%ld  Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed);
  stemp += ctemp;
  AliFemtoString returnThis = stemp;
  return returnThis;}
//__________________

void AliFemtoShareQualityQAPairCut::SetShareQualityMax(Double_t aShareQualityMax) {
  fShareQualityMax = aShareQualityMax;
}

void AliFemtoShareQualityQAPairCut::SetShareQualitymin(Double_t aShareQualitymin) {
  fShareQualitymin = aShareQualitymin;
}

void AliFemtoShareQualityQAPairCut::SetShareQualityQASwitch(bool aSwitch) {
  fShareQualityQASwitch = aSwitch;
}

void AliFemtoShareQualityQAPairCut::SetShareQualityQAExclusionZone(Double_t lo, Double_t hi) {
  fShareQualityQAExclusionZone[0] = lo;
  fShareQualityQAExclusionZone[1] = hi;
}

Double_t AliFemtoShareQualityQAPairCut::GetAliFemtoShareQualityMax() const {
  return fShareQualityMax;
}

void AliFemtoShareQualityQAPairCut::SetShareFractionMax(Double_t aShareFractionMax) {
  fShareFractionMax = aShareFractionMax;
}

void AliFemtoShareQualityQAPairCut::SetShareFractionmin(Double_t aShareFractionmin) {
  fShareFractionmin = aShareFractionmin;
}

void AliFemtoShareQualityQAPairCut::SetShareFractionQASwitch(bool aSwitch) {
  fShareFractionQASwitch = aSwitch;
}

void AliFemtoShareQualityQAPairCut::SetShareFractionQAExclusionZone(Double_t lo, Double_t hi) {
  fShareFractionQAExclusionZone[0] = lo;
  fShareFractionQAExclusionZone[1] = hi;
}

Double_t AliFemtoShareQualityQAPairCut::GetAliFemtoShareFractionMax() const {
  return fShareFractionMax;
}

TList *AliFemtoShareQualityQAPairCut::ListSettings()
{
  // return a list of settings in a writable form
  TList *tListSetttings = new TList();
  char buf[200];
  snprintf(buf, 200, "AliFemtoShareQualityQAPairCut.sharequalitymax=%f", fShareQualityMax);
  snprintf(buf, 200, "AliFemtoShareQualityQAPairCut.sharefractionmax=%f", fShareFractionMax);
  tListSetttings->AddLast(new TObjString(buf));

  return tListSetttings;
}

void     AliFemtoShareQualityQAPairCut::SetRemoveSameLabel(Bool_t aRemove)
{
  fRemoveSameLabel = aRemove;
}
Commit	Line	Data
f9210de7	1	/////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoShareQualityQAPairCut - a pair cut which checks for some pair //
	4	// qualities that attempt to identify slit/doubly reconstructed tracks //
	5	// //
	6	/////////////////////////////////////////////////////////////////////////////
	7	/***************************************************************************
	8	*
	9	* $Id: AliFemtoShareQualityQAPairCut.cxx 24360 2008-03-10 09:48:27Z akisiel $
	10	*
	11	* Author: Adam Kisiel, Ohio State, kisiel@mps.ohio-state.edu
	12	***************************************************************************
	13	*
	14	* Description: part of STAR HBT Framework: AliFemtoMaker package
	15	* a cut to remove "shared" and "split" pairs
	16	*
	17	***************************************************************************
	18	*
	19	*
	20	**************************************************************************/
	21
	22	#include "AliFemtoShareQualityQAPairCut.h"
	23	#include <string>
	24	#include <cstdio>
	25
	26	#ifdef __ROOT__
	27	ClassImp(AliFemtoShareQualityQAPairCut)
	28	#endif
	29
	30	//__________________
	31	AliFemtoShareQualityQAPairCut::AliFemtoShareQualityQAPairCut():
	32	fNPairsPassed(0),
	33	fNPairsFailed(0),
	34	fShareQualityMax(1.0),
	35	fShareQualitymin(-0.5),
	36	fShareFractionMax(1.0),
	37	fShareFractionmin(0.0),
1e79fae4	38	fRemoveSameLabel(0),
	39	fShareQualityQASwitch(0),
	40	fShareFractionQASwitch(0)
f9210de7	41	{
	42	fShareQualityQASwitch = false;
	43	fShareQualityQAExclusionZone[0] = -0.5;
	44	fShareQualityQAExclusionZone[1] = 1.0;
	45	fShareFractionQASwitch = false;
	46	fShareFractionQAExclusionZone[0] = 0.0;
	47	fShareFractionQAExclusionZone[1] = 1.0;
	48	}
	49
	50	//__________________
	51	AliFemtoShareQualityQAPairCut::~AliFemtoShareQualityQAPairCut(){
	52	/* no-op */
	53	}
	54	//__________________
83e74701	55	AliFemtoShareQualityQAPairCut& AliFemtoShareQualityQAPairCut::operator=(const AliFemtoShareQualityQAPairCut& cut)
	56	{
	57	if (this != &cut) {
	58	AliFemtoPairCut::operator=(cut);
	59	fNPairsPassed = 0;
	60	fNPairsFailed = 0;
	61	fShareQualityMax = 1.0;
	62	fShareQualitymin = -0.5;
	63	fShareFractionMax = 1.0;
	64	fShareFractionmin = 0.0;
	65	fRemoveSameLabel = 0;
	66	fShareQualityQASwitch = 0;
	67	fShareFractionQASwitch = 0;
	68	fShareQualityQASwitch = cut.fShareQualityQASwitch;
	69	fShareQualityQAExclusionZone[0] = cut.fShareQualityQAExclusionZone[0];
	70	fShareQualityQAExclusionZone[1] = cut.fShareQualityQAExclusionZone[1];
	71	fShareFractionQASwitch = cut.fShareFractionQASwitch;
	72	fShareFractionQAExclusionZone[0] = cut.fShareFractionQAExclusionZone[0];
	73	fShareFractionQAExclusionZone[1] = cut.fShareFractionQAExclusionZone[1];
	74	}
	75
	76	return *this;
	77
	78	}
	79	//__________________
f9210de7	80	bool AliFemtoShareQualityQAPairCut::Pass(const AliFemtoPair* pair){
	81	// Check for pairs that are possibly shared/double reconstruction
	82	bool pass;
	83
	84	Int_t nh = 0;
	85	Int_t an = 0;
	86	Int_t ns = 0;
	87
	88	for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
	89	// If both have clusters in the same row
	90	if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) &&
	91	pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
	92	// Do they share it ?
	93	if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) &&
	94	pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
	95	{
	96	// cout << "A shared cluster !!!" << endl;
	97	// cout << "imap idx1 idx2 "
	98	// << imap << " "
	99	// << tP1idx[imap] << " " << tP2idx[imap] << endl;
	100	an++;
	101	nh+=2;
	102	ns+=2;
	103	}
	104
	105	// Different hits on the same padrow
	106	else {
	107	an--;
	108	nh+=2;
	109	}
	110	}
	111	else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) \|\|
	112	pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
	113	// One track has a hit, the other does not
	114	an++;
	115	nh++;
	116	}
	117	}
	118
	119	Float_t hsmval = 0.0;
	120	Float_t hsfval = 0.0;
	121
	122	if (nh >0) {
	123	hsmval = an*1.0/nh;
	124	hsfval = ns*1.0/nh;
	125	}
	126	// if (hsmval > -0.4) {
	127	// cout << "Pair quality: " << hsmval << " " << an << " " << nh << " "
	128	// << (pair->Track1()->Track()) << " "
	129	// << (pair->Track2()->Track()) << endl;
	130	// cout << "Bits: " << pair->Track1()->Track()->TPCclusters().GetNbits() << endl;
	131	// }
	132	// if (hsfval > 0.0) {
	133	// cout << "Pair sharity: " << hsfval << " " << ns << " " << nh << " " << hsmval << " " << an << " " << nh << endl;
	134	// }
	135
	136	// Determine if pair pass/fail cuts:
	137	if (fShareQualityQASwitch) {
62029d32	138	pass = (((hsmval >= fShareQualitymin) && (hsmval < fShareQualityQAExclusionZone[0])) \|\|
	139	((hsmval >= fShareQualityQAExclusionZone[1]) && (hsmval < fShareQualityMax))) &&
	140	(hsfval >= fShareFractionmin) && (hsfval < fShareFractionMax);
f9210de7	141	}
f9210de7	142	else if (fShareFractionQASwitch) {
62029d32	143	pass = (((hsfval >= fShareFractionmin) && (hsfval < fShareFractionQAExclusionZone[0])) \|\|
	144	((hsfval >= fShareFractionQAExclusionZone[1]) && (hsfval < fShareFractionMax))) &&
	145	(hsmval >= fShareQualitymin) && (hsmval < fShareQualityMax);
f9210de7	146	}
f9210de7	147	else {
62029d32	148	pass = (hsmval >= fShareQualitymin) && (hsmval < fShareQualityMax) &&
62029d32	149	(hsfval >= fShareFractionmin) && (hsfval < fShareFractionMax);
f9210de7	150	}
	151
	152	if (fRemoveSameLabel) {
	153	if (abs(pair->Track1()->Track()->Label()) == abs(pair->Track2()->Track()->Label())) {
	154	cout << "Found a pair with same label " << pair->Track1()->Track()->Label() << endl;
	155	cout << "Quality Sharity Passed " << hsmval << " " << hsfval << " " << pair->QInv() << " " << pass << endl;
	156	pass = kFALSE;
	157	}
	158	}
	159
	160	pass ? fNPairsPassed++ : fNPairsFailed++;
	161	return pass;
	162	}
	163	//__________________
	164	AliFemtoString AliFemtoShareQualityQAPairCut::Report(){
	165	// Prepare the report from the execution
	166	string stemp = "AliFemtoShareQuality Pair Cut - remove shared and split pairs\n"; char ctemp[100];
adecdc37	167	snprintf(ctemp , 100, "Number of pairs which passed:\t%ld Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed);
f9210de7	168	stemp += ctemp;
	169	AliFemtoString returnThis = stemp;
	170	return returnThis;}
	171	//__________________
	172
	173	void AliFemtoShareQualityQAPairCut::SetShareQualityMax(Double_t aShareQualityMax) {
	174	fShareQualityMax = aShareQualityMax;
	175	}
	176
	177	void AliFemtoShareQualityQAPairCut::SetShareQualitymin(Double_t aShareQualitymin) {
	178	fShareQualitymin = aShareQualitymin;
	179	}
	180
	181	void AliFemtoShareQualityQAPairCut::SetShareQualityQASwitch(bool aSwitch) {
	182	fShareQualityQASwitch = aSwitch;
	183	}
	184
	185	void AliFemtoShareQualityQAPairCut::SetShareQualityQAExclusionZone(Double_t lo, Double_t hi) {
	186	fShareQualityQAExclusionZone[0] = lo;
	187	fShareQualityQAExclusionZone[1] = hi;
	188	}
	189
	190	Double_t AliFemtoShareQualityQAPairCut::GetAliFemtoShareQualityMax() const {
	191	return fShareQualityMax;
	192	}
	193
	194	void AliFemtoShareQualityQAPairCut::SetShareFractionMax(Double_t aShareFractionMax) {
	195	fShareFractionMax = aShareFractionMax;
	196	}
	197
	198	void AliFemtoShareQualityQAPairCut::SetShareFractionmin(Double_t aShareFractionmin) {
	199	fShareFractionmin = aShareFractionmin;
	200	}
	201
	202	void AliFemtoShareQualityQAPairCut::SetShareFractionQASwitch(bool aSwitch) {
	203	fShareFractionQASwitch = aSwitch;
	204	}
	205
	206	void AliFemtoShareQualityQAPairCut::SetShareFractionQAExclusionZone(Double_t lo, Double_t hi) {
	207	fShareFractionQAExclusionZone[0] = lo;
	208	fShareFractionQAExclusionZone[1] = hi;
	209	}
	210
	211	Double_t AliFemtoShareQualityQAPairCut::GetAliFemtoShareFractionMax() const {
	212	return fShareFractionMax;
	213	}
	214
	215	TList *AliFemtoShareQualityQAPairCut::ListSettings()
	216	{
	217	// return a list of settings in a writable form
	218	TList *tListSetttings = new TList();
	219	char buf[200];
	220	snprintf(buf, 200, "AliFemtoShareQualityQAPairCut.sharequalitymax=%f", fShareQualityMax);
	221	snprintf(buf, 200, "AliFemtoShareQualityQAPairCut.sharefractionmax=%f", fShareFractionMax);
	222	tListSetttings->AddLast(new TObjString(buf));
	223
	224	return tListSetttings;
	225	}
	226
	227	void AliFemtoShareQualityQAPairCut::SetRemoveSameLabel(Bool_t aRemove)
	228	{
	229	fRemoveSameLabel = aRemove;
	230	}