[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemto / AliFemtoV0TrackPairCut.cxx

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// AliFemtoShareQualityPairCut - a pair cut which checks for some pair     //
// qualities that attempt to identify slit/doubly reconstructed tracks     //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
/***************************************************************************
 *
 * $Id: AliFemtoShareQualityPairCut.cxx 50722 2011-07-21 15:18:38Z 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 "AliFemtoV0TrackPairCut.h"
#include <string>
#include <cstdio>

#ifdef __ROOT__
ClassImp(AliFemtoV0TrackPairCut)
#endif

//__________________
AliFemtoV0TrackPairCut::AliFemtoV0TrackPairCut():
  fNPairsPassed(0),
  fNPairsFailed(0),
  fV0Max(1.0),
  fShareQualityMax(1.0),
  fShareFractionMax(1.0),
  fRemoveSameLabel(0),
  fTrackTPCOnly(0),
  fDataType(kAOD),
  fDTPCMin(0),
  fDTPCExitMin(0),
  fKstarCut(0),
  fFirstParticleType(kLambda), 
  fSecondParticleType(kProton),
  fMinAvgSepTrackPos(0),
  fMinAvgSepTrackNeg(0)
{
  // Default constructor
  // Nothing to do
}
//__________________
AliFemtoV0TrackPairCut::~AliFemtoV0TrackPairCut(){
  /* no-op */
}
AliFemtoV0TrackPairCut& AliFemtoV0TrackPairCut::operator=(const AliFemtoV0TrackPairCut& cut)
{
  if (this != &cut) {
   
    AliFemtoPairCut::operator=(cut);
    fNPairsPassed = 0;
    fNPairsFailed =0;
    fV0Max = 1.0;
    fShareQualityMax = 1.0;
    fShareFractionMax = 1.0;
    fRemoveSameLabel = 0;
    fTrackTPCOnly = 0;
    fDataType = kAOD;
    fDTPCMin = 0;
    fDTPCExitMin = 0;
    fMinAvgSepTrackPos = 0;
    fMinAvgSepTrackNeg = 0;
  }

  return *this;
}
//__________________
bool AliFemtoV0TrackPairCut::Pass(const AliFemtoPair* pair){
  // Check for pairs that are possibly shared/double reconstruction

  bool temp = true;
  //Track1 - V0
  //Track2 - track

  if(!(pair->Track1()->V0() && pair->Track2()->Track()))
    {
      return false;
    }
  if(fTrackTPCOnly)
    {
      if( (-(pair->Track1()->V0()->IdNeg()+1)) ==pair->Track2()->TrackId() || (-(pair->Track1()->V0()->IdPos()+1)) ==pair->Track2()->TrackId())
	{
	  return false;
	}
    }
  else
    {
      if(pair->Track1()->V0()->IdNeg()==pair->Track2()->TrackId() || pair->Track1()->V0()->IdPos()==pair->Track2()->TrackId())
	{
	  return false;
	}
    }
  
  bool tempTPCEntrancePos = true;
  bool tempTPCEntranceNeg = true;
  bool tempTPCExitPos = true;
  bool tempTPCExitNeg = true;
  if(fDataType==kESD || fDataType==kAOD)
    {
      double distx = pair->Track1()->V0()->NominalTpcEntrancePointPos().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
      double disty = pair->Track1()->V0()->NominalTpcEntrancePointPos().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
      double distz = pair->Track1()->V0()->NominalTpcEntrancePointPos().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
      double distPos = sqrt(distx*distx + disty*disty + distz*distz);

      distx = pair->Track1()->V0()->NominalTpcEntrancePointNeg().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
      disty = pair->Track1()->V0()->NominalTpcEntrancePointNeg().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
      distz = pair->Track1()->V0()->NominalTpcEntrancePointNeg().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
      double distNeg = sqrt(distx*distx + disty*disty + distz*distz);

      double distExitX = pair->Track1()->V0()->NominalTpcExitPointPos().x() - pair->Track2()->Track()->NominalTpcExitPoint().x();
      double distExitY = pair->Track1()->V0()->NominalTpcExitPointPos().y() - pair->Track2()->Track()->NominalTpcExitPoint().y();
      double distExitZ = pair->Track1()->V0()->NominalTpcExitPointPos().z() - pair->Track2()->Track()->NominalTpcExitPoint().z();
      double distExitPos = sqrt(distExitX*distExitX + distExitY*distExitY + distExitZ*distExitZ);

      distExitX = pair->Track1()->V0()->NominalTpcExitPointNeg().x() - pair->Track2()->Track()->NominalTpcExitPoint().x();
      distExitY = pair->Track1()->V0()->NominalTpcExitPointNeg().y() - pair->Track2()->Track()->NominalTpcExitPoint().y();
      distExitZ = pair->Track1()->V0()->NominalTpcExitPointNeg().z() - pair->Track2()->Track()->NominalTpcExitPoint().z();
      double distExitNeg = sqrt(distExitX*distExitX + distExitY*distExitY + distExitZ*distExitZ);

      tempTPCEntrancePos = distPos > fDTPCMin;
      tempTPCEntranceNeg = distNeg > fDTPCMin;
      tempTPCExitPos = distExitPos > fDTPCExitMin;
      tempTPCExitNeg = distExitNeg > fDTPCExitMin;
    }
 
  if(!tempTPCEntrancePos || !tempTPCEntranceNeg || !tempTPCExitPos || !tempTPCExitNeg) return false;

  //reject merged trakcs in TPC
  

  //temp = dist > fDTPCMin;
  //koniec kopii
  //if(!temp)
  //return false;

  //kopia z AliFemtoShareQualityPairCut.cxx
  Int_t nh = 0;
  Int_t an = 0;
  Int_t ns = 0;

 
  if ((fShareFractionMax < 1.0) && ( fShareQualityMax < 1.0)) {
    for (unsigned int imap=0; imap<pair->Track1()->V0()->TPCclustersPos().GetNbits(); imap++) {
      // If both have clusters in the same row
      if (pair->Track1()->V0()->TPCclustersPos().TestBitNumber(imap) && 
	  pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
	// Do they share it ?
	if (pair->Track1()->V0()->TPCsharingPos().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()->V0()->TPCclustersPos().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;
    //   }
    
    temp = (hsmval < fShareQualityMax) && (hsfval < fShareFractionMax);
    if(!temp) return false;

    nh = 0;
    an = 0;
    ns = 0;

    for (unsigned int imap=0; imap<pair->Track1()->V0()->TPCclustersNeg().GetNbits(); imap++) {
      // If both have clusters in the same row
      if (pair->Track1()->V0()->TPCclustersNeg().TestBitNumber(imap) && 
	  pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
	// Do they share it ?
	if (pair->Track1()->V0()->TPCsharingNeg().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()->V0()->TPCclustersNeg().TestBitNumber(imap) ||
	       pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
	// One track has a hit, the other does not
	an++;
	nh++;
      }
    }
    
    hsmval = 0.0;
    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;
    //   }
    
    temp = (hsmval < fShareQualityMax) && (hsfval < fShareFractionMax);


  }
  else
    temp = true;
  //koniec kopii


  //avg sep pair cut
  double avgSep=0;
  AliFemtoThreeVector first, second, tmp;
  for(int i=0; i<8 ;i++)
    {
      tmp = pair->Track1()->V0()->NominalTpcPointPos(i);
      //cout<<"X pos: "<<tmp.x()<<endl;
      first.SetX((double)(tmp.x()));
      first.SetY((double)tmp.y());
      first.SetZ((double)tmp.z());
      
      tmp = pair->Track2()->Track()->NominalTpcPoint(i);
      second.SetX((double)tmp.x());
      second.SetY((double)tmp.y());
      second.SetZ((double)tmp.z()); 

      avgSep += TMath::Sqrt(((double)first.x()-(double)second.x())*((double)first.x()-(double)second.x())+((double)first.y()-(double)second.y())*((double)first.y()-second.y())+((double)first.z()-(double)second.z())*((double)first.z()-(double)second.z()));
    }
  avgSep /= 8;

  if(avgSep<fMinAvgSepTrackPos) return false;

  avgSep = 0;
  
  for(int i=0; i<8 ;i++)
    {
      tmp = pair->Track1()->V0()->NominalTpcPointNeg(i);
      //cout<<"X pos: "<<tmp.x()<<endl;
      first.SetX((double)(tmp.x()));
      first.SetY((double)tmp.y());
      first.SetZ((double)tmp.z());
      
      tmp = pair->Track2()->Track()->NominalTpcPoint(i);
      second.SetX((double)tmp.x());
      second.SetY((double)tmp.y());
      second.SetZ((double)tmp.z()); 

      avgSep += TMath::Sqrt(((double)first.x()-(double)second.x())*((double)first.x()-(double)second.x())+((double)first.y()-(double)second.y())*((double)first.y()-second.y())+((double)first.z()-(double)second.z())*((double)first.z()-(double)second.z()));
    }
  avgSep /= 8;

  if(avgSep<fMinAvgSepTrackNeg) return false;


  //Qinv cut (we are trying to get rid of antiresidual correlation between primary protons)
  if(fKstarCut > 0)
    {

      //2 daughters of first V0
      double mom1PosX = pair->Track1()->V0()->MomPosX();
      double mom1PosY = pair->Track1()->V0()->MomPosY();
      double mom1PosZ = pair->Track1()->V0()->MomPosZ();
      double mom1NegX = pair->Track1()->V0()->MomNegX();
      double mom1NegY = pair->Track1()->V0()->MomNegY();
      double mom1NegZ = pair->Track1()->V0()->MomNegZ();


      //double PionMass = 0.13956995;
      //double KaonMass = 0.493677;
      double ProtonMass = 0.938272;
      //double LambdaMass = 1.115683;

      AliFemtoLorentzVector fFourMomentum1; // Particle momentum
      AliFemtoLorentzVector fFourMomentum2; // Particle momentum

      AliFemtoThreeVector temp1;
      double ener1=0;
      if(fFirstParticleType == 0) //lambda
	{
	  if(fSecondParticleType == 2) //proton
	    {
	      //temp1 = ::sqrt(mom1PosX*mom1PosX+mom1PosY*mom1PosY+mom1PosZ*mom1PosZ);
	      temp1.SetX(mom1PosX);
	      temp1.SetY(mom1PosY);
	      temp1.SetZ(mom1PosZ);
	      ener1 = ::sqrt(temp1.Mag2()+ProtonMass*ProtonMass);	  
	    }
	}
      else if(fFirstParticleType == 1) //antilambda
	{
	  if(fSecondParticleType == 3) //antiproton
	    {
	      //temp1 = ::sqrt(mom1NegX*mom1NegX+mom1NegY*mom1NegY+mom1NegZ*mom1NegZ);
	      temp1.SetX(mom1NegX);
	      temp1.SetY(mom1NegY);
	      temp1.SetZ(mom1NegZ);
	      ener1 = ::sqrt(temp1.Mag2()+ProtonMass*ProtonMass);
	    }
	}
      fFourMomentum1.SetVect(temp1);
      fFourMomentum1.SetE(ener1);

      //AliFemtoLorentzVector fFourMomentum2; // Particle momentum
      AliFemtoThreeVector temp2;
      double ener2=0;

      //2 daughters of second V0
      temp2 = pair->Track2()->Track()->P();
 
 
      if(fSecondParticleType == 2 || fSecondParticleType == 3) //proton
	{
	    ener2 = ::sqrt(temp2.Mag2()+ProtonMass*ProtonMass);
	}

      fFourMomentum2.SetVect(temp2);
      fFourMomentum2.SetE(ener2);

      //QInv calculation
      AliFemtoLorentzVector tDiff = (fFourMomentum1-fFourMomentum2);

      double tQinv = fabs(-1.*tDiff.m());   // note - qInv() will be negative for identical pairs...
   
      //cout<<"tQinv/2: "<<tQinv/2<<endl;
      if(tQinv/2 < fKstarCut)
	{

	  temp = false;
	}
    }
  return temp;
}
//__________________
AliFemtoString AliFemtoV0TrackPairCut::Report(){
  // Prepare the report from the execution
  string stemp = "AliFemtoV0 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 AliFemtoV0TrackPairCut::SetV0Max(Double_t aV0Max) {
  fV0Max = aV0Max;
}

Double_t AliFemtoV0TrackPairCut::GetAliFemtoV0Max() const {
  return fV0Max;
}


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

  return tListSetttings;
}

void     AliFemtoV0TrackPairCut::SetRemoveSameLabel(Bool_t aRemove)
{
  fRemoveSameLabel = aRemove;
}

void AliFemtoV0TrackPairCut::SetTPCOnly(Bool_t tpconly)
{
  fTrackTPCOnly = tpconly;
}

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

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

void AliFemtoV0TrackPairCut::SetDataType(AliFemtoDataType type)
{
  fDataType = type;
}

void AliFemtoV0TrackPairCut::SetTPCEntranceSepMinimum(double dtpc)
{
  fDTPCMin = dtpc;
}

void AliFemtoV0TrackPairCut::SetTPCExitSepMinimum(double dtpc)
{
  fDTPCExitMin = dtpc;
}

void AliFemtoV0TrackPairCut::SetKstarCut(double kstar, AliFemtoParticleType firstParticle, AliFemtoParticleType secondParticle)
{
  fKstarCut = kstar; 
  fFirstParticleType = firstParticle;  //for kstar - first particle type (V0 type) 
  fSecondParticleType = secondParticle;
}

void AliFemtoV0TrackPairCut::SetMinAvgSeparation(int type, double minSep)
{
  if(type == 0) //Track-Pos
    fMinAvgSepTrackPos = minSep;
  else if(type == 1) //Track-Neg
    fMinAvgSepTrackNeg = minSep;
}
Commit	Line	Data
973a91f8	1	/////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoShareQualityPairCut - a pair cut which checks for some pair //
	4	// qualities that attempt to identify slit/doubly reconstructed tracks //
	5	// //
	6	/////////////////////////////////////////////////////////////////////////////
	7	/***************************************************************************
	8	*
	9	* $Id: AliFemtoShareQualityPairCut.cxx 50722 2011-07-21 15:18:38Z 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 "AliFemtoV0TrackPairCut.h"
	23	#include <string>
	24	#include <cstdio>
	25
	26	#ifdef __ROOT__
	27	ClassImp(AliFemtoV0TrackPairCut)
	28	#endif
	29
	30	//__________________
	31	AliFemtoV0TrackPairCut::AliFemtoV0TrackPairCut():
	32	fNPairsPassed(0),
	33	fNPairsFailed(0),
	34	fV0Max(1.0),
	35	fShareQualityMax(1.0),
	36	fShareFractionMax(1.0),
	37	fRemoveSameLabel(0),
1ae130d9	38	fTrackTPCOnly(0),
	39	fDataType(kAOD),
	40	fDTPCMin(0),
ce7b3d98	41	fDTPCExitMin(0),
	42	fKstarCut(0),
	43	fFirstParticleType(kLambda),
6691ea4f	44	fSecondParticleType(kProton),
	45	fMinAvgSepTrackPos(0),
	46	fMinAvgSepTrackNeg(0)
973a91f8	47	{
	48	// Default constructor
	49	// Nothing to do
	50	}
	51	//__________________
	52	AliFemtoV0TrackPairCut::~AliFemtoV0TrackPairCut(){
	53	/* no-op */
	54	}
1ae130d9	55	AliFemtoV0TrackPairCut& AliFemtoV0TrackPairCut::operator=(const AliFemtoV0TrackPairCut& cut)
578b62b9	56	{
	57	if (this != &cut) {
	58
	59	AliFemtoPairCut::operator=(cut);
	60	fNPairsPassed = 0;
	61	fNPairsFailed =0;
	62	fV0Max = 1.0;
	63	fShareQualityMax = 1.0;
	64	fShareFractionMax = 1.0;
	65	fRemoveSameLabel = 0;
1ae130d9	66	fTrackTPCOnly = 0;
	67	fDataType = kAOD;
	68	fDTPCMin = 0;
	69	fDTPCExitMin = 0;
6691ea4f	70	fMinAvgSepTrackPos = 0;
6691ea4f	71	fMinAvgSepTrackNeg = 0;
578b62b9	72	}
	73
	74	return *this;
	75	}
973a91f8	76	//__________________
	77	bool AliFemtoV0TrackPairCut::Pass(const AliFemtoPair* pair){
	78	// Check for pairs that are possibly shared/double reconstruction
	79
	80	bool temp = true;
	81	//Track1 - V0
	82	//Track2 - track
	83
973a91f8	84	if(!(pair->Track1()->V0() && pair->Track2()->Track()))
	85	{
	86	return false;
	87	}
	88	if(fTrackTPCOnly)
	89	{
	90	if( (-(pair->Track1()->V0()->IdNeg()+1)) ==pair->Track2()->TrackId() \|\| (-(pair->Track1()->V0()->IdPos()+1)) ==pair->Track2()->TrackId())
	91	{
	92	return false;
	93	}
	94	}
	95	else
	96	{
	97	if(pair->Track1()->V0()->IdNeg()==pair->Track2()->TrackId() \|\| pair->Track1()->V0()->IdPos()==pair->Track2()->TrackId())
	98	{
	99	return false;
	100	}
	101	}
	102
1ae130d9	103	bool tempTPCEntrancePos = true;
	104	bool tempTPCEntranceNeg = true;
	105	bool tempTPCExitPos = true;
	106	bool tempTPCExitNeg = true;
	107	if(fDataType==kESD \|\| fDataType==kAOD)
	108	{
	109	double distx = pair->Track1()->V0()->NominalTpcEntrancePointPos().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
	110	double disty = pair->Track1()->V0()->NominalTpcEntrancePointPos().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
	111	double distz = pair->Track1()->V0()->NominalTpcEntrancePointPos().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
	112	double distPos = sqrt(distxdistx + distydisty + distz*distz);
	113
	114	distx = pair->Track1()->V0()->NominalTpcEntrancePointNeg().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
	115	disty = pair->Track1()->V0()->NominalTpcEntrancePointNeg().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
	116	distz = pair->Track1()->V0()->NominalTpcEntrancePointNeg().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
	117	double distNeg = sqrt(distxdistx + distydisty + distz*distz);
	118
	119	double distExitX = pair->Track1()->V0()->NominalTpcExitPointPos().x() - pair->Track2()->Track()->NominalTpcExitPoint().x();
	120	double distExitY = pair->Track1()->V0()->NominalTpcExitPointPos().y() - pair->Track2()->Track()->NominalTpcExitPoint().y();
	121	double distExitZ = pair->Track1()->V0()->NominalTpcExitPointPos().z() - pair->Track2()->Track()->NominalTpcExitPoint().z();
	122	double distExitPos = sqrt(distExitXdistExitX + distExitYdistExitY + distExitZ*distExitZ);
	123
	124	distExitX = pair->Track1()->V0()->NominalTpcExitPointNeg().x() - pair->Track2()->Track()->NominalTpcExitPoint().x();
	125	distExitY = pair->Track1()->V0()->NominalTpcExitPointNeg().y() - pair->Track2()->Track()->NominalTpcExitPoint().y();
	126	distExitZ = pair->Track1()->V0()->NominalTpcExitPointNeg().z() - pair->Track2()->Track()->NominalTpcExitPoint().z();
	127	double distExitNeg = sqrt(distExitXdistExitX + distExitYdistExitY + distExitZ*distExitZ);
	128
	129	tempTPCEntrancePos = distPos > fDTPCMin;
	130	tempTPCEntranceNeg = distNeg > fDTPCMin;
	131	tempTPCExitPos = distExitPos > fDTPCExitMin;
	132	tempTPCExitNeg = distExitNeg > fDTPCExitMin;
	133	}
	134
	135	if(!tempTPCEntrancePos \|\| !tempTPCEntranceNeg \|\| !tempTPCExitPos \|\| !tempTPCExitNeg) return false;
973a91f8	136
	137	//reject merged trakcs in TPC
	138
	139
	140	//temp = dist > fDTPCMin;
	141	//koniec kopii
	142	//if(!temp)
	143	//return false;
	144
	145	//kopia z AliFemtoShareQualityPairCut.cxx
	146	Int_t nh = 0;
	147	Int_t an = 0;
	148	Int_t ns = 0;
	149
	150
	151
	152	if ((fShareFractionMax < 1.0) && ( fShareQualityMax < 1.0)) {
	153	for (unsigned int imap=0; imap<pair->Track1()->V0()->TPCclustersPos().GetNbits(); imap++) {
	154	// If both have clusters in the same row
	155	if (pair->Track1()->V0()->TPCclustersPos().TestBitNumber(imap) &&
	156	pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
	157	// Do they share it ?
	158	if (pair->Track1()->V0()->TPCsharingPos().TestBitNumber(imap) &&
	159	pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
	160	{
	161	// cout << "A shared cluster !!!" << endl;
	162	// cout << "imap idx1 idx2 "
	163	// << imap << " "
	164	// << tP1idx[imap] << " " << tP2idx[imap] << endl;
	165	an++;
	166	nh+=2;
	167	ns+=2;
	168	}
	169
	170	// Different hits on the same padrow
	171	else {
	172	an--;
	173	nh+=2;
	174	}
	175	}
	176	else if (pair->Track1()->V0()->TPCclustersPos().TestBitNumber(imap) \|\|
	177	pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
	178	// One track has a hit, the other does not
	179	an++;
	180	nh++;
	181	}
	182	}
	183
	184	Float_t hsmval = 0.0;
	185	Float_t hsfval = 0.0;
	186
	187	if (nh >0) {
	188	hsmval = an*1.0/nh;
	189	hsfval = ns*1.0/nh;
	190	}
	191	// if (hsmval > -0.4) {
	192	// cout << "Pair quality: " << hsmval << " " << an << " " << nh << " "
	193	// << (pair->Track1()->Track()) << " "
	194	// << (pair->Track2()->Track()) << endl;
	195	// cout << "Bits: " << pair->Track1()->Track()->TPCclusters().GetNbits() << endl;
	196	// }
	197	// if (hsfval > 0.0) {
	198	// cout << "Pair sharity: " << hsfval << " " << ns << " " << nh << " " << hsmval << " " << an << " " << nh << endl;
	199	// }
200
201	temp = (hsmval < fShareQualityMax) && (hsfval < fShareFractionMax);
202	if(!temp) return false;
203
204	nh = 0;
205	an = 0;
206	ns = 0;
207
208	for (unsigned int imap=0; imap<pair->Track1()->V0()->TPCclustersNeg().GetNbits(); imap++) {
209	// If both have clusters in the same row
210	if (pair->Track1()->V0()->TPCclustersNeg().TestBitNumber(imap) &&
211	pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
212	// Do they share it ?
213	if (pair->Track1()->V0()->TPCsharingNeg().TestBitNumber(imap) &&
214	pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
215	{
216	// cout << "A shared cluster !!!" << endl;
217	// cout << "imap idx1 idx2 "
218	// << imap << " "
219	// << tP1idx[imap] << " " << tP2idx[imap] << endl;
220	an++;
221	nh+=2;
222	ns+=2;
223	}
224
225	// Different hits on the same padrow
226	else {
227	an--;
228	nh+=2;
229	}
230	}
231	else if (pair->Track1()->V0()->TPCclustersNeg().TestBitNumber(imap) \|\|
232	pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
233	// One track has a hit, the other does not
234	an++;
235	nh++;
236	}
237	}
238
239	hsmval = 0.0;
240	hsfval = 0.0;
241
242	if (nh >0) {
243	hsmval = an*1.0/nh;
244	hsfval = ns*1.0/nh;
245	}
246	// if (hsmval > -0.4) {
247	// cout << "Pair quality: " << hsmval << " " << an << " " << nh << " "
248	// << (pair->Track1()->Track()) << " "
249	// << (pair->Track2()->Track()) << endl;
250	// cout << "Bits: " << pair->Track1()->Track()->TPCclusters().GetNbits() << endl;
251	// }
252	// if (hsfval > 0.0) {
253	// cout << "Pair sharity: " << hsfval << " " << ns << " " << nh << " " << hsmval << " " << an << " " << nh << endl;
254	// }
255
256	temp = (hsmval < fShareQualityMax) && (hsfval < fShareFractionMax);
257
258
259	}
260	else
261	temp = true;
262	//koniec kopii
263
264
6691ea4f	265	//avg sep pair cut
	266	double avgSep=0;
	267	AliFemtoThreeVector first, second, tmp;
	268	for(int i=0; i<8 ;i++)
	269	{
	270	tmp = pair->Track1()->V0()->NominalTpcPointPos(i);
	271	//cout<<"X pos: "<<tmp.x()<<endl;
	272	first.SetX((double)(tmp.x()));
	273	first.SetY((double)tmp.y());
	274	first.SetZ((double)tmp.z());
	275
	276	tmp = pair->Track2()->Track()->NominalTpcPoint(i);
	277	second.SetX((double)tmp.x());
	278	second.SetY((double)tmp.y());
	279	second.SetZ((double)tmp.z());
	280
	281	avgSep += TMath::Sqrt(((double)first.x()-(double)second.x())((double)first.x()-(double)second.x())+((double)first.y()-(double)second.y())((double)first.y()-second.y())+((double)first.z()-(double)second.z())*((double)first.z()-(double)second.z()));
	282	}
	283	avgSep /= 8;
	284
	285	if(avgSep<fMinAvgSepTrackPos) return false;
	286
	287	avgSep = 0;
	288
	289	for(int i=0; i<8 ;i++)
	290	{
	291	tmp = pair->Track1()->V0()->NominalTpcPointNeg(i);
	292	//cout<<"X pos: "<<tmp.x()<<endl;
	293	first.SetX((double)(tmp.x()));
	294	first.SetY((double)tmp.y());
	295	first.SetZ((double)tmp.z());
	296
	297	tmp = pair->Track2()->Track()->NominalTpcPoint(i);
	298	second.SetX((double)tmp.x());
	299	second.SetY((double)tmp.y());
	300	second.SetZ((double)tmp.z());
	301
	302	avgSep += TMath::Sqrt(((double)first.x()-(double)second.x())((double)first.x()-(double)second.x())+((double)first.y()-(double)second.y())((double)first.y()-second.y())+((double)first.z()-(double)second.z())*((double)first.z()-(double)second.z()));
	303	}
	304	avgSep /= 8;
	305
	306	if(avgSep<fMinAvgSepTrackNeg) return false;
	307
	308
	309
	310
ce7b3d98	311	//Qinv cut (we are trying to get rid of antiresidual correlation between primary protons)
	312	if(fKstarCut > 0)
	313	{
	314
	315	//2 daughters of first V0
	316	double mom1PosX = pair->Track1()->V0()->MomPosX();
	317	double mom1PosY = pair->Track1()->V0()->MomPosY();
	318	double mom1PosZ = pair->Track1()->V0()->MomPosZ();
	319	double mom1NegX = pair->Track1()->V0()->MomNegX();
	320	double mom1NegY = pair->Track1()->V0()->MomNegY();
	321	double mom1NegZ = pair->Track1()->V0()->MomNegZ();
	322
	323
	324	//double PionMass = 0.13956995;
	325	//double KaonMass = 0.493677;
	326	double ProtonMass = 0.938272;
	327	//double LambdaMass = 1.115683;
	328
	329	AliFemtoLorentzVector fFourMomentum1; // Particle momentum
	330	AliFemtoLorentzVector fFourMomentum2; // Particle momentum
	331
	332	AliFemtoThreeVector temp1;
	333	double ener1=0;
	334	if(fFirstParticleType == 0) //lambda
	335	{
	336	if(fSecondParticleType == 2) //proton
	337	{
	338	//temp1 = ::sqrt(mom1PosXmom1PosX+mom1PosYmom1PosY+mom1PosZ*mom1PosZ);
	339	temp1.SetX(mom1PosX);
	340	temp1.SetY(mom1PosY);
	341	temp1.SetZ(mom1PosZ);
	342	ener1 = ::sqrt(temp1.Mag2()+ProtonMass*ProtonMass);
	343	}
	344	}
	345	else if(fFirstParticleType == 1) //antilambda
	346	{
	347	if(fSecondParticleType == 3) //antiproton
	348	{
	349	//temp1 = ::sqrt(mom1NegXmom1NegX+mom1NegYmom1NegY+mom1NegZ*mom1NegZ);
	350	temp1.SetX(mom1NegX);
	351	temp1.SetY(mom1NegY);
	352	temp1.SetZ(mom1NegZ);
	353	ener1 = ::sqrt(temp1.Mag2()+ProtonMass*ProtonMass);
	354	}
	355	}
	356	fFourMomentum1.SetVect(temp1);
	357	fFourMomentum1.SetE(ener1);
	358
	359	//AliFemtoLorentzVector fFourMomentum2; // Particle momentum
	360	AliFemtoThreeVector temp2;
	361	double ener2=0;
	362
	363	//2 daughters of second V0
	364	temp2 = pair->Track2()->Track()->P();
	365
	366
	367	if(fSecondParticleType == 2 \|\| fSecondParticleType == 3) //proton
	368	{
	369	ener2 = ::sqrt(temp2.Mag2()+ProtonMass*ProtonMass);
	370	}
	371
	372	fFourMomentum2.SetVect(temp2);
	373	fFourMomentum2.SetE(ener2);
973a91f8	374
ce7b3d98	375	//QInv calculation
ce7b3d98	376	AliFemtoLorentzVector tDiff = (fFourMomentum1-fFourMomentum2);
973a91f8	377
ce7b3d98	378	double tQinv = fabs(-1.*tDiff.m()); // note - qInv() will be negative for identical pairs...
	379
	380	//cout<<"tQinv/2: "<<tQinv/2<<endl;
	381	if(tQinv/2 < fKstarCut)
	382	{
973a91f8	383
ce7b3d98	384	temp = false;
	385	}
	386	}
973a91f8	387	return temp;
	388	}
	389	//__________________
	390	AliFemtoString AliFemtoV0TrackPairCut::Report(){
	391	// Prepare the report from the execution
	392	string stemp = "AliFemtoV0 Pair Cut - remove shared and split pairs\n"; char ctemp[100];
	393	snprintf(ctemp , 100, "Number of pairs which passed:\t%ld Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed);
	394	stemp += ctemp;
	395	AliFemtoString returnThis = stemp;
	396	return returnThis;}
	397	//__________________
	398
	399	void AliFemtoV0TrackPairCut::SetV0Max(Double_t aV0Max) {
	400	fV0Max = aV0Max;
	401	}
	402
	403	Double_t AliFemtoV0TrackPairCut::GetAliFemtoV0Max() const {
	404	return fV0Max;
	405	}
	406
	407
	408	TList *AliFemtoV0TrackPairCut::ListSettings()
	409	{
	410	// return a list of settings in a writable form
	411	TList *tListSetttings = new TList();
	412	char buf[200];
	413	snprintf(buf, 200, "AliFemtoV0TrackPairCut.sharequalitymax=%f", fV0Max);
	414	snprintf(buf, 200, "AliFemtoV0TrackPairCut.sharefractionmax=%f", fShareFractionMax);
	415	tListSetttings->AddLast(new TObjString(buf));
	416
	417	return tListSetttings;
	418	}
	419
	420	void AliFemtoV0TrackPairCut::SetRemoveSameLabel(Bool_t aRemove)
	421	{
	422	fRemoveSameLabel = aRemove;
	423	}
	424
	425	void AliFemtoV0TrackPairCut::SetTPCOnly(Bool_t tpconly)
	426	{
	427	fTrackTPCOnly = tpconly;
	428	}
	429
	430	void AliFemtoV0TrackPairCut::SetShareQualityMax(Double_t aShareQualityMax) {
	431	fShareQualityMax = aShareQualityMax;
	432	}
	433
	434	void AliFemtoV0TrackPairCut::SetShareFractionMax(Double_t aShareFractionMax) {
	435	fShareFractionMax = aShareFractionMax;
	436	}
1ae130d9	437
	438	void AliFemtoV0TrackPairCut::SetDataType(AliFemtoDataType type)
	439	{
	440	fDataType = type;
	441	}
	442
	443	void AliFemtoV0TrackPairCut::SetTPCEntranceSepMinimum(double dtpc)
	444	{
	445	fDTPCMin = dtpc;
	446	}
	447
	448	void AliFemtoV0TrackPairCut::SetTPCExitSepMinimum(double dtpc)
	449	{
	450	fDTPCExitMin = dtpc;
	451	}
ce7b3d98	452
	453	void AliFemtoV0TrackPairCut::SetKstarCut(double kstar, AliFemtoParticleType firstParticle, AliFemtoParticleType secondParticle)
	454	{
	455	fKstarCut = kstar;
	456	fFirstParticleType = firstParticle; //for kstar - first particle type (V0 type)
	457	fSecondParticleType = secondParticle;
	458	}
6691ea4f	459
	460	void AliFemtoV0TrackPairCut::SetMinAvgSeparation(int type, double minSep)
	461	{
	462	if(type == 0) //Track-Pos
	463	fMinAvgSepTrackPos = minSep;
	464	else if(type == 1) //Track-Neg
	465	fMinAvgSepTrackNeg = minSep;
	466	}