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

//
// Class AliRsnCutV0
//
// 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: Massimo Venaruzzo (massimo.venaruzzo@ts.infn.it)
// modified: Enrico Fragiacomo (enrico.fragiacomo@ts.infn.it)
//

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

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

#include "AliRsnEvent.h"
#include "AliRsnDaughter.h"
#include "AliRsnCutV0.h"

ClassImp(AliRsnCutV0)

//_________________________________________________________________________________________________
//AliRsnCutV0::AliRsnCutV0(const char *name, Int_t hypothesis) :
AliRsnCutV0::AliRsnCutV0(const char *name, Int_t hypothesis, AliPID::EParticleType pid, AliPID::EParticleType pid2) :
   AliRsnCut(name, AliRsnTarget::kDaughter),
   fHypothesis(0),
   fMass(0.0),
   fTolerance(0.01),
   fMaxDCAVertex(0.3),
   fMinCosPointAngle(0.95),
   fMaxDaughtersDCA(0.5),
   fMinTPCcluster(70),
   fMaxRapidity(0.8),
   fPID(pid),
   fPID2(pid2),
   fPIDCut1(0),
   fPIDCut2(0),
   fPIDCut3(0),
   fESDtrackCuts(0x0),
   fCutQuality(Form("%sDaughtersQuality", name)),
   fAODTestFilterBit(5)
{
//
// Default constructor.
// Initializes all cuts in such a way that all of them are disabled.
//

   SetHypothesis(hypothesis);
}

//_________________________________________________________________________________________________
AliRsnCutV0::AliRsnCutV0(const AliRsnCutV0 &copy) :
   AliRsnCut(copy),
   fHypothesis(copy.fHypothesis),
   fMass(copy.fMass),
   fTolerance(copy.fTolerance),
   fMaxDCAVertex(copy.fMaxDCAVertex),
   fMinCosPointAngle(copy.fMinCosPointAngle),
   fMaxDaughtersDCA(copy.fMaxDaughtersDCA),
   fMinTPCcluster(copy.fMinTPCcluster),
   fMaxRapidity(copy.fMaxRapidity),
   fPID(copy.fPID),
   fPID2(copy.fPID2),
   fPIDCut1(copy.fPIDCut1),
   fPIDCut2(copy.fPIDCut2),
   fPIDCut3(copy.fPIDCut3),
   fESDtrackCuts(copy.fESDtrackCuts),
   fCutQuality(copy.fCutQuality),
   fAODTestFilterBit(copy.fAODTestFilterBit)
{
//
// Copy constructor.
// Just copy all data member values.:IsSelected: Object is not a V0 (RESONANCES/AliRsnCutV0.cxx:149)

//
   fCutQuality.SetPtRange(0.15, 1E+20);
   fCutQuality.SetEtaRange(-0.8, 0.8);
   fCutQuality.SetDCARmax(0.05);
   //fCutQuality.SetDCARPtFormula("0.0182+0.0350/pt^1.01");
   fCutQuality.SetDCAZmax(2.0);
   fCutQuality.SetSPDminNClusters(1);
   fCutQuality.SetITSminNClusters(0);
   fCutQuality.SetITSmaxChi2(1E+20);
   fCutQuality.SetTPCminNClusters(fMinTPCcluster);
   fCutQuality.SetTPCmaxChi2(4.0);
   fCutQuality.SetRejectKinkDaughters();
   fCutQuality.SetAODTestFilterBit(5);

}

//_________________________________________________________________________________________________
AliRsnCutV0 &AliRsnCutV0::operator=(const AliRsnCutV0 &copy)
{
//
// Assignment operator.
// Just copy all data member values.
//
   if (this == &copy)
     return *this;
   fHypothesis = copy.fHypothesis;
   fMass = copy.fMass;
   fTolerance = copy.fTolerance;
   fMaxDCAVertex = copy.fMaxDCAVertex;
   fMinCosPointAngle = copy.fMinCosPointAngle;
   fMaxDaughtersDCA = copy.fMaxDaughtersDCA;
   fMinTPCcluster = copy.fMinTPCcluster;
   fMaxRapidity = copy.fMaxRapidity;
   fPID = copy.fPID;
   fPID2 = copy.fPID2;
   fPIDCut1 = copy.fPIDCut1;
   fPIDCut2 = copy.fPIDCut2;
   fPIDCut3 = copy.fPIDCut3;
   fESDtrackCuts = copy.fESDtrackCuts;
   fCutQuality = copy.fCutQuality;
   fAODTestFilterBit = copy.fAODTestFilterBit;

   return (*this);
}

//_________________________________________________________________________________________________
Bool_t AliRsnCutV0::IsSelected(TObject *object)
{
//:IsSelected: Object is not a V0 (RESONANCES/AliRsnCutV0.cxx:149)

// 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;

   // check cast
   AliESDv0 *v0esd = fDaughter->Ref2ESDv0();
   AliAODv0 *v0aod = fDaughter->Ref2AODv0();
   //cout << fDaughter->GetRef()->ClassName() << ' ' << v0esd << ' ' << v0aod << endl;

   // operate depending on cast:IsSelected: Object is not a V0 (RESONANCES/AliRsnCutV0.cxx:149)

   if (v0esd) {
      return CheckESD(v0esd);
   } else if (v0aod) {
      return CheckAOD(v0aod);
   } else {
      AliDebugClass(1, "Object is not a V0");
      return kFALSE;
   }
}

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

   AliDebugClass(1, "Check ESD");
   if (v0->GetOnFlyStatus()) {
      AliDebugClass(1, "Rejecting V0 in 'on fly' status");
      return kFALSE; // if kTRUE, then this V0 is recontructed
   }

   // retrieve pointer to owner event
   AliESDEvent *lESDEvent = fEvent->GetRefESD();
   Double_t xPrimaryVertex = lESDEvent->GetPrimaryVertex()->GetX();
   Double_t yPrimaryVertex = lESDEvent->GetPrimaryVertex()->GetY();
   Double_t zPrimaryVertex = lESDEvent->GetPrimaryVertex()->GetZ();
   AliDebugClass(2, Form("Primary vertex: %f %f %f", xPrimaryVertex, yPrimaryVertex, zPrimaryVertex));

   // retrieve the V0 daughters
   UInt_t lIdxPos      = (UInt_t) TMath::Abs(v0->GetPindex());
   UInt_t lIdxNeg      = (UInt_t) TMath::Abs(v0->GetNindex());
   AliESDtrack *pTrack = lESDEvent->GetTrack(lIdxPos);
   AliESDtrack *nTrack = lESDEvent->GetTrack(lIdxNeg);

   // check quality cuts
   if (fESDtrackCuts) {
      AliDebugClass(2, "Checking quality cuts");
      if (!fESDtrackCuts->IsSelected(pTrack)) {
         AliDebugClass(2, "Positive daughter failed quality cuts");
         return kFALSE;
      }
      if (!fESDtrackCuts->IsSelected(nTrack)) {
         AliDebugClass(2, "Negative daughter failed quality cuts");
         return kFALSE;
      }
   }

   // filter like-sign V0
   if ( TMath::Abs( ((pTrack->GetSign()) - (nTrack->GetSign())) ) < 0.1) {
      AliDebugClass(2, "Failed like-sign V0 check");
      return kFALSE;
   }


   // check compatibility with expected species hypothesis
   v0->ChangeMassHypothesis(fHypothesis);
   if ((TMath::Abs(v0->GetEffMass() - fMass)) > fTolerance) {
      AliDebugClass(2, "V0 is not in the expected inv mass range");
      return kFALSE;
   }

   // topological checks
   if (TMath::Abs(v0->GetD(xPrimaryVertex, yPrimaryVertex, zPrimaryVertex)) > fMaxDCAVertex) {
      AliDebugClass(2, "Failed check on DCA to primary vertes");
      return kFALSE;
   }
   if (TMath::Abs(v0->GetV0CosineOfPointingAngle()) < fMinCosPointAngle) {
      AliDebugClass(2, "Failed check on cosine of pointing angle");
      return kFALSE;
   }
   if (TMath::Abs(v0->GetDcaV0Daughters()) > fMaxDaughtersDCA) {
      AliDebugClass(2, "Failed check on DCA between daughters");
      return kFALSE;
   }
   if (TMath::Abs(v0->Y(fHypothesis)) > fMaxRapidity) {
      AliDebugClass(2, "Failed check on V0 rapidity");
      return kFALSE;
   }


   // check PID on proton or antiproton from V0

   // check initialization of PID object
   AliPIDResponse *pid = fEvent->GetPIDResponse();
   if (!pid) {
      AliFatal("NULL PID response");
      return kFALSE;
   }

   // check if TOF is matched
   // and computes all values used in the PID cut
   //Bool_t   isTOFpos  = MatchTOF(ptrack);
   //Bool_t   isTOFneg  = MatchTOF(ntrack);
   Double_t pospTPC   = pTrack->GetTPCmomentum();
   Double_t negpTPC   = nTrack->GetTPCmomentum();
   //Double_t posp      = pTrack->P();
   //Double_t negp      = nTrack->P();
   Double_t posnsTPC   = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID));
   Double_t posnsTPC2  = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID2));
   //Double_t posnsTOF  = TMath::Abs(pid->NumberOfSigmasTOF(ptrack, fPID));
   Double_t negnsTPC   = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID));
   Double_t negnsTPC2  = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID2));
   //Double_t negnsTOF  = TMath::Abs(pid->NumberOfSigmasTOF(ntrack, fPID));
   Double_t maxTPC = 1E20;
   Double_t maxTPC2 = 1E20;
   //Double_t maxTOF = 1E20;

   // applies the cut differently depending on the PID and the momentum

   if(fHypothesis==kLambda0) {
      //if (isTOFpos) {
      // TPC: 5sigma cut for all
      //if (posnsTPC > 5.0) return kFALSE;
      // TOF: 3sigma
      // maxTOF = 3.0;
      //return (posnsTOF <= maxTOF);
      //} else {
      // TPC:
      // below 600 MeV: 4sigma
      // above 600 MeV: 3sigma

      if (pospTPC <= 0.6 && fPID==AliPID::kProton)
         maxTPC = fPIDCut1;
      else if (pospTPC > 0.6 && fPID==AliPID::kProton)
         maxTPC = fPIDCut2;
      else
         return kFALSE;

      maxTPC2 = fPIDCut3;

      if (! ((posnsTPC <= maxTPC) && (negnsTPC2 <= maxTPC2)) ) {
         AliDebugClass(2, "Failed check on V0 PID");
         return kFALSE;
      }
   }

   //}

   if(fHypothesis==kLambda0Bar) {
      //if (isTOFneg) {
      // TPC: 5sigma cut for all
      //if (negnsTPC > 5.0) return kFALSE;
      // TOF: 3sigma
      // maxTOF = 3.0;
      //return (negnsTOF <= maxTOF);
      //} else {
      // TPC:
      // below 600 MeV: 4sigma
      // above 600 MeV: 3sigma

      if (negpTPC <= 0.6 && fPID==AliPID::kProton)
         maxTPC = fPIDCut1;
      else if (negpTPC > 0.6 && fPID==AliPID::kProton)
         maxTPC = fPIDCut2;
      else
         return kFALSE;

      maxTPC2 = fPIDCut3;

      if(! ((negnsTPC <= maxTPC) && (posnsTPC2 <= maxTPC2)) ) {
         AliDebugClass(2, "Failed check on V0 PID");
         return kFALSE;
      }
   }
   //}


   // if we reach this point, all checks were successful
   AliDebugClass(2, "Good V0 (hallelujah)");
   return kTRUE;
}

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

   AliDebugClass(2, "Check AOD");
   if (v0->GetOnFlyStatus()) {
      AliDebugClass(2, "Rejecting V0 in 'on fly' status");
      return kFALSE; // if kTRUE, then this V0 is recontructed
   }


   // retrieve pointer to owner event
   AliAODEvent *lAODEvent = fEvent->GetRefAOD();
   Double_t xPrimaryVertex = lAODEvent->GetPrimaryVertex()->GetX();
   Double_t yPrimaryVertex = lAODEvent->GetPrimaryVertex()->GetY();
   Double_t zPrimaryVertex = lAODEvent->GetPrimaryVertex()->GetZ();
   AliDebugClass(2, Form("Primary vertex: %f %f %f", xPrimaryVertex, yPrimaryVertex, zPrimaryVertex));

   // retrieve the V0 daughters
   AliAODTrack *pTrack = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
   AliAODTrack *nTrack = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));


   // check quality cuts
   UInt_t  filtermapP = 9999;
   UInt_t  filtermapN = 9999;
   filtermapP = pTrack->GetFilterMap();
   filtermapN = nTrack->GetFilterMap();


   if ( !pTrack->TestFilterBit(fAODTestFilterBit)   ) {
      AliDebugClass(2, Form("Positive daughter failed quality cuts filtermapP=%d",filtermapP));
      return kFALSE;
   }
   if ( !nTrack->TestFilterBit(fAODTestFilterBit)   ) {
      AliDebugClass(2, Form("Negative daughter failed quality cuts filtermapN=%d",filtermapN));
      return kFALSE;
   }


   /*AliDebugClass(1, Form("fESDtrackCuts=%p",fESDtrackCuts));
   if (fESDtrackCuts) { // use fESDtrackCuts to retrieve cuts values


      AliDebugClass(2, "Checking quality cuts");


      const Bool_t  CutAcceptKinkDaughters  = fESDtrackCuts->GetAcceptKinkDaughters(); // 0 = kFalse
      const Float_t CutMaxChi2PerClusterTPC = fESDtrackCuts->GetMaxChi2PerClusterTPC();
      const Int_t   CutMinNClusterTPC       = fESDtrackCuts->GetMinNClusterTPC();
      const Bool_t  CutRequireTPCRefit      = fESDtrackCuts->GetRequireTPCRefit();
      //AliDebugClass(2, Form("accept kink=%d maxchi2=%f minnclSTPC=%d requireTPCrefit=%d", CutAcceptKinkDaughters,
      //        CutMaxChi2PerClusterTPC, CutMinNClusterTPC, CutRequireTPCRefit));
      //AliDebugClass(2, Form("pTrack->TPCNcls=%d", pTrack->GetTPCNcls() ));
      //AliDebugClass(2, Form("nTrack->TPCNcls=%d", nTrack->GetTPCNcls() ));


      if(pTrack->GetTPCNcls() < CutMinNClusterTPC) {AliDebugClass(2, "Positive daughter not MinNclsTPC"); return kFALSE;}
      if(nTrack->GetTPCNcls() < CutMinNClusterTPC) {AliDebugClass(2, "Negative daughter not MinNclsTPC"); return kFALSE;}

   }*/

   // filter like-sign V0
   if ( TMath::Abs( ((pTrack->Charge()) - (nTrack->Charge())) ) < 0.1) {
      AliDebugClass(2, "Failed like-sign V0 check");
      return kFALSE;
   }

   // check compatibility with expected species hypothesis
   Double_t mass = 0.0;
   if(fHypothesis==kLambda0) {
      mass = v0->MassLambda();
   }
   else if (fHypothesis==kLambda0Bar) {
      mass = v0->MassAntiLambda();
   }
   if ((TMath::Abs(mass - fMass)) > fTolerance) {
      AliDebugClass(2, Form("V0 is not in the expected inv mass range  Mass: %d %f %f", fHypothesis, fMass, mass));
      return kFALSE;
   }
   AliDebugClass(2, Form("Mass: %d %f %f", fHypothesis, fMass, mass));


   // topological checks
   if (TMath::Abs(v0->DcaV0ToPrimVertex()) > fMaxDCAVertex) {
      AliDebugClass(2, Form("Failed check on DCA to primary vertes dca=%f maxdca=%f",TMath::Abs(v0->DcaV0ToPrimVertex()),fMaxDCAVertex));
      return kFALSE;
   }
   if (TMath::Abs( TMath::Cos(v0->OpenAngleV0()) ) < fMinCosPointAngle) {
      AliDebugClass(2, "Failed check on cosine of pointing angle");
      return kFALSE;
   }
   if (TMath::Abs(v0->DcaV0Daughters()) > fMaxDaughtersDCA) {
      AliDebugClass(2, "Failed check on DCA between daughters");
      return kFALSE;
   }
   if (TMath::Abs(v0->RapLambda()) > fMaxRapidity) {
      AliDebugClass(2, "Failed check on V0 rapidity");
      return kFALSE;
   }


   // check initialization of PID object
   AliPIDResponse *pid = fEvent->GetPIDResponse();
   if (!pid) {
      AliFatal("NULL PID response");
      return kFALSE;
   }

   // check if TOF is matched
   // and computes all values used in the PID cut
   //Bool_t   isTOFpos  = MatchTOF(ptrack);
   //Bool_t   isTOFneg  = MatchTOF(ntrack);
   Double_t pospTPC   = pTrack->GetTPCmomentum();
   Double_t negpTPC   = nTrack->GetTPCmomentum();
   //Double_t posp      = pTrack->P();
   //Double_t negp      = nTrack->P();
   Double_t posnsTPC   = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID));
   Double_t posnsTPC2  = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID2));
   //Double_t posnsTOF  = TMath::Abs(pid->NumberOfSigmasTOF(ptrack, fPID));
   Double_t negnsTPC   = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID));
   Double_t negnsTPC2  = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID2));
   //Double_t negnsTOF  = TMath::Abs(pid->NumberOfSigmasTOF(ntrack, fPID));
   Double_t maxTPC = 1E20;
   Double_t maxTPC2 = 1E20;
   //Double_t maxTOF = 1E20;

   // applies the cut differently depending on the PID and the momentum

   if(fHypothesis==kLambda0) {


      //if (isTOFpos) {
      // TPC: 5sigma cut for all
      //if (posnsTPC > 5.0) return kFALSE;
      // TOF: 3sigma
      // maxTOF = 3.0;
      //return (posnsTOF <= maxTOF);
      //} else {
      // TPC:
      // below 600 MeV: 4sigma
      // above 600 MeV: 3sigma

      if (pospTPC <= 0.6 && fPID==AliPID::kProton)
         maxTPC = fPIDCut1; // EF safer value to run on MC
      else if (pospTPC > 0.6 && fPID==AliPID::kProton)
         maxTPC = fPIDCut2; // EF safer value to run on MC
      else
         return kFALSE;

      maxTPC2 = fPIDCut3; // EF safer value to run on MC

      if (! ((posnsTPC <= maxTPC) && (negnsTPC2 <= maxTPC2)) ) {
         AliDebugClass(2, "Failed check on V0 PID");
         return kFALSE;
      }
   }

   if(fHypothesis==kLambda0Bar) {

      //if (isTOFneg) {
      // TPC: 5sigma cut for all
      //if (negnsTPC > 5.0) return kFALSE;
      // TOF: 3sigma
      // maxTOF = 3.0;
      //return (negnsTOF <= maxTOF);
      //} else {
      // TPC:
      // below 600 MeV: 4sigma
      // above 600 MeV: 3sigma

      if (negpTPC <= 0.6 && fPID==AliPID::kProton)
         maxTPC = fPIDCut1; // EF safer value to run on MC
      else if (negpTPC > 0.6 && fPID==AliPID::kProton)
         maxTPC = fPIDCut2; // EF safer value to run on MC
      else
         return kFALSE;

      maxTPC2 = fPIDCut3; // EF safer value to run on MC

      if(! ((negnsTPC <= maxTPC) && (posnsTPC2 <= maxTPC2)) ) {
         AliDebugClass(2, "Failed check on V0 PID");
         return kFALSE;
      }
   }


   // if we reach this point, all checks were successful
   AliDebugClass(1, "Good AOD V0 (hallelujah)");
   AliDebugClass(1, Form("Mass: %d %f %f %d %d", fHypothesis, fMass, mass, filtermapP, filtermapN));
   return kTRUE;

}

//_________________________________________________________________________________________________
void AliRsnCutV0::Print(const Option_t *) const
{
//
// Print information on this cut
//
}
Commit	Line	Data
8082e538	1	//
	2	// Class AliRsnCutV0
	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	//
76d89202	15	// authors: Massimo Venaruzzo (massimo.venaruzzo@ts.infn.it)
76d89202	16	// modified: Enrico Fragiacomo (enrico.fragiacomo@ts.infn.it)
8082e538	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 "AliRsnCutV0.h"
	29
	30	ClassImp(AliRsnCutV0)
	31
	32	//_________________________________________________________________________________________________
3da8cef7	33	//AliRsnCutV0::AliRsnCutV0(const char *name, Int_t hypothesis) :
3da8cef7	34	AliRsnCutV0::AliRsnCutV0(const char *name, Int_t hypothesis, AliPID::EParticleType pid, AliPID::EParticleType pid2) :
8082e538	35	AliRsnCut(name, AliRsnTarget::kDaughter),
	36	fHypothesis(0),
	37	fMass(0.0),
3da8cef7	38	fTolerance(0.01),
8082e538	39	fMaxDCAVertex(0.3),
8082e538	40	fMinCosPointAngle(0.95),
3da8cef7	41	fMaxDaughtersDCA(0.5),
9d4bb6d8	42	fMinTPCcluster(70),
3da8cef7	43	fMaxRapidity(0.8),
	44	fPID(pid),
	45	fPID2(pid2),
76d89202	46	fPIDCut1(0),
	47	fPIDCut2(0),
	48	fPIDCut3(0),
	49	fESDtrackCuts(0x0),
	50	fCutQuality(Form("%sDaughtersQuality", name)),
	51	fAODTestFilterBit(5)
8082e538	52	{
	53	//
	54	// Default constructor.
	55	// Initializes all cuts in such a way that all of them are disabled.
	56	//
	57
	58	SetHypothesis(hypothesis);
	59	}
	60
	61	//_________________________________________________________________________________________________
	62	AliRsnCutV0::AliRsnCutV0(const AliRsnCutV0 &copy) :
	63	AliRsnCut(copy),
	64	fHypothesis(copy.fHypothesis),
	65	fMass(copy.fMass),
	66	fTolerance(copy.fTolerance),
	67	fMaxDCAVertex(copy.fMaxDCAVertex),
	68	fMinCosPointAngle(copy.fMinCosPointAngle),
	69	fMaxDaughtersDCA(copy.fMaxDaughtersDCA),
9d4bb6d8	70	fMinTPCcluster(copy.fMinTPCcluster),
e6952ec7	71	fMaxRapidity(copy.fMaxRapidity),
3da8cef7	72	fPID(copy.fPID),
	73	fPID2(copy.fPID2),
	74	fPIDCut1(copy.fPIDCut1),
	75	fPIDCut2(copy.fPIDCut2),
	76	fPIDCut3(copy.fPIDCut3),
76d89202	77	fESDtrackCuts(copy.fESDtrackCuts),
	78	fCutQuality(copy.fCutQuality),
	79	fAODTestFilterBit(copy.fAODTestFilterBit)
8082e538	80	{
	81	//
	82	// Copy constructor.
76d89202	83	// Just copy all data member values.:IsSelected: Object is not a V0 (RESONANCES/AliRsnCutV0.cxx:149)
76d89202	84
8082e538	85	//
76d89202	86	fCutQuality.SetPtRange(0.15, 1E+20);
	87	fCutQuality.SetEtaRange(-0.8, 0.8);
	88	fCutQuality.SetDCARmax(0.05);
	89	//fCutQuality.SetDCARPtFormula("0.0182+0.0350/pt^1.01");
	90	fCutQuality.SetDCAZmax(2.0);
	91	fCutQuality.SetSPDminNClusters(1);
	92	fCutQuality.SetITSminNClusters(0);
	93	fCutQuality.SetITSmaxChi2(1E+20);
9d4bb6d8	94	fCutQuality.SetTPCminNClusters(fMinTPCcluster);
76d89202	95	fCutQuality.SetTPCmaxChi2(4.0);
	96	fCutQuality.SetRejectKinkDaughters();
	97	fCutQuality.SetAODTestFilterBit(5);
	98
8082e538	99	}
	100
	101	//_________________________________________________________________________________________________
e6952ec7	102	AliRsnCutV0 &AliRsnCutV0::operator=(const AliRsnCutV0 &copy)
8082e538	103	{
	104	//
	105	// Assignment operator.
	106	// Just copy all data member values.
	107	//
4cc7af7d	108	if (this == &copy)
4cc7af7d	109	return *this;
61f275d1	110	fHypothesis = copy.fHypothesis;
	111	fMass = copy.fMass;
	112	fTolerance = copy.fTolerance;
	113	fMaxDCAVertex = copy.fMaxDCAVertex;
	114	fMinCosPointAngle = copy.fMinCosPointAngle;
	115	fMaxDaughtersDCA = copy.fMaxDaughtersDCA;
9d4bb6d8	116	fMinTPCcluster = copy.fMinTPCcluster;
3da8cef7	117	fMaxRapidity = copy.fMaxRapidity;
	118	fPID = copy.fPID;
	119	fPID2 = copy.fPID2;
	120	fPIDCut1 = copy.fPIDCut1;
	121	fPIDCut2 = copy.fPIDCut2;
	122	fPIDCut3 = copy.fPIDCut3;
61f275d1	123	fESDtrackCuts = copy.fESDtrackCuts;
76d89202	124	fCutQuality = copy.fCutQuality;
76d89202	125	fAODTestFilterBit = copy.fAODTestFilterBit;
61f275d1	126
8082e538	127	return (*this);
e6952ec7	128	}
8082e538	129
	130	//_________________________________________________________________________________________________
	131	Bool_t AliRsnCutV0::IsSelected(TObject *object)
	132	{
76d89202	133	//:IsSelected: Object is not a V0 (RESONANCES/AliRsnCutV0.cxx:149)
76d89202	134
8082e538	135	// Cut checker.
	136	// Checks the type of object being evaluated
	137	// and then calls the appropriate sub-function (for ESD or AOD)
	138	//
	139
	140	// coherence check
	141	if (!TargetOK(object)) return kFALSE;
e6952ec7	142
8082e538	143	// check cast
	144	AliESDv0 *v0esd = fDaughter->Ref2ESDv0();
	145	AliAODv0 *v0aod = fDaughter->Ref2AODv0();
	146	//cout << fDaughter->GetRef()->ClassName() << ' ' << v0esd << ' ' << v0aod << endl;
e6952ec7	147
76d89202	148	// operate depending on cast:IsSelected: Object is not a V0 (RESONANCES/AliRsnCutV0.cxx:149)
61f275d1	149
8082e538	150	if (v0esd) {
	151	return CheckESD(v0esd);
	152	} else if (v0aod) {
	153	return CheckAOD(v0aod);
	154	} else {
	155	AliDebugClass(1, "Object is not a V0");
	156	return kFALSE;
	157	}
	158	}
	159
	160	//_________________________________________________________________________________________________
	161	Bool_t AliRsnCutV0::CheckESD(AliESDv0 *v0)
	162	{
	163	//
	164	// Check an ESD V0.
	165	// This is done using the default track checker for ESD.
	166	// It is declared static, not to recreate it every time.
	167	//
	168
	169	AliDebugClass(1, "Check ESD");
	170	if (v0->GetOnFlyStatus()) {
	171	AliDebugClass(1, "Rejecting V0 in 'on fly' status");
	172	return kFALSE; // if kTRUE, then this V0 is recontructed
	173	}
e6952ec7	174
8082e538	175	// retrieve pointer to owner event
	176	AliESDEvent *lESDEvent = fEvent->GetRefESD();
	177	Double_t xPrimaryVertex = lESDEvent->GetPrimaryVertex()->GetX();
	178	Double_t yPrimaryVertex = lESDEvent->GetPrimaryVertex()->GetY();
	179	Double_t zPrimaryVertex = lESDEvent->GetPrimaryVertex()->GetZ();
	180	AliDebugClass(2, Form("Primary vertex: %f %f %f", xPrimaryVertex, yPrimaryVertex, zPrimaryVertex));
e6952ec7	181
8082e538	182	// retrieve the V0 daughters
	183	UInt_t lIdxPos = (UInt_t) TMath::Abs(v0->GetPindex());
	184	UInt_t lIdxNeg = (UInt_t) TMath::Abs(v0->GetNindex());
	185	AliESDtrack *pTrack = lESDEvent->GetTrack(lIdxPos);
	186	AliESDtrack *nTrack = lESDEvent->GetTrack(lIdxNeg);
e6952ec7	187
8082e538	188	// check quality cuts
	189	if (fESDtrackCuts) {
	190	AliDebugClass(2, "Checking quality cuts");
	191	if (!fESDtrackCuts->IsSelected(pTrack)) {
	192	AliDebugClass(2, "Positive daughter failed quality cuts");
	193	return kFALSE;
	194	}
	195	if (!fESDtrackCuts->IsSelected(nTrack)) {
	196	AliDebugClass(2, "Negative daughter failed quality cuts");
	197	return kFALSE;
	198	}
	199	}
e6952ec7	200
8082e538	201	// filter like-sign V0
e6952ec7	202	if ( TMath::Abs( ((pTrack->GetSign()) - (nTrack->GetSign())) ) < 0.1) {
	203	AliDebugClass(2, "Failed like-sign V0 check");
	204	return kFALSE;
	205	}
	206
61f275d1	207
8082e538	208	// check compatibility with expected species hypothesis
	209	v0->ChangeMassHypothesis(fHypothesis);
	210	if ((TMath::Abs(v0->GetEffMass() - fMass)) > fTolerance) {
	211	AliDebugClass(2, "V0 is not in the expected inv mass range");
	212	return kFALSE;
	213	}
e6952ec7	214
8082e538	215	// topological checks
	216	if (TMath::Abs(v0->GetD(xPrimaryVertex, yPrimaryVertex, zPrimaryVertex)) > fMaxDCAVertex) {
	217	AliDebugClass(2, "Failed check on DCA to primary vertes");
	218	return kFALSE;
	219	}
	220	if (TMath::Abs(v0->GetV0CosineOfPointingAngle()) < fMinCosPointAngle) {
	221	AliDebugClass(2, "Failed check on cosine of pointing angle");
	222	return kFALSE;
	223	}
	224	if (TMath::Abs(v0->GetDcaV0Daughters()) > fMaxDaughtersDCA) {
	225	AliDebugClass(2, "Failed check on DCA between daughters");
	226	return kFALSE;
	227	}
3da8cef7	228	if (TMath::Abs(v0->Y(fHypothesis)) > fMaxRapidity) {
e6952ec7	229	AliDebugClass(2, "Failed check on V0 rapidity");
e6952ec7	230	return kFALSE;
76d89202	231	}
e6952ec7	232
e6952ec7	233
76d89202	234	// check PID on proton or antiproton from V0
e6952ec7	235
	236	// check initialization of PID object
	237	AliPIDResponse *pid = fEvent->GetPIDResponse();
	238	if (!pid) {
	239	AliFatal("NULL PID response");
	240	return kFALSE;
	241	}
	242
	243	// check if TOF is matched
	244	// and computes all values used in the PID cut
	245	//Bool_t isTOFpos = MatchTOF(ptrack);
	246	//Bool_t isTOFneg = MatchTOF(ntrack);
	247	Double_t pospTPC = pTrack->GetTPCmomentum();
	248	Double_t negpTPC = nTrack->GetTPCmomentum();
	249	//Double_t posp = pTrack->P();
	250	//Double_t negp = nTrack->P();
	251	Double_t posnsTPC = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID));
	252	Double_t posnsTPC2 = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID2));
	253	//Double_t posnsTOF = TMath::Abs(pid->NumberOfSigmasTOF(ptrack, fPID));
	254	Double_t negnsTPC = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID));
	255	Double_t negnsTPC2 = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID2));
	256	//Double_t negnsTOF = TMath::Abs(pid->NumberOfSigmasTOF(ntrack, fPID));
	257	Double_t maxTPC = 1E20;
	258	Double_t maxTPC2 = 1E20;
	259	//Double_t maxTOF = 1E20;
	260
	261	// applies the cut differently depending on the PID and the momentum
	262
	263	if(fHypothesis==kLambda0) {
	264	//if (isTOFpos) {
	265	// TPC: 5sigma cut for all
	266	//if (posnsTPC > 5.0) return kFALSE;
	267	// TOF: 3sigma
	268	// maxTOF = 3.0;
	269	//return (posnsTOF <= maxTOF);
	270	//} else {
	271	// TPC:
	272	// below 600 MeV: 4sigma
	273	// above 600 MeV: 3sigma
	274
	275	if (pospTPC <= 0.6 && fPID==AliPID::kProton)
	276	maxTPC = fPIDCut1;
	277	else if (pospTPC > 0.6 && fPID==AliPID::kProton)
	278	maxTPC = fPIDCut2;
	279	else
	280	return kFALSE;
	281
	282	maxTPC2 = fPIDCut3;
	283
	284	if (! ((posnsTPC <= maxTPC) && (negnsTPC2 <= maxTPC2)) ) {
	285	AliDebugClass(2, "Failed check on V0 PID");
	286	return kFALSE;
	287	}
	288	}
	289
	290	//}
	291
	292	if(fHypothesis==kLambda0Bar) {
	293	//if (isTOFneg) {
	294	// TPC: 5sigma cut for all
	295	//if (negnsTPC > 5.0) return kFALSE;
	296	// TOF: 3sigma
	297	// maxTOF = 3.0;
	298	//return (negnsTOF <= maxTOF);
299	//} else {
300	// TPC:
301	// below 600 MeV: 4sigma
302	// above 600 MeV: 3sigma
303
304	if (negpTPC <= 0.6 && fPID==AliPID::kProton)
305	maxTPC = fPIDCut1;
306	else if (negpTPC > 0.6 && fPID==AliPID::kProton)
307	maxTPC = fPIDCut2;
308	else
309	return kFALSE;
310
311	maxTPC2 = fPIDCut3;
312
313	if(! ((negnsTPC <= maxTPC) && (posnsTPC2 <= maxTPC2)) ) {
314	AliDebugClass(2, "Failed check on V0 PID");
315	return kFALSE;
316	}
317	}
318	//}
319
320
76d89202	321	// if we reach this point, all checks were successful
76d89202	322	AliDebugClass(2, "Good V0 (hallelujah)");
e6952ec7	323	return kTRUE;
76d89202	324	}
	325
	326	//_________________________________________________________________________________________________
	327	Bool_t AliRsnCutV0::CheckAOD(AliAODv0 *v0)
	328	{
	329	//
	330	// Check an AOD V0.
	331	// This is done doing directly all checks, since there is not
	332	// an equivalend checker for AOD tracks
	333	//
	334
	335	AliDebugClass(2, "Check AOD");
	336	if (v0->GetOnFlyStatus()) {
	337	AliDebugClass(2, "Rejecting V0 in 'on fly' status");
	338	return kFALSE; // if kTRUE, then this V0 is recontructed
	339	}
	340
e6952ec7	341
e6952ec7	342
76d89202	343	// retrieve pointer to owner event
	344	AliAODEvent *lAODEvent = fEvent->GetRefAOD();
	345	Double_t xPrimaryVertex = lAODEvent->GetPrimaryVertex()->GetX();
	346	Double_t yPrimaryVertex = lAODEvent->GetPrimaryVertex()->GetY();
	347	Double_t zPrimaryVertex = lAODEvent->GetPrimaryVertex()->GetZ();
	348	AliDebugClass(2, Form("Primary vertex: %f %f %f", xPrimaryVertex, yPrimaryVertex, zPrimaryVertex));
	349
	350	// retrieve the V0 daughters
	351	AliAODTrack pTrack = (AliAODTrack ) (v0->GetSecondaryVtx()->GetDaughter(0));
	352	AliAODTrack nTrack = (AliAODTrack ) (v0->GetSecondaryVtx()->GetDaughter(1));
	353
	354
	355	// check quality cuts
	356	UInt_t filtermapP = 9999;
	357	UInt_t filtermapN = 9999;
	358	filtermapP = pTrack->GetFilterMap();
	359	filtermapN = nTrack->GetFilterMap();
e6952ec7	360
e6952ec7	361
76d89202	362	if ( !pTrack->TestFilterBit(fAODTestFilterBit) ) {
e6952ec7	363	AliDebugClass(2, Form("Positive daughter failed quality cuts filtermapP=%d",filtermapP));
e6952ec7	364	return kFALSE;
76d89202	365	}
76d89202	366	if ( !nTrack->TestFilterBit(fAODTestFilterBit) ) {
e6952ec7	367	AliDebugClass(2, Form("Negative daughter failed quality cuts filtermapN=%d",filtermapN));
e6952ec7	368	return kFALSE;
76d89202	369	}
e6952ec7	370
e6952ec7	371
76d89202	372
1c7bcf66	373	/*AliDebugClass(1, Form("fESDtrackCuts=%p",fESDtrackCuts));
76d89202	374	if (fESDtrackCuts) { // use fESDtrackCuts to retrieve cuts values
76d89202	375
e6952ec7	376
	377	AliDebugClass(2, "Checking quality cuts");
	378
	379
	380
	381	const Bool_t CutAcceptKinkDaughters = fESDtrackCuts->GetAcceptKinkDaughters(); // 0 = kFalse
	382	const Float_t CutMaxChi2PerClusterTPC = fESDtrackCuts->GetMaxChi2PerClusterTPC();
	383	const Int_t CutMinNClusterTPC = fESDtrackCuts->GetMinNClusterTPC();
	384	const Bool_t CutRequireTPCRefit = fESDtrackCuts->GetRequireTPCRefit();
	385	//AliDebugClass(2, Form("accept kink=%d maxchi2=%f minnclSTPC=%d requireTPCrefit=%d", CutAcceptKinkDaughters,
	386	// CutMaxChi2PerClusterTPC, CutMinNClusterTPC, CutRequireTPCRefit));
	387	//AliDebugClass(2, Form("pTrack->TPCNcls=%d", pTrack->GetTPCNcls() ));
	388	//AliDebugClass(2, Form("nTrack->TPCNcls=%d", nTrack->GetTPCNcls() ));
	389
	390
76d89202	391	if(pTrack->GetTPCNcls() < CutMinNClusterTPC) {AliDebugClass(2, "Positive daughter not MinNclsTPC"); return kFALSE;}
	392	if(nTrack->GetTPCNcls() < CutMinNClusterTPC) {AliDebugClass(2, "Negative daughter not MinNclsTPC"); return kFALSE;}
	393
1c7bcf66	394	}*/
76d89202	395
76d89202	396	// filter like-sign V0
e6952ec7	397	if ( TMath::Abs( ((pTrack->Charge()) - (nTrack->Charge())) ) < 0.1) {
	398	AliDebugClass(2, "Failed like-sign V0 check");
	399	return kFALSE;
	400	}
76d89202	401
76d89202	402	// check compatibility with expected species hypothesis
e6952ec7	403	Double_t mass = 0.0;
	404	if(fHypothesis==kLambda0) {
	405	mass = v0->MassLambda();
	406	}
	407	else if (fHypothesis==kLambda0Bar) {
	408	mass = v0->MassAntiLambda();
	409	}
	410	if ((TMath::Abs(mass - fMass)) > fTolerance) {
	411	AliDebugClass(2, Form("V0 is not in the expected inv mass range Mass: %d %f %f", fHypothesis, fMass, mass));
	412	return kFALSE;
	413	}
76d89202	414	AliDebugClass(2, Form("Mass: %d %f %f", fHypothesis, fMass, mass));
e6952ec7	415
	416
	417	// topological checks
	418	if (TMath::Abs(v0->DcaV0ToPrimVertex()) > fMaxDCAVertex) {
	419	AliDebugClass(2, Form("Failed check on DCA to primary vertes dca=%f maxdca=%f",TMath::Abs(v0->DcaV0ToPrimVertex()),fMaxDCAVertex));
	420	return kFALSE;
	421	}
	422	if (TMath::Abs( TMath::Cos(v0->OpenAngleV0()) ) < fMinCosPointAngle) {
76d89202	423	AliDebugClass(2, "Failed check on cosine of pointing angle");
3da8cef7	424	return kFALSE;
3da8cef7	425	}
76d89202	426	if (TMath::Abs(v0->DcaV0Daughters()) > fMaxDaughtersDCA) {
	427	AliDebugClass(2, "Failed check on DCA between daughters");
	428	return kFALSE;
	429	}
	430	if (TMath::Abs(v0->RapLambda()) > fMaxRapidity) {
e6952ec7	431	AliDebugClass(2, "Failed check on V0 rapidity");
e6952ec7	432	return kFALSE;
76d89202	433	}
3da8cef7	434
3da8cef7	435
76d89202	436
e6952ec7	437
3da8cef7	438
	439	// check initialization of PID object
	440	AliPIDResponse *pid = fEvent->GetPIDResponse();
	441	if (!pid) {
e6952ec7	442	AliFatal("NULL PID response");
e6952ec7	443	return kFALSE;
3da8cef7	444	}
	445
	446	// check if TOF is matched
	447	// and computes all values used in the PID cut
	448	//Bool_t isTOFpos = MatchTOF(ptrack);
	449	//Bool_t isTOFneg = MatchTOF(ntrack);
	450	Double_t pospTPC = pTrack->GetTPCmomentum();
	451	Double_t negpTPC = nTrack->GetTPCmomentum();
	452	//Double_t posp = pTrack->P();
	453	//Double_t negp = nTrack->P();
	454	Double_t posnsTPC = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID));
	455	Double_t posnsTPC2 = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID2));
	456	//Double_t posnsTOF = TMath::Abs(pid->NumberOfSigmasTOF(ptrack, fPID));
	457	Double_t negnsTPC = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID));
	458	Double_t negnsTPC2 = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID2));
	459	//Double_t negnsTOF = TMath::Abs(pid->NumberOfSigmasTOF(ntrack, fPID));
	460	Double_t maxTPC = 1E20;
	461	Double_t maxTPC2 = 1E20;
	462	//Double_t maxTOF = 1E20;
e6952ec7	463
3da8cef7	464	// applies the cut differently depending on the PID and the momentum
e6952ec7	465
	466	if(fHypothesis==kLambda0) {
	467
	468
	469	//if (isTOFpos) {
	470	// TPC: 5sigma cut for all
	471	//if (posnsTPC > 5.0) return kFALSE;
	472	// TOF: 3sigma
	473	// maxTOF = 3.0;
	474	//return (posnsTOF <= maxTOF);
	475	//} else {
	476	// TPC:
	477	// below 600 MeV: 4sigma
	478	// above 600 MeV: 3sigma
	479
	480	if (pospTPC <= 0.6 && fPID==AliPID::kProton)
	481	maxTPC = fPIDCut1; // EF safer value to run on MC
	482	else if (pospTPC > 0.6 && fPID==AliPID::kProton)
	483	maxTPC = fPIDCut2; // EF safer value to run on MC
	484	else
	485	return kFALSE;
	486
	487	maxTPC2 = fPIDCut3; // EF safer value to run on MC
	488
	489	if (! ((posnsTPC <= maxTPC) && (negnsTPC2 <= maxTPC2)) ) {
	490	AliDebugClass(2, "Failed check on V0 PID");
	491	return kFALSE;
	492	}
3da8cef7	493	}
e6952ec7	494
	495	if(fHypothesis==kLambda0Bar) {
	496
	497	//if (isTOFneg) {
	498	// TPC: 5sigma cut for all
	499	//if (negnsTPC > 5.0) return kFALSE;
	500	// TOF: 3sigma
	501	// maxTOF = 3.0;
	502	//return (negnsTOF <= maxTOF);
	503	//} else {
	504	// TPC:
	505	// below 600 MeV: 4sigma
	506	// above 600 MeV: 3sigma
	507
	508	if (negpTPC <= 0.6 && fPID==AliPID::kProton)
	509	maxTPC = fPIDCut1; // EF safer value to run on MC
	510	else if (negpTPC > 0.6 && fPID==AliPID::kProton)
	511	maxTPC = fPIDCut2; // EF safer value to run on MC
	512	else
	513	return kFALSE;
	514
	515	maxTPC2 = fPIDCut3; // EF safer value to run on MC
	516
	517	if(! ((negnsTPC <= maxTPC) && (posnsTPC2 <= maxTPC2)) ) {
	518	AliDebugClass(2, "Failed check on V0 PID");
	519	return kFALSE;
	520	}
3da8cef7	521	}
e6952ec7	522
e6952ec7	523
61f275d1	524
8082e538	525	// if we reach this point, all checks were successful
76d89202	526	AliDebugClass(1, "Good AOD V0 (hallelujah)");
76d89202	527	AliDebugClass(1, Form("Mass: %d %f %f %d %d", fHypothesis, fMass, mass, filtermapP, filtermapN));
e6952ec7	528	return kTRUE;
61f275d1	529
8082e538	530	}
	531
	532	//_________________________________________________________________________________________________
	533	void AliRsnCutV0::Print(const Option_t *) const
	534	{
	535	//
	536	// Print information on this cut
	537	//
	538	}