#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, 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),
fPIDCutProton(0),
fPIDCutPion(0),
fESDtrackCuts(0x0),
fCutQuality(Form("%sDaughtersQuality", name)),
fAODTestFilterBit(5)
{
SetHypothesis(hypothesis);
}
AliRsnCutV0::AliRsnCutV0(const AliRsnCutV0 ©) :
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),
fPIDCutProton(copy.fPIDCutProton),
fPIDCutPion(copy.fPIDCutPion),
fESDtrackCuts(copy.fESDtrackCuts),
fCutQuality(copy.fCutQuality),
fAODTestFilterBit(copy.fAODTestFilterBit)
{
fCutQuality.SetPtRange(0.15, 1E+20);
fCutQuality.SetEtaRange(-0.8, 0.8);
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 ©)
{
if (this == ©)
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;
fCutQuality = copy.fCutQuality;
fPID = copy.fPID;
fPID2 = copy.fPID2;
fPIDCutProton = copy.fPIDCutProton;
fPIDCutPion = copy.fPIDCutPion;
fESDtrackCuts = copy.fESDtrackCuts;
fCutQuality = copy.fCutQuality;
fAODTestFilterBit = copy.fAODTestFilterBit;
return (*this);
}
Bool_t AliRsnCutV0::IsSelected(TObject *object)
{
if (!TargetOK(object)) return kFALSE;
AliESDv0 *v0esd = fDaughter->Ref2ESDv0();
AliAODv0 *v0aod = fDaughter->Ref2AODv0();
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)
{
AliDebugClass(1, "Check ESD");
if (v0->GetOnFlyStatus()) {
AliDebugClass(1, "Rejecting V0 in 'on fly' status");
return kFALSE;
}
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));
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);
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;
}
}
if ( TMath::Abs( ((pTrack->GetSign()) - (nTrack->GetSign())) ) < 0.1) {
AliDebugClass(2, "Failed like-sign V0 check");
return kFALSE;
}
v0->ChangeMassHypothesis(fHypothesis);
if ((TMath::Abs(v0->GetEffMass() - fMass)) > fTolerance) {
AliDebugClass(2, "V0 is not in the expected inv mass range");
return kFALSE;
}
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;
}
Double_t v0Position[3];
v0->GetXYZ(v0Position[0],v0Position[1],v0Position[2]);
Double_t radius = TMath::Sqrt(TMath::Power(v0Position[0],2) + TMath::Power(v0Position[1],2));
if ( ( radius < 0.8 ) || ( radius > 100 ) ) {
AliDebugClass(2, "Failed fiducial volume");
return kFALSE;
}
AliPIDResponse *pid = fEvent->GetPIDResponse();
if (!pid) {
AliFatal("NULL PID response");
return kFALSE;
}
Double_t posnsTPC = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID));
Double_t posnsTPC2 = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID2));
Double_t negnsTPC = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID));
Double_t negnsTPC2 = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID2));
Double_t maxTPC = 1E20;
Double_t maxTPC2 = 1E20;
if(fHypothesis==kLambda0) {
maxTPC = fPIDCutProton;
maxTPC2 = fPIDCutPion;
if (! ((posnsTPC <= maxTPC) && (negnsTPC2 <= maxTPC2)) ) {
AliDebugClass(2, "Failed check on V0 PID");
return kFALSE;
}
}
if(fHypothesis==kLambda0Bar) {
maxTPC = fPIDCutProton;
maxTPC2 = fPIDCutPion;
if(! ((negnsTPC <= maxTPC) && (posnsTPC2 <= maxTPC2)) ) {
AliDebugClass(2, "Failed check on V0 PID");
return kFALSE;
}
}
AliDebugClass(2, "Good V0");
return kTRUE;
}
Bool_t AliRsnCutV0::CheckAOD(AliAODv0 *v0)
{
AliDebugClass(2, "Check AOD");
if (v0->GetOnFlyStatus()) {
AliDebugClass(2, "Rejecting V0 in 'on fly' status");
return kFALSE;
}
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));
AliAODTrack *pTrack = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
AliAODTrack *nTrack = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
UInt_t filtermapP = 9999;
UInt_t filtermapN = 9999;
filtermapP = pTrack->GetFilterMap();
filtermapN = nTrack->GetFilterMap();
if ( TMath::Abs( ((pTrack->Charge()) - (nTrack->Charge())) ) < 0.1) {
AliDebugClass(2, "Failed like-sign V0 check");
return kFALSE;
}
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));
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;
}
AliAODVertex *vertex = lAODEvent->GetPrimaryVertex();
Double_t cospointangle = v0->CosPointingAngle(vertex);
if (TMath::Abs( cospointangle ) < 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;
}
Double_t radius = v0->RadiusV0();
if ( ( radius < 0.8 ) || ( radius > 100 ) ) {
AliDebugClass(2, "Failed fiducial volume");
return kFALSE;
}
AliPIDResponse *pid = fEvent->GetPIDResponse();
if (!pid) {
AliFatal("NULL PID response");
return kFALSE;
}
Double_t posnsTPC = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID));
Double_t posnsTPC2 = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID2));
Double_t negnsTPC = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID));
Double_t negnsTPC2 = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID2));
Double_t maxTPC = 1E20;
Double_t maxTPC2 = 1E20;
if(fHypothesis==kLambda0) {
maxTPC = fPIDCutProton;
maxTPC2 = fPIDCutPion;
if (! ((posnsTPC <= maxTPC) && (negnsTPC2 <= maxTPC2)) ) {
AliDebugClass(2, "Failed check on V0 PID");
return kFALSE;
}
}
if(fHypothesis==kLambda0Bar) {
maxTPC = fPIDCutProton;
maxTPC2 = fPIDCutPion;
if(! ((negnsTPC <= maxTPC) && (posnsTPC2 <= maxTPC2)) ) {
AliDebugClass(2, "Failed check on V0 PID");
return kFALSE;
}
}
AliDebugClass(1, "Good AOD V0");
AliDebugClass(1, Form("Mass: %d %f %f %d %d", fHypothesis, fMass, mass, filtermapP, filtermapN));
return kTRUE;
}
void AliRsnCutV0::Print(const Option_t *) const
{
}