[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / AliRsnCutTrackQuality.cxx

//
// Class AliRsnCutTrackQuality
//
// General implementation of a single cut strategy, which can be:
// - a value contained in a given interval  [--> IsBetween()   ]
// - a value equal to a given reference     [--> MatchesValue()]
//
// In all cases, the reference value(s) is (are) given as data members
// and each kind of cut requires a given value type (Int, UInt, Double),
// but the cut check procedure is then automatized and chosen thanks to
// an enumeration of the implemented cut types.
// At the end, the user (or any other point which uses this object) has
// to use the method IsSelected() to check if this cut has been passed.
//
// authors: Martin Vala (martin.vala@cern.ch)
//          Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
//

#include <Riostream.h>
#include <TFormula.h>
#include <TBits.h>

#include "AliLog.h"
#include "AliESDtrackCuts.h"

#include "AliRsnEvent.h"
#include "AliRsnDaughter.h"
#include "AliRsnCutTrackQuality.h"

ClassImp(AliRsnCutTrackQuality)

//_________________________________________________________________________________________________
AliRsnCutTrackQuality::AliRsnCutTrackQuality(const char *name) :
   AliRsnCut(name, AliRsnTarget::kDaughter, 0.0, 0.0),
   fFlagsOn(0x0),
   fFlagsOff(0x0),
   fRejectKinkDaughters(kTRUE),
   fDCARfixed(kTRUE),
   fDCARptFormula(""),
   fDCARmax(1E20),
   fDCAZfixed(kTRUE),
   fDCAZptFormula(""),
   fDCAZmax(1E20),
   fSPDminNClusters(0),
   fITSminNClusters(0),
   fITSmaxChi2(1E20),
   fTPCminNClusters(0),
   fTPCmaxChi2(1E20),
   fCutMaxChi2TPCConstrainedVsGlobal(1E20),
   fAODTestFilterBit(-1),
   fESDtrackCuts(0x0)
{
//
// Default constructor.
// Initializes all cuts in such a way that all of them are disabled.
//

   SetPtRange(0.0, 1E20);
   SetEtaRange(-1E20, 1E20);
}

//_________________________________________________________________________________________________
AliRsnCutTrackQuality::AliRsnCutTrackQuality(const AliRsnCutTrackQuality &copy) :
   AliRsnCut(copy),
   fFlagsOn(copy.fFlagsOn),
   fFlagsOff(copy.fFlagsOff),
   fRejectKinkDaughters(copy.fRejectKinkDaughters),
   fDCARfixed(copy.fDCARfixed),
   fDCARptFormula(copy.fDCARptFormula),
   fDCARmax(copy.fDCARmax),
   fDCAZfixed(copy.fDCAZfixed),
   fDCAZptFormula(copy.fDCAZptFormula),
   fDCAZmax(copy.fDCAZmax),
   fSPDminNClusters(copy.fSPDminNClusters),
   fITSminNClusters(copy.fITSminNClusters),
   fITSmaxChi2(copy.fITSmaxChi2),
   fTPCminNClusters(copy.fTPCminNClusters),
   fTPCmaxChi2(copy.fTPCmaxChi2),
   fCutMaxChi2TPCConstrainedVsGlobal(copy.fCutMaxChi2TPCConstrainedVsGlobal),
   fAODTestFilterBit(copy.fAODTestFilterBit),
   fESDtrackCuts(copy.fESDtrackCuts)
{
//
// Copy constructor.
// Just copy all data member values.
//

   SetPtRange(copy.fPt[0], copy.fPt[1]);
   SetEtaRange(copy.fEta[0], copy.fEta[1]);
}

//_________________________________________________________________________________________________
AliRsnCutTrackQuality &AliRsnCutTrackQuality::operator=(const AliRsnCutTrackQuality &copy)
{
//
// Assignment operator.
// Just copy all data member values.
//

   if (this == &copy)
      return *this;

   fFlagsOn = copy.fFlagsOn;
   fFlagsOff = copy.fFlagsOff;
   fRejectKinkDaughters = copy.fRejectKinkDaughters;
   fDCARfixed = copy.fDCARfixed;
   fDCARptFormula = copy.fDCARptFormula;
   fDCARmax = copy.fDCARmax;
   fDCAZfixed = copy.fDCAZfixed;
   fDCAZptFormula = copy.fDCAZptFormula;
   fDCAZmax = copy.fDCAZmax;
   fSPDminNClusters = copy.fSPDminNClusters;
   fITSminNClusters = copy.fITSminNClusters;
   fITSmaxChi2 = copy.fITSmaxChi2;
   fTPCminNClusters = copy.fTPCminNClusters;
   fTPCmaxChi2 = copy.fTPCmaxChi2;
   fAODTestFilterBit = copy.fAODTestFilterBit;
   fESDtrackCuts = copy.fESDtrackCuts;
   SetPtRange(copy.fPt[0], copy.fPt[1]);
   SetEtaRange(copy.fEta[0], copy.fEta[1]);

   return (*this);
}

//_________________________________________________________________________________________________
void AliRsnCutTrackQuality::DisableAll()
{
//
// Disable all cuts
//

   fFlagsOn = 0x0;
   fFlagsOff = 0x0;
   fRejectKinkDaughters = kFALSE;
   fDCARfixed = kTRUE;
   fDCARptFormula = "";
   fDCARmax = 1E20;
   fDCAZfixed = kTRUE;
   fDCAZptFormula = "";
   fDCAZmax = 1E20;
   fSPDminNClusters = 0;
   fITSminNClusters = 0;
   fITSmaxChi2 = 1E20;
   fTPCminNClusters = 0;
   fTPCmaxChi2 = 1E20;
   fAODTestFilterBit = -1;
   if (fESDtrackCuts) {
      const char *cutsName = fESDtrackCuts->GetName();
      const char *cutsTitle = fESDtrackCuts->GetTitle();
      delete fESDtrackCuts;
      fESDtrackCuts = new AliESDtrackCuts(cutsName,cutsTitle);
   }
   SetPtRange(0.0, 1E20);
   SetEtaRange(-1E20, 1E20);
}

//_________________________________________________________________________________________________
Bool_t AliRsnCutTrackQuality::IsSelected(TObject *object)
{
//
// Cut checker.
// Checks the type of object being evaluated
// and then calls the appropriate sub-function (for ESD or AOD)
//

   // coherence check
   if (!TargetOK(object)) return kFALSE;

   // status is checked in the same way for all tracks, using AliVTrack
   // as a convention, if a the collection of 'on' flags is '0x0', it
   // is assumed that no flags are required, and this check is skipped;
   // for the collection of 'off' flags this is not needed
   AliVTrack *vtrack = fDaughter->Ref2Vtrack();
   if (!vtrack) {
      AliDebug(AliLog::kDebug + 2, "This object is not either an ESD nor AOD track");
      return kFALSE;
   }
   ULong_t status   = (ULong_t)vtrack->GetStatus();
   ULong_t checkOn  = status & fFlagsOn;
   ULong_t checkOff = status & fFlagsOff;
   if (fFlagsOn != 0x0 && checkOn != fFlagsOn) {
      AliDebug(AliLog::kDebug + 2, Form("Failed flag check: required  %s", Binary(fFlagsOn)));
      AliDebug(AliLog::kDebug + 2, Form("                   track has %s", Binary(status  )));
      return kFALSE;
   }
   if (checkOff != 0) {
      AliDebug(AliLog::kDebug + 2, Form("Failed flag check: forbidden %s", Binary(fFlagsOff)));
      AliDebug(AliLog::kDebug + 2, Form("                   track has %s", Binary(status  )));
      return kFALSE;
   }
   AliDebug(AliLog::kDebug + 3, Form("Flag check OK: required  %s", Binary(fFlagsOn)));
   AliDebug(AliLog::kDebug + 3, Form("               forbidden %s", Binary(fFlagsOff)));
   AliDebug(AliLog::kDebug + 3, Form("               track has %s", Binary(status  )));

   // retrieve real object type
   AliESDtrack *esdTrack = fDaughter->Ref2ESDtrack();
   AliAODTrack *aodTrack = fDaughter->Ref2AODtrack();
   if (esdTrack) {
      AliDebug(AliLog::kDebug + 2, "Checking an ESD track");
      if (fESDtrackCuts)
         return fESDtrackCuts->IsSelected(esdTrack);
      else
         return CheckESD(esdTrack);
   } else if (aodTrack) {
      AliDebug(AliLog::kDebug + 2, "Checking an AOD track");
      return CheckAOD(aodTrack);
   } else {
      AliDebug(AliLog::kDebug + 2, Form("This object is not either an ESD nor AOD track, it is an %s", fDaughter->GetRef()->ClassName()));
      return kFALSE;
   }
}

//_________________________________________________________________________________________________
Bool_t AliRsnCutTrackQuality::CheckESD(AliESDtrack *track)
{
//
// Check an ESD track.
// This is done using the default track checker for ESD.
// It is declared static, not to recreate it every time.
//

   static AliESDtrackCuts cuts;

   // general acceptance/pt cuts
   cuts.SetPtRange(fPt[0], fPt[1]);
   cuts.SetEtaRange(fEta[0], fEta[1]);

   // transverse DCA cuts
   if (fDCARfixed)
      cuts.SetMaxDCAToVertexXY(fDCARmax);
   else
      cuts.SetMaxDCAToVertexXYPtDep(fDCARptFormula.Data());

   // longitudinal DCA cuts
   if (fDCAZfixed)
      cuts.SetMaxDCAToVertexZ(fDCAZmax);
   else
      cuts.SetMaxDCAToVertexZPtDep(fDCAZptFormula.Data());

   // these options are always disabled in current version
   cuts.SetDCAToVertex2D(kFALSE);
   cuts.SetRequireSigmaToVertex(kFALSE);

   // TPC related cuts for TPC+ITS tracks
   cuts.SetMinNClustersTPC(fTPCminNClusters);
   cuts.SetMaxChi2PerClusterTPC(fTPCmaxChi2);
   cuts.SetAcceptKinkDaughters(!fRejectKinkDaughters);
   cuts.SetMaxChi2TPCConstrainedGlobal(fCutMaxChi2TPCConstrainedVsGlobal);

   // ITS related cuts for TPC+ITS tracks
   if (fSPDminNClusters > 0)
      cuts.SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
   cuts.SetMaxChi2PerClusterITS(fITSmaxChi2);

   // now that all is initialized, do the check
   return cuts.IsSelected(track);
}

//_________________________________________________________________________________________________
Bool_t AliRsnCutTrackQuality::CheckAOD(AliAODTrack *track)
{
//
// Check an AOD track.
// This is done doing directly all checks, since there is not
// an equivalend checker for AOD tracks
//

   // if a test bit is used, check it and skip the following
   if (fAODTestFilterBit >= 0) {
      UInt_t bit = 1 << fAODTestFilterBit;
      AliDebugClass(2, Form("Required a test filter bit for AOD check: %u (result: %s)", bit, (track->TestFilterBit(bit) ? "accept" : "reject")));
      if (!track->TestFilterBit(bit))
         return kFALSE;
      else {
         if (track->Pt() < fPt[0] || track->Pt() > fPt[1]) return kFALSE;
         if (track->Eta() < fEta[0] || track->Eta() > fEta[1]) return kFALSE;
         return kTRUE;
      }
   }

   // try to retrieve the reference AOD event
   AliAODEvent *aodEvent = 0x0;
   if (fEvent) aodEvent = fEvent->GetRefAOD();
   if (!aodEvent) {
      AliError("AOD reference event is not initialized!");
      return kFALSE;
   }

   // step #0: check SPD and ITS clusters
   Int_t nSPD = 0;
   nSPD  = TESTBIT(track->GetITSClusterMap(), 0);
   nSPD += TESTBIT(track->GetITSClusterMap(), 1);
   if (nSPD < fSPDminNClusters) {
      AliDebug(AliLog::kDebug + 2, "Not enough SPD clusters in this track. Rejected");
      return kFALSE;
   }


   // step #1: check number of clusters in TPC
   if (track->GetTPCNcls() < fTPCminNClusters) {
      AliDebug(AliLog::kDebug + 2, "Too few TPC clusters. Rejected");
      return kFALSE;
   }
   if (track->GetITSNcls() < fITSminNClusters) {
      AliDebug(AliLog::kDebug + 2, "Too few ITS clusters. Rejected");
      return kFALSE;
   }

   // step #2: check chi square
   if (track->Chi2perNDF() > fTPCmaxChi2) {
      AliDebug(AliLog::kDebug + 2, "Bad chi2. Rejected");
      return kFALSE;
   }
   if (track->Chi2perNDF() > fITSmaxChi2) {
      AliDebug(AliLog::kDebug + 2, "Bad chi2. Rejected");
      return kFALSE;
   }

   // step #3: reject kink daughters
   AliAODVertex *vertex = track->GetProdVertex();
   if (vertex && fRejectKinkDaughters) {
      if (vertex->GetType() == AliAODVertex::kKink) {
         AliDebug(AliLog::kDebug + 2, "Kink daughter. Rejected");
         return kFALSE;
      }
   }

   // step #4: DCA cut (transverse)
   // --> reject all tracks not ITS refitted
   Double_t b[2], cov[3];
   vertex = aodEvent->GetPrimaryVertex();
   if (!vertex) {
      AliDebug(AliLog::kDebug + 2, "NULL vertex");
      return kFALSE;
   }
   if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
      AliDebug(AliLog::kDebug + 2, "Not ITS refitted");
      return kFALSE;
   }
   if (!track->PropagateToDCA(vertex, aodEvent->GetMagneticField(), kVeryBig, b, cov)) {
      AliDebug(AliLog::kDebug + 2, "Failed propagation to vertex");
      return kFALSE;
   }
   // if the DCA cut is not fixed, compute current value
   if (!fDCARfixed) {
      static TString str(fDCARptFormula);
      str.ReplaceAll("pt", "x");
      static const TFormula dcaXY(Form("%s_dcaXY", GetName()), str.Data());
      fDCARmax = dcaXY.Eval(track->Pt());
   }
   // check the cut
   if (TMath::Abs(b[0]) > fDCARmax) {
      AliDebug(AliLog::kDebug + 2, "Too large transverse DCA");
      return kFALSE;
   }

   // step #5: DCA cut (longitudinal)
   // the DCA has already been computed above
   // if the DCA cut is not fixed, compute current value
   if (!fDCAZfixed) {
      static TString str(fDCAZptFormula);
      str.ReplaceAll("pt", "x");
      static const TFormula dcaZ(Form("%s_dcaXY", GetName()), str.Data());
      fDCAZmax = dcaZ.Eval(track->Pt());
   }
   // check the cut
   if (TMath::Abs(b[1]) > fDCAZmax) {
      AliDebug(AliLog::kDebug + 2, "Too large longitudinal DCA");
      return kFALSE;
   }

   // step #6: check eta/pt range
   if (track->Eta() < fEta[0] || track->Eta() > fEta[1]) {
      AliDebug(AliLog::kDebug + 2, "Outside ETA acceptance");
      return kFALSE;
   }
   if (track->Pt() < fPt[0] || track->Pt() > fPt[1]) {
      AliDebug(AliLog::kDebug + 2, "Outside PT acceptance");
      return kFALSE;
   }

   // if we are here, all cuts were passed and no exit point was got
   AliDebug(AliLog::kDebug + 2, "============================= ACCEPTED TRACK =====================================================");
   return kTRUE;
}

//_________________________________________________________________________________________________
void AliRsnCutTrackQuality::Print(const Option_t *) const
{
//
// Print information on this cut
//

   AliInfo(Form("Cut name                : %s", GetName()));
   AliInfo(Form("Required flags (off, on): %lx %lx", fFlagsOn, fFlagsOff));
   AliInfo(Form("Ranges in eta, pt       : %.2f - %.2f, %.2f - %.2f", fEta[0], fEta[1], fPt[0], fPt[1]));
   AliInfo(Form("Kink daughters are      : %s", (fRejectKinkDaughters ? "rejected" : "accepted")));
   AliInfo(Form("TPC requirements        : min. cluster = %d, max chi2 = %f", fTPCminNClusters, fTPCmaxChi2));
   AliInfo(Form("ITS requirements        : min. cluster = %d (all), %d (SPD), max chi2 = %f", fITSminNClusters, fSPDminNClusters, fITSmaxChi2));

   if (fDCARfixed) {
      AliInfo(Form("DCA r cut               : fixed to %f cm", fDCARmax));
   } else {
      AliInfo(Form("DCA r cut formula       : %s", fDCARptFormula.Data()));
   }

   if (fDCAZfixed) {
      AliInfo(Form("DCA z cut               : fixed to %f cm", fDCAZmax));
   } else {
      AliInfo(Form("DCA z cut formula       : %s", fDCAZptFormula.Data()));
   }
}
//__________________________________________________________________________________________________
void AliRsnCutTrackQuality::SetDefaults2010()
{
//
// Default settings for cuts used in 2010
//
   AddStatusFlag(AliESDtrack::kTPCin   , kTRUE);
   AddStatusFlag(AliESDtrack::kTPCrefit, kTRUE);
   AddStatusFlag(AliESDtrack::kITSrefit, kTRUE);
   SetPtRange(0.15, 1E+20);
   SetEtaRange(-0.8, 0.8);
   SetDCARPtFormula("0.0182+0.0350/pt^1.01");
   SetDCAZmax(2.0);
   SetSPDminNClusters(1);
   SetITSminNClusters(0);
   SetITSmaxChi2(36);
   SetMaxChi2TPCConstrainedGlobal(36);
   SetTPCminNClusters(70);
   SetTPCmaxChi2(4.0);
   SetRejectKinkDaughters();
   SetAODTestFilterBit(5);
}

//__________________________________________________________________________________________________
const char *AliRsnCutTrackQuality::Binary(UInt_t number)
{
//
// Convert an integer in binary
//

   static char b[50];
   b[0] = '\0';

   UInt_t z;
   for (z = 512; z > 0; z >>= 1)
      strncat(b, ((number & z) == z) ? "1" : "0", 1);

   return b;
}
Commit	Line	Data
	1	//
	2	// Class AliRsnCutTrackQuality
	3	//
	4	// General implementation of a single cut strategy, which can be:
	5	// - a value contained in a given interval [--> IsBetween() ]
	6	// - a value equal to a given reference [--> MatchesValue()]
	7	//
	8	// In all cases, the reference value(s) is (are) given as data members
	9	// and each kind of cut requires a given value type (Int, UInt, Double),
	10	// but the cut check procedure is then automatized and chosen thanks to
	11	// an enumeration of the implemented cut types.
	12	// At the end, the user (or any other point which uses this object) has
	13	// to use the method IsSelected() to check if this cut has been passed.
	14	//
	15	// authors: Martin Vala (martin.vala@cern.ch)
	16	// Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
	17	//
	18
	19	#include <Riostream.h>
	20	#include <TFormula.h>
	21	#include <TBits.h>
	22
	23	#include "AliLog.h"
	24	#include "AliESDtrackCuts.h"
	25
	26	#include "AliRsnEvent.h"
	27	#include "AliRsnDaughter.h"
	28	#include "AliRsnCutTrackQuality.h"
	29
	30	ClassImp(AliRsnCutTrackQuality)
	31
	32	//_________________________________________________________________________________________________
	33	AliRsnCutTrackQuality::AliRsnCutTrackQuality(const char *name) :
	34	AliRsnCut(name, AliRsnTarget::kDaughter, 0.0, 0.0),
	35	fFlagsOn(0x0),
	36	fFlagsOff(0x0),
	37	fRejectKinkDaughters(kTRUE),
	38	fDCARfixed(kTRUE),
	39	fDCARptFormula(""),
	40	fDCARmax(1E20),
	41	fDCAZfixed(kTRUE),
	42	fDCAZptFormula(""),
	43	fDCAZmax(1E20),
	44	fSPDminNClusters(0),
	45	fITSminNClusters(0),
	46	fITSmaxChi2(1E20),
	47	fTPCminNClusters(0),
	48	fTPCmaxChi2(1E20),
	49	fCutMaxChi2TPCConstrainedVsGlobal(1E20),
	50	fAODTestFilterBit(-1),
	51	fESDtrackCuts(0x0)
	52	{
	53	//
	54	// Default constructor.
	55	// Initializes all cuts in such a way that all of them are disabled.
	56	//
	57
	58	SetPtRange(0.0, 1E20);
	59	SetEtaRange(-1E20, 1E20);
	60	}
	61
	62	//_________________________________________________________________________________________________
	63	AliRsnCutTrackQuality::AliRsnCutTrackQuality(const AliRsnCutTrackQuality &copy) :
	64	AliRsnCut(copy),
	65	fFlagsOn(copy.fFlagsOn),
	66	fFlagsOff(copy.fFlagsOff),
	67	fRejectKinkDaughters(copy.fRejectKinkDaughters),
	68	fDCARfixed(copy.fDCARfixed),
	69	fDCARptFormula(copy.fDCARptFormula),
	70	fDCARmax(copy.fDCARmax),
	71	fDCAZfixed(copy.fDCAZfixed),
	72	fDCAZptFormula(copy.fDCAZptFormula),
	73	fDCAZmax(copy.fDCAZmax),
	74	fSPDminNClusters(copy.fSPDminNClusters),
	75	fITSminNClusters(copy.fITSminNClusters),
	76	fITSmaxChi2(copy.fITSmaxChi2),
	77	fTPCminNClusters(copy.fTPCminNClusters),
	78	fTPCmaxChi2(copy.fTPCmaxChi2),
	79	fCutMaxChi2TPCConstrainedVsGlobal(copy.fCutMaxChi2TPCConstrainedVsGlobal),
	80	fAODTestFilterBit(copy.fAODTestFilterBit),
	81	fESDtrackCuts(copy.fESDtrackCuts)
	82	{
	83	//
	84	// Copy constructor.
	85	// Just copy all data member values.
	86	//
	87
	88	SetPtRange(copy.fPt[0], copy.fPt[1]);
	89	SetEtaRange(copy.fEta[0], copy.fEta[1]);
	90	}
	91
	92	//_________________________________________________________________________________________________
	93	AliRsnCutTrackQuality &AliRsnCutTrackQuality::operator=(const AliRsnCutTrackQuality &copy)
	94	{
	95	//
	96	// Assignment operator.
	97	// Just copy all data member values.
	98	//
	99
	100	if (this == &copy)
	101	return *this;
	102
	103	fFlagsOn = copy.fFlagsOn;
	104	fFlagsOff = copy.fFlagsOff;
	105	fRejectKinkDaughters = copy.fRejectKinkDaughters;
	106	fDCARfixed = copy.fDCARfixed;
	107	fDCARptFormula = copy.fDCARptFormula;
	108	fDCARmax = copy.fDCARmax;
	109	fDCAZfixed = copy.fDCAZfixed;
	110	fDCAZptFormula = copy.fDCAZptFormula;
	111	fDCAZmax = copy.fDCAZmax;
	112	fSPDminNClusters = copy.fSPDminNClusters;
	113	fITSminNClusters = copy.fITSminNClusters;
	114	fITSmaxChi2 = copy.fITSmaxChi2;
	115	fTPCminNClusters = copy.fTPCminNClusters;
	116	fTPCmaxChi2 = copy.fTPCmaxChi2;
	117	fAODTestFilterBit = copy.fAODTestFilterBit;
	118	fESDtrackCuts = copy.fESDtrackCuts;
	119	SetPtRange(copy.fPt[0], copy.fPt[1]);
	120	SetEtaRange(copy.fEta[0], copy.fEta[1]);
	121
	122	return (*this);
	123	}
	124
	125	//_________________________________________________________________________________________________
	126	void AliRsnCutTrackQuality::DisableAll()
	127	{
	128	//
	129	// Disable all cuts
	130	//
	131
	132	fFlagsOn = 0x0;
	133	fFlagsOff = 0x0;
	134	fRejectKinkDaughters = kFALSE;
	135	fDCARfixed = kTRUE;
	136	fDCARptFormula = "";
	137	fDCARmax = 1E20;
	138	fDCAZfixed = kTRUE;
	139	fDCAZptFormula = "";
	140	fDCAZmax = 1E20;
	141	fSPDminNClusters = 0;
	142	fITSminNClusters = 0;
	143	fITSmaxChi2 = 1E20;
	144	fTPCminNClusters = 0;
	145	fTPCmaxChi2 = 1E20;
	146	fAODTestFilterBit = -1;
	147	if (fESDtrackCuts) {
	148	const char *cutsName = fESDtrackCuts->GetName();
	149	const char *cutsTitle = fESDtrackCuts->GetTitle();
	150	delete fESDtrackCuts;
	151	fESDtrackCuts = new AliESDtrackCuts(cutsName,cutsTitle);
	152	}
	153	SetPtRange(0.0, 1E20);
	154	SetEtaRange(-1E20, 1E20);
	155	}
	156
	157	//_________________________________________________________________________________________________
	158	Bool_t AliRsnCutTrackQuality::IsSelected(TObject *object)
	159	{
	160	//
	161	// Cut checker.
	162	// Checks the type of object being evaluated
	163	// and then calls the appropriate sub-function (for ESD or AOD)
	164	//
	165
	166	// coherence check
	167	if (!TargetOK(object)) return kFALSE;
	168
	169	// status is checked in the same way for all tracks, using AliVTrack
	170	// as a convention, if a the collection of 'on' flags is '0x0', it
	171	// is assumed that no flags are required, and this check is skipped;
	172	// for the collection of 'off' flags this is not needed
	173	AliVTrack *vtrack = fDaughter->Ref2Vtrack();
	174	if (!vtrack) {
	175	AliDebug(AliLog::kDebug + 2, "This object is not either an ESD nor AOD track");
	176	return kFALSE;
	177	}
	178	ULong_t status = (ULong_t)vtrack->GetStatus();
	179	ULong_t checkOn = status & fFlagsOn;
	180	ULong_t checkOff = status & fFlagsOff;
	181	if (fFlagsOn != 0x0 && checkOn != fFlagsOn) {
	182	AliDebug(AliLog::kDebug + 2, Form("Failed flag check: required %s", Binary(fFlagsOn)));
	183	AliDebug(AliLog::kDebug + 2, Form(" track has %s", Binary(status )));
	184	return kFALSE;
	185	}
	186	if (checkOff != 0) {
	187	AliDebug(AliLog::kDebug + 2, Form("Failed flag check: forbidden %s", Binary(fFlagsOff)));
	188	AliDebug(AliLog::kDebug + 2, Form(" track has %s", Binary(status )));
	189	return kFALSE;
	190	}
	191	AliDebug(AliLog::kDebug + 3, Form("Flag check OK: required %s", Binary(fFlagsOn)));
	192	AliDebug(AliLog::kDebug + 3, Form(" forbidden %s", Binary(fFlagsOff)));
	193	AliDebug(AliLog::kDebug + 3, Form(" track has %s", Binary(status )));
	194
	195	// retrieve real object type
	196	AliESDtrack *esdTrack = fDaughter->Ref2ESDtrack();
	197	AliAODTrack *aodTrack = fDaughter->Ref2AODtrack();
	198	if (esdTrack) {
	199	AliDebug(AliLog::kDebug + 2, "Checking an ESD track");
	200	if (fESDtrackCuts)
	201	return fESDtrackCuts->IsSelected(esdTrack);
	202	else
	203	return CheckESD(esdTrack);
	204	} else if (aodTrack) {
	205	AliDebug(AliLog::kDebug + 2, "Checking an AOD track");
	206	return CheckAOD(aodTrack);
	207	} else {
	208	AliDebug(AliLog::kDebug + 2, Form("This object is not either an ESD nor AOD track, it is an %s", fDaughter->GetRef()->ClassName()));
	209	return kFALSE;
	210	}
	211	}
	212
	213	//_________________________________________________________________________________________________
	214	Bool_t AliRsnCutTrackQuality::CheckESD(AliESDtrack *track)
	215	{
	216	//
	217	// Check an ESD track.
	218	// This is done using the default track checker for ESD.
	219	// It is declared static, not to recreate it every time.
	220	//
	221
	222	static AliESDtrackCuts cuts;
	223
	224	// general acceptance/pt cuts
	225	cuts.SetPtRange(fPt[0], fPt[1]);
	226	cuts.SetEtaRange(fEta[0], fEta[1]);
	227
	228	// transverse DCA cuts
	229	if (fDCARfixed)
	230	cuts.SetMaxDCAToVertexXY(fDCARmax);
	231	else
	232	cuts.SetMaxDCAToVertexXYPtDep(fDCARptFormula.Data());
	233
	234	// longitudinal DCA cuts
	235	if (fDCAZfixed)
	236	cuts.SetMaxDCAToVertexZ(fDCAZmax);
	237	else
	238	cuts.SetMaxDCAToVertexZPtDep(fDCAZptFormula.Data());
	239
	240	// these options are always disabled in current version
	241	cuts.SetDCAToVertex2D(kFALSE);
	242	cuts.SetRequireSigmaToVertex(kFALSE);
	243
	244	// TPC related cuts for TPC+ITS tracks
	245	cuts.SetMinNClustersTPC(fTPCminNClusters);
	246	cuts.SetMaxChi2PerClusterTPC(fTPCmaxChi2);
	247	cuts.SetAcceptKinkDaughters(!fRejectKinkDaughters);
	248	cuts.SetMaxChi2TPCConstrainedGlobal(fCutMaxChi2TPCConstrainedVsGlobal);
	249
	250	// ITS related cuts for TPC+ITS tracks
	251	if (fSPDminNClusters > 0)
	252	cuts.SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
	253	cuts.SetMaxChi2PerClusterITS(fITSmaxChi2);
	254
	255	// now that all is initialized, do the check
	256	return cuts.IsSelected(track);
	257	}
	258
	259	//_________________________________________________________________________________________________
	260	Bool_t AliRsnCutTrackQuality::CheckAOD(AliAODTrack *track)
	261	{
	262	//
	263	// Check an AOD track.
	264	// This is done doing directly all checks, since there is not
	265	// an equivalend checker for AOD tracks
	266	//
	267
	268	// if a test bit is used, check it and skip the following
	269	if (fAODTestFilterBit >= 0) {
	270	UInt_t bit = 1 << fAODTestFilterBit;
	271	AliDebugClass(2, Form("Required a test filter bit for AOD check: %u (result: %s)", bit, (track->TestFilterBit(bit) ? "accept" : "reject")));
	272	if (!track->TestFilterBit(bit))
	273	return kFALSE;
	274	else {
	275	if (track->Pt() < fPt[0] \|\| track->Pt() > fPt[1]) return kFALSE;
	276	if (track->Eta() < fEta[0] \|\| track->Eta() > fEta[1]) return kFALSE;
	277	return kTRUE;
	278	}
	279	}
	280
	281	// try to retrieve the reference AOD event
	282	AliAODEvent *aodEvent = 0x0;
	283	if (fEvent) aodEvent = fEvent->GetRefAOD();
	284	if (!aodEvent) {
	285	AliError("AOD reference event is not initialized!");
	286	return kFALSE;
	287	}
	288
	289	// step #0: check SPD and ITS clusters
	290	Int_t nSPD = 0;
	291	nSPD = TESTBIT(track->GetITSClusterMap(), 0);
	292	nSPD += TESTBIT(track->GetITSClusterMap(), 1);
	293	if (nSPD < fSPDminNClusters) {
	294	AliDebug(AliLog::kDebug + 2, "Not enough SPD clusters in this track. Rejected");
	295	return kFALSE;
	296	}
	297
	298
	299	// step #1: check number of clusters in TPC
	300	if (track->GetTPCNcls() < fTPCminNClusters) {
	301	AliDebug(AliLog::kDebug + 2, "Too few TPC clusters. Rejected");
	302	return kFALSE;
	303	}
	304	if (track->GetITSNcls() < fITSminNClusters) {
	305	AliDebug(AliLog::kDebug + 2, "Too few ITS clusters. Rejected");
	306	return kFALSE;
	307	}
	308
	309	// step #2: check chi square
	310	if (track->Chi2perNDF() > fTPCmaxChi2) {
	311	AliDebug(AliLog::kDebug + 2, "Bad chi2. Rejected");
	312	return kFALSE;
	313	}
	314	if (track->Chi2perNDF() > fITSmaxChi2) {
	315	AliDebug(AliLog::kDebug + 2, "Bad chi2. Rejected");
	316	return kFALSE;
	317	}
	318
	319	// step #3: reject kink daughters
	320	AliAODVertex *vertex = track->GetProdVertex();
	321	if (vertex && fRejectKinkDaughters) {
	322	if (vertex->GetType() == AliAODVertex::kKink) {
	323	AliDebug(AliLog::kDebug + 2, "Kink daughter. Rejected");
	324	return kFALSE;
	325	}
	326	}
	327
	328	// step #4: DCA cut (transverse)
	329	// --> reject all tracks not ITS refitted
	330	Double_t b[2], cov[3];
	331	vertex = aodEvent->GetPrimaryVertex();
	332	if (!vertex) {
	333	AliDebug(AliLog::kDebug + 2, "NULL vertex");
	334	return kFALSE;
	335	}
	336	if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
	337	AliDebug(AliLog::kDebug + 2, "Not ITS refitted");
	338	return kFALSE;
	339	}
	340	if (!track->PropagateToDCA(vertex, aodEvent->GetMagneticField(), kVeryBig, b, cov)) {
	341	AliDebug(AliLog::kDebug + 2, "Failed propagation to vertex");
	342	return kFALSE;
	343	}
	344	// if the DCA cut is not fixed, compute current value
	345	if (!fDCARfixed) {
	346	static TString str(fDCARptFormula);
	347	str.ReplaceAll("pt", "x");
	348	static const TFormula dcaXY(Form("%s_dcaXY", GetName()), str.Data());
	349	fDCARmax = dcaXY.Eval(track->Pt());
	350	}
	351	// check the cut
	352	if (TMath::Abs(b[0]) > fDCARmax) {
	353	AliDebug(AliLog::kDebug + 2, "Too large transverse DCA");
	354	return kFALSE;
	355	}
	356
	357	// step #5: DCA cut (longitudinal)
	358	// the DCA has already been computed above
	359	// if the DCA cut is not fixed, compute current value
	360	if (!fDCAZfixed) {
	361	static TString str(fDCAZptFormula);
	362	str.ReplaceAll("pt", "x");
	363	static const TFormula dcaZ(Form("%s_dcaXY", GetName()), str.Data());
	364	fDCAZmax = dcaZ.Eval(track->Pt());
	365	}
	366	// check the cut
	367	if (TMath::Abs(b[1]) > fDCAZmax) {
	368	AliDebug(AliLog::kDebug + 2, "Too large longitudinal DCA");
	369	return kFALSE;
	370	}
	371
	372	// step #6: check eta/pt range
	373	if (track->Eta() < fEta[0] \|\| track->Eta() > fEta[1]) {
	374	AliDebug(AliLog::kDebug + 2, "Outside ETA acceptance");
	375	return kFALSE;
	376	}
	377	if (track->Pt() < fPt[0] \|\| track->Pt() > fPt[1]) {
	378	AliDebug(AliLog::kDebug + 2, "Outside PT acceptance");
	379	return kFALSE;
	380	}
	381
	382	// if we are here, all cuts were passed and no exit point was got
	383	AliDebug(AliLog::kDebug + 2, "============================= ACCEPTED TRACK =====================================================");
	384	return kTRUE;
	385	}
	386
	387	//_________________________________________________________________________________________________
	388	void AliRsnCutTrackQuality::Print(const Option_t *) const
	389	{
	390	//
	391	// Print information on this cut
	392	//
	393
	394	AliInfo(Form("Cut name : %s", GetName()));
	395	AliInfo(Form("Required flags (off, on): %lx %lx", fFlagsOn, fFlagsOff));
	396	AliInfo(Form("Ranges in eta, pt : %.2f - %.2f, %.2f - %.2f", fEta[0], fEta[1], fPt[0], fPt[1]));
	397	AliInfo(Form("Kink daughters are : %s", (fRejectKinkDaughters ? "rejected" : "accepted")));
	398	AliInfo(Form("TPC requirements : min. cluster = %d, max chi2 = %f", fTPCminNClusters, fTPCmaxChi2));
	399	AliInfo(Form("ITS requirements : min. cluster = %d (all), %d (SPD), max chi2 = %f", fITSminNClusters, fSPDminNClusters, fITSmaxChi2));
	400
	401	if (fDCARfixed) {
	402	AliInfo(Form("DCA r cut : fixed to %f cm", fDCARmax));
	403	} else {
	404	AliInfo(Form("DCA r cut formula : %s", fDCARptFormula.Data()));
	405	}
	406
	407	if (fDCAZfixed) {
	408	AliInfo(Form("DCA z cut : fixed to %f cm", fDCAZmax));
	409	} else {
	410	AliInfo(Form("DCA z cut formula : %s", fDCAZptFormula.Data()));
	411	}
	412	}
	413	//__________________________________________________________________________________________________
	414	void AliRsnCutTrackQuality::SetDefaults2010()
	415	{
	416	//
	417	// Default settings for cuts used in 2010
	418	//
	419	AddStatusFlag(AliESDtrack::kTPCin , kTRUE);
	420	AddStatusFlag(AliESDtrack::kTPCrefit, kTRUE);
	421	AddStatusFlag(AliESDtrack::kITSrefit, kTRUE);
	422	SetPtRange(0.15, 1E+20);
	423	SetEtaRange(-0.8, 0.8);
	424	SetDCARPtFormula("0.0182+0.0350/pt^1.01");
	425	SetDCAZmax(2.0);
	426	SetSPDminNClusters(1);
	427	SetITSminNClusters(0);
	428	SetITSmaxChi2(36);
	429	SetMaxChi2TPCConstrainedGlobal(36);
	430	SetTPCminNClusters(70);
	431	SetTPCmaxChi2(4.0);
	432	SetRejectKinkDaughters();
	433	SetAODTestFilterBit(5);
	434	}
	435
	436	//__________________________________________________________________________________________________
	437	const char *AliRsnCutTrackQuality::Binary(UInt_t number)
	438	{
	439	//
	440	// Convert an integer in binary
	441	//
	442
	443	static char b[50];
	444	b[0] = '\0';
	445
	446	UInt_t z;
	447	for (z = 512; z > 0; z >>= 1)
	448	strncat(b, ((number & z) == z) ? "1" : "0", 1);
	449
	450	return b;
	451	}