[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAODRecoCascadeHF.cxx

/**************************************************************************
 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

/////////////////////////////////////////////////////////////
//
// Class for AOD reconstructed heavy-flavour cascades
//
// Author: X-M. Zhang, zhangxm@iopp.ccnu.edu.cn
/////////////////////////////////////////////////////////////

#include <TVector3.h>
#include <TDatabasePDG.h>
#include <TClonesArray.h>
#include "AliAODMCParticle.h"
#include "AliAODRecoDecay.h"
#include "AliAODVertex.h"
#include "AliAODRecoDecayHF2Prong.h"
#include "AliAODRecoCascadeHF.h"

ClassImp(AliAODRecoCascadeHF)
//-----------------------------------------------------------------------------

AliAODRecoCascadeHF::AliAODRecoCascadeHF() :
  AliAODRecoDecayHF2Prong()
{
  //
  // Default Constructor
  //
}
//-----------------------------------------------------------------------------
AliAODRecoCascadeHF::AliAODRecoCascadeHF(AliAODVertex *vtx2, Short_t charge,
					 Double_t *px, Double_t *py, Double_t *pz,
					 Double_t *d0, Double_t *d0err, Double_t dca) :
  AliAODRecoDecayHF2Prong(vtx2, px, py, pz, d0, d0err, dca)
{
  //
  //  Constructor with AliAODVertex for decay vertex
  //
  SetCharge(charge);
}
//-----------------------------------------------------------------------------
AliAODRecoCascadeHF::AliAODRecoCascadeHF(AliAODVertex *vtx2, Short_t charge,
					 Double_t *d0, Double_t *d0err, Double_t dca) :
  AliAODRecoDecayHF2Prong(vtx2, d0, d0err, dca)
{
  //
  //  Constructor with decay vertex and without prongs momenta
  //
  SetCharge(charge);
}
//-----------------------------------------------------------------------------
AliAODRecoCascadeHF::AliAODRecoCascadeHF(const AliAODRecoCascadeHF &source) :
  AliAODRecoDecayHF2Prong(source)
{
  //
  // Copy constructor
  //
}
//-----------------------------------------------------------------------------
AliAODRecoCascadeHF &AliAODRecoCascadeHF::operator=(const AliAODRecoCascadeHF &source)
{
  //
  // assignment operator
  //
  if(&source == this) return *this;

  AliAODRecoDecayHF2Prong::operator=(source);

  return *this;
}
//-----------------------------------------------------------------------------
AliAODRecoCascadeHF::~AliAODRecoCascadeHF()
{
  //
  // Default Destructor
  //
}
//-----------------------------------------------------------------------------
Double_t AliAODRecoCascadeHF::InvMassDstarKpipi() const 
{
  //
  // 3 prong invariant mass of the D0 daughters and the soft pion
  //
  Double_t e[3];
  if (Charge()>0){
    e[0]=Get2Prong()->EProng(0,211);
    e[1]=Get2Prong()->EProng(1,321);
  }else{
    e[0]=Get2Prong()->EProng(0,321);
    e[1]=Get2Prong()->EProng(1,211);
  }
  e[2]=EProng(0,211);

  Double_t esum = e[0]+e[1]+e[2];
  Double_t minv = TMath::Sqrt(esum*esum-P()*P());

  return minv; 
}
//----------------------------------------------------------------------------
Int_t AliAODRecoCascadeHF::MatchToMC(Int_t pdgabs,Int_t pdgabs2prong,
                                     Int_t *pdgDg,Int_t *pdgDg2prong,
				     TClonesArray *mcArray, Bool_t isV0) const
{
  //
  // Check if this candidate is matched to a MC signal
  // If no, return -1
  // If yes, return label (>=0) of the AliAODMCParticle
  // 

  Int_t ndg=GetNDaughters();
  if(ndg==0) {
    AliError("No daughters available");
    return -1;
  }

  if ( isV0 &&
       ( (pdgDg[1]==2212 && pdgDg[0]==310) ||
	 (pdgDg[1]==211 && pdgDg[0]==3122) ) ) {
    AliWarning("Please, pay attention: first element in AliAODRecoCascadeHF object must be the bachelor and second one V0. Skipping!");
    return -1;
  }

  Int_t lab2Prong = -1;

  if (!isV0) {
    AliAODRecoDecayHF2Prong *the2Prong = Get2Prong();
    lab2Prong = the2Prong->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong);
  } else {
    AliAODv0 *theV0 = dynamic_cast<AliAODv0*>(Getv0());
    lab2Prong = theV0->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong); // the V0
  }

  if(lab2Prong<0) return -1;

  Int_t dgLabels[10]={0,0,0,0,0,0,0,0,0,0};

  if (!isV0) {
    // loop on daughters and write labels
    for(Int_t i=0; i<ndg; i++) {
      AliVTrack *trk = dynamic_cast<AliVTrack*>(GetDaughter(i));
      if(!trk) continue;
      Int_t lab = trk->GetLabel();
      if(lab==-1) { // this daughter is the 2prong
	lab=lab2Prong;
      } else if(lab<-1) continue;
      dgLabels[i] = lab;
    }
  } else {
    AliVTrack *trk = dynamic_cast<AliVTrack*>(GetBachelor()); // the bachelor
    if (!trk) return -1;
    dgLabels[0] = trk->GetLabel();//TMath::Abs(trk->GetLabel());
    dgLabels[1] = lab2Prong;
  }

  Int_t finalLabel = AliAODRecoDecay::MatchToMC(pdgabs,mcArray,dgLabels,2,2,pdgDg);

  if (finalLabel>=0){
    // debug printouts for Lc->V0 bachelor case

    if ( isV0 && (dgLabels[0]!=-1 && dgLabels[1]!=-1) ) {
      AliAODv0 *theV0 = dynamic_cast<AliAODv0*>(Getv0());
      Bool_t onTheFly = theV0->GetOnFlyStatus();
      if (pdgDg[0]==2212 && pdgDg[1]==310) {
	AliAODMCParticle*k0s = dynamic_cast<AliAODMCParticle*>(mcArray->At(lab2Prong));
	if(k0s){
	  Int_t labK0 = k0s->GetMother();	
	  AliAODMCParticle*k0bar = dynamic_cast<AliAODMCParticle*>(mcArray->At(labK0));
	  if(k0bar){
	    AliDebug(1,Form(" (onTheFly=%1d) LabelV0=%d (%d) -> LabelK0S=%d (%d -> %d %d)",onTheFly,labK0,k0bar->GetPdgCode(),lab2Prong,pdgabs2prong,pdgDg2prong[0],pdgDg2prong[1]));
	    AliDebug(1,Form(" LabelLc=%d (%d) -> LabelBachelor=%d (%d) LabelV0=%d (%d)",
			    finalLabel,pdgabs,
			    dgLabels[0],pdgDg[0],dgLabels[1],pdgDg[1]));
	  }
	}
      } else if (pdgDg[0]==211 && pdgDg[1]==3122) {
	AliDebug(1,Form(" (onTheFly=%1d) LabelV0=%d (%d -> %d %d)",onTheFly,lab2Prong,pdgabs2prong,pdgDg2prong[0],pdgDg2prong[1]));
	AliDebug(1,Form(" LabelLc=%d (%d) -> LabelBachelor=%d (%d) LabelV0=%d (%d)",
			finalLabel,pdgabs,
		      dgLabels[0],pdgDg[0],dgLabels[1],pdgDg[1]));
      }

    }
  }

  return finalLabel;

}
//-----------------------------------------------------------------------------
Bool_t AliAODRecoCascadeHF::SelectDstar(const Double_t *cutsDstar,
					const Double_t *cutsD0,
					Bool_t testD0) const
{
  //
  // cutsDstar[0] = inv. mass half width of D* [GeV]
  // cutsDstar[1] = half width of (M_Kpipi-M_D0) [GeV]
  // cutsDstar[2] = PtMin of pi_s [GeV/c]
  // cutsDstar[3] = PtMax of pi_s [GeV/c]
  // cutsDstar[4] = theta, angle between the pi_s and decay plane of the D0 [rad]
  //
  // cutsD0[0] = inv. mass half width [GeV]   
  // cutsD0[1] = dca [cm]
  // cutsD0[2] = cosThetaStar 
  // cutsD0[3] = pTK [GeV/c]
  // cutsD0[4] = pTPi [GeV/c]
  // cutsD0[5] = d0K [cm]   upper limit!
  // cutsD0[6] = d0Pi [cm]  upper limit!
  // cutsD0[7] = d0d0 [cm^2]
  // cutsD0[8] = cosThetaPoint


  // check that the D0 passes the cuts
  // (if we have a D*+, it has to pass as D0, 
  //  if we have a D*-, it has to pass as D0bar)

  if(testD0) {
    Int_t okD0=0,okD0bar=0;
    Get2Prong()->SelectD0(cutsD0,okD0,okD0bar);
    if((Charge()==+1 && !okD0) || (Charge()==-1 && !okD0bar)) return kFALSE; 
  }
 
  if( (PtProng(0)<cutsDstar[2]) || (PtProng(0)>cutsDstar[3]) ) return kFALSE;

  Double_t mDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
  Double_t invmDstar = InvMassDstarKpipi();
  if(TMath::Abs(mDstar-invmDstar)>cutsDstar[0]) return kFALSE;

  Double_t mD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
  if(TMath::Abs((mDstar-mD0)-DeltaInvMass())>cutsDstar[1]) return kFALSE;

  Double_t theta = AngleD0dkpPisoft(); 
  if(theta>cutsDstar[4]) return kFALSE;
  
  return kTRUE;
}
//-----------------------------------------------------------------------------
Bool_t AliAODRecoCascadeHF::SelectLctoV0(const Double_t *cutsLctoV0, 
					 Bool_t okLck0sp, Bool_t okLcLpi, Bool_t okLcLbarpi) const 
{
  // cuts on Lambdac candidates to V0+bachelor
  // (to be passed to AliAODRecoDecayHF3Prong::SelectLctoV0())
  // 0 = inv. mass half width in K0s hypothesis [GeV]   
  // 1 = inv. mass half width in Lambda hypothesis [GeV]   
  // 2 = inv. mass V0 in K0s hypothesis half width [GeV]   
  // 3 = inv. mass V0 in Lambda hypothesis half width [GeV]   
  // 4 = pT min Bachelor track [GeV/c]
  // 5 = pT min V0-Positive track [GeV/c]
  // 6 = pT min V0-Negative track [GeV/c]
  // 7 = dca cut on the cascade (cm)
  // 8 = dca cut on the V0 (cm)

  //   if ( !Getv0() || !Getv0PositiveTrack() || !Getv0NegativeTrack() ) 
  //     { AliInfo(Form("Not adapted for ESDv0s, return true...")); return false; }

  Double_t mLck0sp,mLcLpi;
  okLck0sp=1; okLcLpi=1; okLcLbarpi=1;
  
  Double_t mLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
  Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
  Double_t mLPDG = TDatabasePDG::Instance()->GetParticle(3122)->Mass();

  // k0s + p
  double mk0s = Getv0()->MassK0Short();
  mLck0sp = InvMassLctoK0sP();

  // lambda + pi 
  double mlambda = Getv0()->MassLambda();
  double malambda = Getv0()->MassAntiLambda();
  mLcLpi = InvMassLctoLambdaPi();

  // cut on Lc mass
  //   with k0s p hypothesis
  if(TMath::Abs(mLck0sp-mLcPDG)>cutsLctoV0[0]) okLck0sp = 0;
  //   with Lambda pi hypothesis
  if(TMath::Abs(mLcLpi-mLcPDG)>cutsLctoV0[1]) okLcLpi = 0;
  okLcLbarpi = okLcLpi;
 
  // cuts on the v0 mass
  if( TMath::Abs(mk0s-mk0sPDG)>cutsLctoV0[2]) okLck0sp = 0;
  //if( TMath::Abs(mlambda-mLPDG)>cutsLctoV0[3] && 
  //TMath::Abs(malambda-mLPDG)>cutsLctoV0[3] ) okLcLpi = 0;
  if( !(GetBachelor()->Charge()==+1 && TMath::Abs(mlambda-mLPDG)<=cutsLctoV0[3]) ) okLcLpi = 0;
  if( !(GetBachelor()->Charge()==-1 && TMath::Abs(malambda-mLPDG)<=cutsLctoV0[3]) ) okLcLbarpi = 0;
  
  if(!okLck0sp && !okLcLpi && !okLcLbarpi) return 0;
  
  // cuts on the minimum pt of the tracks 
  if(TMath::Abs(GetBachelor()->Pt()) < cutsLctoV0[4]) return 0;
  if(TMath::Abs(Getv0PositiveTrack()->Pt()) < cutsLctoV0[5]) return 0;
  if(TMath::Abs(Getv0NegativeTrack()->Pt()) < cutsLctoV0[6]) return 0;
  
  // cut on the cascade dca
  if( TMath::Abs(GetDCA(0))>cutsLctoV0[7] //||
      //TMath::Abs(Getv0()->DcaPosToPrimVertex())>cutsLctoV0[7] ||
      //TMath::Abs(Getv0()->DcaNegToPrimVertex())>cutsLctoV0[7]
      ) return 0;
  
  // cut on the v0 dca
  if(TMath::Abs(Getv0()->DcaV0Daughters()) > cutsLctoV0[8]) return 0;

  // cut on V0 cosine of pointing angle wrt PV
  if (CosV0PointingAngle() < cutsLctoV0[9]) { // cosine of V0 pointing angle wrt primary vertex
    AliDebug(4,Form(" V0 cosine of pointing angle doesn't pass the cut"));
    return 0;
  }

  // cut on bachelor transverse impact parameter wrt PV
  if (TMath::Abs(Getd0Prong(0)) > cutsLctoV0[10]) { // bachelor transverse impact parameter wrt PV
    AliDebug(4,Form(" bachelor transverse impact parameter doesn't pass the cut"));
    return 0;
  }

  // cut on V0 transverse impact parameter wrt PV
  if (TMath::Abs(Getd0Prong(1)) > cutsLctoV0[11]) { // V0 transverse impact parameter wrt PV
    AliDebug(4,Form(" V0 transverse impact parameter doesn't pass the cut"));
    return 0;
  }

  // cut on K0S invariant mass veto
  if (TMath::Abs(Getv0()->MassK0Short()-mk0sPDG) < cutsLctoV0[12]) { // K0S invariant mass veto
    AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
    return 0;
  }

  // cut on Lambda/LambdaBar invariant mass veto
  if (TMath::Abs(Getv0()->MassLambda()-mLPDG) < cutsLctoV0[13] ||
      TMath::Abs(Getv0()->MassAntiLambda()-mLPDG) < cutsLctoV0[13] ) { // Lambda/LambdaBar invariant mass veto
    AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
    return 0;
  }

  // cut on gamma invariant mass veto                                                                                                                      
  if (Getv0()->InvMass2Prongs(0,1,11,11) < cutsLctoV0[14]) { // K0S invariant mass veto
    AliDebug(4,Form(" veto on gamma invariant mass doesn't pass the cut"));
    return 0;
  }

  // cut on V0 pT min                                                                                                                                      
  if (Getv0()->Pt() < cutsLctoV0[15]) { // V0 pT min                                                                                         
    AliDebug(4,Form(" V0 track Pt=%2.2e > %2.2e",Getv0()->Pt(),cutsLctoV0[15]));
    return 0;
  }
  
  return true; 

}
//-----------------------------------------------------------------------------
Double_t AliAODRecoCascadeHF::AngleD0dkpPisoft() const {
  //
  // Angle of soft pion to D0 decay plane
  // 

  TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
  TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
  TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s

  TVector3 perp = p3Trk0.Cross(p3Trk1);
  Double_t theta = p3Trk2.Angle(perp);
  if(theta>(TMath::Pi()-theta)) theta = TMath::Pi() - theta;
  theta = TMath::Pi()/2. - theta;

  return theta;
}
//-----------------------------------------------------------------------------
Bool_t AliAODRecoCascadeHF::TrigonometricalCut() const {
  //  
  // Trigonometrical constraint
  //
  TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
  TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
  TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s

  Double_t alpha = p3Trk0.Angle(p3Trk2);
  Double_t beta = p3Trk1.Angle(p3Trk2);

  Double_t cosphi01 = TMath::Cos(alpha) / TMath::Cos(AngleD0dkpPisoft());
  Double_t cosphi02 = TMath::Cos(beta) / TMath::Cos(AngleD0dkpPisoft());

  Double_t phi01 = TMath::ACos(cosphi01);
  Double_t phi02 = TMath::ACos(cosphi02);
  Double_t phi00 = p3Trk0.Angle(p3Trk1);

  if((phi01>phi00) || (phi02>phi00)) return kFALSE;
  return kTRUE;
}

//-----------------------------------------------------------------------------
Double_t AliAODRecoCascadeHF::DecayLengthV0() const
{
  //
  // Returns V0 decay length wrt primary vertex
  //

  AliAODv0 *v0 = (AliAODv0*)Getv0();

  if (!v0) 
    return -1.;
  AliAODVertex *vtxPrimary = GetPrimaryVtx();
  Double_t posVtx[3] = {0.,0.,0.};
  vtxPrimary->GetXYZ(posVtx);
  return v0->DecayLengthV0(posVtx);

}
//-----------------------------------------------------------------------------
Double_t AliAODRecoCascadeHF::DecayLengthXYV0() const
{
  //
  // Returns transverse V0 decay length wrt primary vertex
  //
  AliAODv0 *v0 = (AliAODv0*)Getv0();

  if (!v0) 
    return -1.;
  AliAODVertex *vtxPrimary = GetPrimaryVtx();
  Double_t posVtx[3] = {0.,0.,0.};
  vtxPrimary->GetXYZ(posVtx);
  return v0->DecayLengthXY(posVtx);

}
//-----------------------------------------------------------------------------
Double_t AliAODRecoCascadeHF::CosV0PointingAngle() const 
{
  //
  // Returns cosine of V0 pointing angle wrt primary vertex
  //

  AliAODv0 *v0 = (AliAODv0*)Getv0();

  if (!v0) 
    return -999.;

  AliAODVertex *vtxPrimary = GetPrimaryVtx();
  Double_t posVtx[3] = {0.,0.,0.};
  vtxPrimary->GetXYZ(posVtx);
  return v0->CosPointingAngle(posVtx);

}
//-----------------------------------------------------------------------------
Double_t AliAODRecoCascadeHF::CosV0PointingAngleXY() const 
{
  //
  // Returns XY cosine of V0 pointing angle wrt primary vertex
  //

  AliAODv0 *v0 = (AliAODv0*)Getv0();

  if (!v0) 
    return -999.;

  AliAODVertex *vtxPrimary = GetPrimaryVtx();
  Double_t posVtx[3] = {0.,0.,0.};
  vtxPrimary->GetXYZ(posVtx);
  return v0->CosPointingAngleXY(posVtx);

}
//-----------------------------------------------------------------------------
Double_t AliAODRecoCascadeHF::NormalizedV0DecayLength() const
{
  //
  // Returns V0 normalized decay length wrt primary vertex
  //

  AliAODv0 *v0 = (AliAODv0*)Getv0();

  if (!v0) 
    return -1.;
  //AliAODVertex *vtxPrimary = GetPrimaryVtx();
  //Double_t posVtx[3] = {0.,0.,0.};
  //vtxPrimary->GetXYZ(posVtx);
  //return v0->NormalizedDecayLength(posVtx);
  return v0->NormalizedDecayLength(GetPrimaryVtx());

}
//-----------------------------------------------------------------------------
Double_t AliAODRecoCascadeHF::NormalizedV0DecayLengthXY() const
{
  //
  // Returns transverse V0 normalized decay length wrt primary vertex
  //
  AliAODv0 *v0 = (AliAODv0*)Getv0();

  if (!v0) 
    return -1.;
  //AliAODVertex *vtxPrimary = GetPrimaryVtx();
  //Double_t posVtx[3] = {0.,0.,0.};
  //vtxPrimary->GetXYZ(posVtx);
  //return v0->NormalizedDecayLengthXY(posVtx);
  return v0->NormalizedDecayLengthXY(GetPrimaryVtx());

}
Commit	Line	Data
05d80dd6	1	/**************************************************************************
	2	* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
27de2dfb	16	/* $Id$ */
27de2dfb	17
05d80dd6	18	/////////////////////////////////////////////////////////////
	19	//
	20	// Class for AOD reconstructed heavy-flavour cascades
	21	//
c3fe1eed	22	// Author: X-M. Zhang, zhangxm@iopp.ccnu.edu.cn
05d80dd6	23	/////////////////////////////////////////////////////////////
	24
	25	#include <TVector3.h>
	26	#include <TDatabasePDG.h>
b3999999	27	#include <TClonesArray.h>
b3999999	28	#include "AliAODMCParticle.h"
05d80dd6	29	#include "AliAODRecoDecay.h"
	30	#include "AliAODVertex.h"
	31	#include "AliAODRecoDecayHF2Prong.h"
	32	#include "AliAODRecoCascadeHF.h"
	33
	34	ClassImp(AliAODRecoCascadeHF)
	35	//-----------------------------------------------------------------------------
	36
	37	AliAODRecoCascadeHF::AliAODRecoCascadeHF() :
dcfa35b3	38	AliAODRecoDecayHF2Prong()
05d80dd6	39	{
	40	//
	41	// Default Constructor
	42	//
	43	}
	44	//-----------------------------------------------------------------------------
	45	AliAODRecoCascadeHF::AliAODRecoCascadeHF(AliAODVertex *vtx2, Short_t charge,
	46	Double_t px, Double_t py, Double_t *pz,
	47	Double_t d0, Double_t d0err, Double_t dca) :
dcfa35b3	48	AliAODRecoDecayHF2Prong(vtx2, px, py, pz, d0, d0err, dca)
05d80dd6	49	{
	50	//
	51	// Constructor with AliAODVertex for decay vertex
	52	//
	53	SetCharge(charge);
	54	}
	55	//-----------------------------------------------------------------------------
	56	AliAODRecoCascadeHF::AliAODRecoCascadeHF(AliAODVertex *vtx2, Short_t charge,
	57	Double_t d0, Double_t d0err, Double_t dca) :
dcfa35b3	58	AliAODRecoDecayHF2Prong(vtx2, d0, d0err, dca)
05d80dd6	59	{
	60	//
	61	// Constructor with decay vertex and without prongs momenta
	62	//
	63	SetCharge(charge);
	64	}
	65	//-----------------------------------------------------------------------------
	66	AliAODRecoCascadeHF::AliAODRecoCascadeHF(const AliAODRecoCascadeHF &source) :
dcfa35b3	67	AliAODRecoDecayHF2Prong(source)
05d80dd6	68	{
	69	//
	70	// Copy constructor
	71	//
	72	}
	73	//-----------------------------------------------------------------------------
	74	AliAODRecoCascadeHF &AliAODRecoCascadeHF::operator=(const AliAODRecoCascadeHF &source)
	75	{
	76	//
	77	// assignment operator
	78	//
	79	if(&source == this) return *this;
	80
	81	AliAODRecoDecayHF2Prong::operator=(source);
	82
05d80dd6	83	return *this;
	84	}
	85	//-----------------------------------------------------------------------------
	86	AliAODRecoCascadeHF::~AliAODRecoCascadeHF()
	87	{
	88	//
	89	// Default Destructor
	90	//
	91	}
	92	//-----------------------------------------------------------------------------
	93	Double_t AliAODRecoCascadeHF::InvMassDstarKpipi() const
	94	{
	95	//
	96	// 3 prong invariant mass of the D0 daughters and the soft pion
	97	//
78489ccc	98	Double_t e[3];
	99	if (Charge()>0){
	100	e[0]=Get2Prong()->EProng(0,211);
	101	e[1]=Get2Prong()->EProng(1,321);
	102	}else{
	103	e[0]=Get2Prong()->EProng(0,321);
	104	e[1]=Get2Prong()->EProng(1,211);
	105	}
	106	e[2]=EProng(0,211);
05d80dd6	107
78489ccc	108	Double_t esum = e[0]+e[1]+e[2];
78489ccc	109	Double_t minv = TMath::Sqrt(esumesum-P()P());
05d80dd6	110
78489ccc	111	return minv;
05d80dd6	112	}
b3999999	113	//----------------------------------------------------------------------------
b3999999	114	Int_t AliAODRecoCascadeHF::MatchToMC(Int_t pdgabs,Int_t pdgabs2prong,
383937f0	115	Int_t pdgDg,Int_t pdgDg2prong,
36462b08	116	TClonesArray *mcArray, Bool_t isV0) const
b3999999	117	{
	118	//
	119	// Check if this candidate is matched to a MC signal
	120	// If no, return -1
	121	// If yes, return label (>=0) of the AliAODMCParticle
	122	//
	123
8862e6bf	124	Int_t ndg=GetNDaughters();
53c91ed7	125	if(ndg==0) {
b3999999	126	AliError("No daughters available");
	127	return -1;
	128	}
e7af8919	129
85036182	130	if ( isV0 &&
	131	( (pdgDg[1]==2212 && pdgDg[0]==310) \|\|
	132	(pdgDg[1]==211 && pdgDg[0]==3122) ) ) {
	133	AliWarning("Please, pay attention: first element in AliAODRecoCascadeHF object must be the bachelor and second one V0. Skipping!");
	134	return -1;
	135	}
	136
36462b08	137	Int_t lab2Prong = -1;
	138
	139	if (!isV0) {
	140	AliAODRecoDecayHF2Prong *the2Prong = Get2Prong();
	141	lab2Prong = the2Prong->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong);
	142	} else {
ff12b981	143	AliAODv0 theV0 = dynamic_cast<AliAODv0>(Getv0());
c75edf3d	144	lab2Prong = theV0->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong); // the V0
36462b08	145	}
b3999999	146
	147	if(lab2Prong<0) return -1;
	148
e11ae259	149	Int_t dgLabels[10]={0,0,0,0,0,0,0,0,0,0};
b3999999	150
c75edf3d	151	if (!isV0) {
	152	// loop on daughters and write labels
	153	for(Int_t i=0; i<ndg; i++) {
	154	AliVTrack trk = dynamic_cast<AliVTrack>(GetDaughter(i));
	155	if(!trk) continue;
	156	Int_t lab = trk->GetLabel();
	157	if(lab==-1) { // this daughter is the 2prong
	158	lab=lab2Prong;
	159	} else if(lab<-1) continue;
	160	dgLabels[i] = lab;
	161	}
	162	} else {
	163	AliVTrack trk = dynamic_cast<AliVTrack>(GetBachelor()); // the bachelor
	164	if (!trk) return -1;
	165	dgLabels[0] = trk->GetLabel();//TMath::Abs(trk->GetLabel());
	166	dgLabels[1] = lab2Prong;
b3999999	167	}
b3999999	168
e7af8919	169	Int_t finalLabel = AliAODRecoDecay::MatchToMC(pdgabs,mcArray,dgLabels,2,2,pdgDg);
9ff3fa8f	170
	171	if (finalLabel>=0){
	172	// debug printouts for Lc->V0 bachelor case
	173
	174	if ( isV0 && (dgLabels[0]!=-1 && dgLabels[1]!=-1) ) {
ff12b981	175	AliAODv0 theV0 = dynamic_cast<AliAODv0>(Getv0());
9ff3fa8f	176	Bool_t onTheFly = theV0->GetOnFlyStatus();
9ff3fa8f	177	if (pdgDg[0]==2212 && pdgDg[1]==310) {
ff12b981	178	AliAODMCParticlek0s = dynamic_cast<AliAODMCParticle>(mcArray->At(lab2Prong));
4184976e	179	if(k0s){
	180	Int_t labK0 = k0s->GetMother();
	181	AliAODMCParticlek0bar = dynamic_cast<AliAODMCParticle>(mcArray->At(labK0));
	182	if(k0bar){
	183	AliDebug(1,Form(" (onTheFly=%1d) LabelV0=%d (%d) -> LabelK0S=%d (%d -> %d %d)",onTheFly,labK0,k0bar->GetPdgCode(),lab2Prong,pdgabs2prong,pdgDg2prong[0],pdgDg2prong[1]));
	184	AliDebug(1,Form(" LabelLc=%d (%d) -> LabelBachelor=%d (%d) LabelV0=%d (%d)",
	185	finalLabel,pdgabs,
	186	dgLabels[0],pdgDg[0],dgLabels[1],pdgDg[1]));
	187	}
	188	}
9ff3fa8f	189	} else if (pdgDg[0]==211 && pdgDg[1]==3122) {
	190	AliDebug(1,Form(" (onTheFly=%1d) LabelV0=%d (%d -> %d %d)",onTheFly,lab2Prong,pdgabs2prong,pdgDg2prong[0],pdgDg2prong[1]));
	191	AliDebug(1,Form(" LabelLc=%d (%d) -> LabelBachelor=%d (%d) LabelV0=%d (%d)",
	192	finalLabel,pdgabs,
e7af8919	193	dgLabels[0],pdgDg[0],dgLabels[1],pdgDg[1]));
9ff3fa8f	194	}
9ff3fa8f	195
e7af8919	196	}
	197	}
	198
	199	return finalLabel;
	200
b3999999	201	}
05d80dd6	202	//-----------------------------------------------------------------------------
	203	Bool_t AliAODRecoCascadeHF::SelectDstar(const Double_t *cutsDstar,
	204	const Double_t *cutsD0,
	205	Bool_t testD0) const
	206	{
	207	//
	208	// cutsDstar[0] = inv. mass half width of D* [GeV]
	209	// cutsDstar[1] = half width of (M_Kpipi-M_D0) [GeV]
	210	// cutsDstar[2] = PtMin of pi_s [GeV/c]
	211	// cutsDstar[3] = PtMax of pi_s [GeV/c]
	212	// cutsDstar[4] = theta, angle between the pi_s and decay plane of the D0 [rad]
	213	//
	214	// cutsD0[0] = inv. mass half width [GeV]
	215	// cutsD0[1] = dca [cm]
	216	// cutsD0[2] = cosThetaStar
	217	// cutsD0[3] = pTK [GeV/c]
	218	// cutsD0[4] = pTPi [GeV/c]
	219	// cutsD0[5] = d0K [cm] upper limit!
	220	// cutsD0[6] = d0Pi [cm] upper limit!
	221	// cutsD0[7] = d0d0 [cm^2]
	222	// cutsD0[8] = cosThetaPoint
	223
	224
	225	// check that the D0 passes the cuts
	226	// (if we have a D*+, it has to pass as D0,
	227	// if we have a D*-, it has to pass as D0bar)
	228
	229	if(testD0) {
	230	Int_t okD0=0,okD0bar=0;
	231	Get2Prong()->SelectD0(cutsD0,okD0,okD0bar);
	232	if((Charge()==+1 && !okD0) \|\| (Charge()==-1 && !okD0bar)) return kFALSE;
	233	}
	234
	235	if( (PtProng(0)<cutsDstar[2]) \|\| (PtProng(0)>cutsDstar[3]) ) return kFALSE;
	236
	237	Double_t mDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
	238	Double_t invmDstar = InvMassDstarKpipi();
	239	if(TMath::Abs(mDstar-invmDstar)>cutsDstar[0]) return kFALSE;
	240
	241	Double_t mD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
	242	if(TMath::Abs((mDstar-mD0)-DeltaInvMass())>cutsDstar[1]) return kFALSE;
	243
3f6adfc9	244	Double_t theta = AngleD0dkpPisoft();
05d80dd6	245	if(theta>cutsDstar[4]) return kFALSE;
05d80dd6	246
	247	return kTRUE;
	248	}
a07ad8e0	249	//-----------------------------------------------------------------------------
a07ad8e0	250	Bool_t AliAODRecoCascadeHF::SelectLctoV0(const Double_t *cutsLctoV0,
9ff3fa8f	251	Bool_t okLck0sp, Bool_t okLcLpi, Bool_t okLcLbarpi) const
a07ad8e0	252	{
	253	// cuts on Lambdac candidates to V0+bachelor
	254	// (to be passed to AliAODRecoDecayHF3Prong::SelectLctoV0())
	255	// 0 = inv. mass half width in K0s hypothesis [GeV]
	256	// 1 = inv. mass half width in Lambda hypothesis [GeV]
	257	// 2 = inv. mass V0 in K0s hypothesis half width [GeV]
	258	// 3 = inv. mass V0 in Lambda hypothesis half width [GeV]
	259	// 4 = pT min Bachelor track [GeV/c]
	260	// 5 = pT min V0-Positive track [GeV/c]
	261	// 6 = pT min V0-Negative track [GeV/c]
9ff3fa8f	262	// 7 = dca cut on the cascade (cm)
9ff3fa8f	263	// 8 = dca cut on the V0 (cm)
a07ad8e0	264
9ff3fa8f	265	// if ( !Getv0() \|\| !Getv0PositiveTrack() \|\| !Getv0NegativeTrack() )
9ff3fa8f	266	// { AliInfo(Form("Not adapted for ESDv0s, return true...")); return false; }
a07ad8e0	267
a07ad8e0	268	Double_t mLck0sp,mLcLpi;
9ff3fa8f	269	okLck0sp=1; okLcLpi=1; okLcLbarpi=1;
a07ad8e0	270
	271	Double_t mLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
	272	Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
	273	Double_t mLPDG = TDatabasePDG::Instance()->GetParticle(3122)->Mass();
	274
	275	// k0s + p
	276	double mk0s = Getv0()->MassK0Short();
	277	mLck0sp = InvMassLctoK0sP();
	278
	279	// lambda + pi
	280	double mlambda = Getv0()->MassLambda();
	281	double malambda = Getv0()->MassAntiLambda();
	282	mLcLpi = InvMassLctoLambdaPi();
	283
	284	// cut on Lc mass
	285	// with k0s p hypothesis
	286	if(TMath::Abs(mLck0sp-mLcPDG)>cutsLctoV0[0]) okLck0sp = 0;
	287	// with Lambda pi hypothesis
	288	if(TMath::Abs(mLcLpi-mLcPDG)>cutsLctoV0[1]) okLcLpi = 0;
9ff3fa8f	289	okLcLbarpi = okLcLpi;
9ff3fa8f	290
a07ad8e0	291	// cuts on the v0 mass
9ff3fa8f	292	if( TMath::Abs(mk0s-mk0sPDG)>cutsLctoV0[2]) okLck0sp = 0;
	293	//if( TMath::Abs(mlambda-mLPDG)>cutsLctoV0[3] &&
	294	//TMath::Abs(malambda-mLPDG)>cutsLctoV0[3] ) okLcLpi = 0;
	295	if( !(GetBachelor()->Charge()==+1 && TMath::Abs(mlambda-mLPDG)<=cutsLctoV0[3]) ) okLcLpi = 0;
	296	if( !(GetBachelor()->Charge()==-1 && TMath::Abs(malambda-mLPDG)<=cutsLctoV0[3]) ) okLcLbarpi = 0;
a07ad8e0	297
9ff3fa8f	298	if(!okLck0sp && !okLcLpi && !okLcLbarpi) return 0;
a07ad8e0	299
	300	// cuts on the minimum pt of the tracks
	301	if(TMath::Abs(GetBachelor()->Pt()) < cutsLctoV0[4]) return 0;
	302	if(TMath::Abs(Getv0PositiveTrack()->Pt()) < cutsLctoV0[5]) return 0;
	303	if(TMath::Abs(Getv0NegativeTrack()->Pt()) < cutsLctoV0[6]) return 0;
	304
a07ad8e0	305	// cut on the cascade dca
9ff3fa8f	306	if( TMath::Abs(GetDCA(0))>cutsLctoV0[7] //\|\|
	307	//TMath::Abs(Getv0()->DcaPosToPrimVertex())>cutsLctoV0[7] \|\|
	308	//TMath::Abs(Getv0()->DcaNegToPrimVertex())>cutsLctoV0[7]
	309	) return 0;
	310
	311	// cut on the v0 dca
	312	if(TMath::Abs(Getv0()->DcaV0Daughters()) > cutsLctoV0[8]) return 0;
	313
	314	// cut on V0 cosine of pointing angle wrt PV
	315	if (CosV0PointingAngle() < cutsLctoV0[9]) { // cosine of V0 pointing angle wrt primary vertex
	316	AliDebug(4,Form(" V0 cosine of pointing angle doesn't pass the cut"));
	317	return 0;
	318	}
	319
	320	// cut on bachelor transverse impact parameter wrt PV
	321	if (TMath::Abs(Getd0Prong(0)) > cutsLctoV0[10]) { // bachelor transverse impact parameter wrt PV
	322	AliDebug(4,Form(" bachelor transverse impact parameter doesn't pass the cut"));
	323	return 0;
	324	}
	325
	326	// cut on V0 transverse impact parameter wrt PV
	327	if (TMath::Abs(Getd0Prong(1)) > cutsLctoV0[11]) { // V0 transverse impact parameter wrt PV
	328	AliDebug(4,Form(" V0 transverse impact parameter doesn't pass the cut"));
	329	return 0;
	330	}
	331
	332	// cut on K0S invariant mass veto
	333	if (TMath::Abs(Getv0()->MassK0Short()-mk0sPDG) < cutsLctoV0[12]) { // K0S invariant mass veto
	334	AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
	335	return 0;
	336	}
	337
	338	// cut on Lambda/LambdaBar invariant mass veto
	339	if (TMath::Abs(Getv0()->MassLambda()-mLPDG) < cutsLctoV0[13] \|\|
	340	TMath::Abs(Getv0()->MassAntiLambda()-mLPDG) < cutsLctoV0[13] ) { // Lambda/LambdaBar invariant mass veto
	341	AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
	342	return 0;
	343	}
	344
	345	// cut on gamma invariant mass veto
	346	if (Getv0()->InvMass2Prongs(0,1,11,11) < cutsLctoV0[14]) { // K0S invariant mass veto
	347	AliDebug(4,Form(" veto on gamma invariant mass doesn't pass the cut"));
	348	return 0;
	349	}
	350
	351	// cut on V0 pT min
	352	if (Getv0()->Pt() < cutsLctoV0[15]) { // V0 pT min
	353	AliDebug(4,Form(" V0 track Pt=%2.2e > %2.2e",Getv0()->Pt(),cutsLctoV0[15]));
	354	return 0;
	355	}
a07ad8e0	356
	357	return true;
	358
	359	}
	360	//-----------------------------------------------------------------------------
3f6adfc9	361	Double_t AliAODRecoCascadeHF::AngleD0dkpPisoft() const {
	362	//
	363	// Angle of soft pion to D0 decay plane
	364	//
	365
	366	TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
	367	TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
	368	TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s
	369
	370	TVector3 perp = p3Trk0.Cross(p3Trk1);
	371	Double_t theta = p3Trk2.Angle(perp);
	372	if(theta>(TMath::Pi()-theta)) theta = TMath::Pi() - theta;
	373	theta = TMath::Pi()/2. - theta;
	374
	375	return theta;
	376	}
ccfae9f3	377	//-----------------------------------------------------------------------------
	378	Bool_t AliAODRecoCascadeHF::TrigonometricalCut() const {
	379	//
	380	// Trigonometrical constraint
	381	//
	382	TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
	383	TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
	384	TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s
	385
	386	Double_t alpha = p3Trk0.Angle(p3Trk2);
	387	Double_t beta = p3Trk1.Angle(p3Trk2);
	388
	389	Double_t cosphi01 = TMath::Cos(alpha) / TMath::Cos(AngleD0dkpPisoft());
	390	Double_t cosphi02 = TMath::Cos(beta) / TMath::Cos(AngleD0dkpPisoft());
	391
	392	Double_t phi01 = TMath::ACos(cosphi01);
	393	Double_t phi02 = TMath::ACos(cosphi02);
	394	Double_t phi00 = p3Trk0.Angle(p3Trk1);
	395
	396	if((phi01>phi00) \|\| (phi02>phi00)) return kFALSE;
	397	return kTRUE;
	398	}
8cb66396	399
c9283e92	400	//-----------------------------------------------------------------------------
	401	Double_t AliAODRecoCascadeHF::DecayLengthV0() const
	402	{
	403	//
	404	// Returns V0 decay length wrt primary vertex
	405	//
	406
	407	AliAODv0 v0 = (AliAODv0)Getv0();
	408
	409	if (!v0)
	410	return -1.;
	411	AliAODVertex *vtxPrimary = GetPrimaryVtx();
	412	Double_t posVtx[3] = {0.,0.,0.};
	413	vtxPrimary->GetXYZ(posVtx);
	414	return v0->DecayLengthV0(posVtx);
	415
	416	}
9ff3fa8f	417	//-----------------------------------------------------------------------------
	418	Double_t AliAODRecoCascadeHF::DecayLengthXYV0() const
	419	{
	420	//
	421	// Returns transverse V0 decay length wrt primary vertex
	422	//
	423	AliAODv0 v0 = (AliAODv0)Getv0();
	424
	425	if (!v0)
	426	return -1.;
	427	AliAODVertex *vtxPrimary = GetPrimaryVtx();
	428	Double_t posVtx[3] = {0.,0.,0.};
	429	vtxPrimary->GetXYZ(posVtx);
	430	return v0->DecayLengthXY(posVtx);
	431
	432	}
	433	//-----------------------------------------------------------------------------
	434	Double_t AliAODRecoCascadeHF::CosV0PointingAngle() const
	435	{
	436	//
	437	// Returns cosine of V0 pointing angle wrt primary vertex
	438	//
	439
	440	AliAODv0 v0 = (AliAODv0)Getv0();
	441
	442	if (!v0)
	443	return -999.;
	444
	445	AliAODVertex *vtxPrimary = GetPrimaryVtx();
	446	Double_t posVtx[3] = {0.,0.,0.};
	447	vtxPrimary->GetXYZ(posVtx);
	448	return v0->CosPointingAngle(posVtx);
	449
	450	}
ff12b981	451	//-----------------------------------------------------------------------------
	452	Double_t AliAODRecoCascadeHF::CosV0PointingAngleXY() const
	453	{
	454	//
	455	// Returns XY cosine of V0 pointing angle wrt primary vertex
	456	//
	457
	458	AliAODv0 v0 = (AliAODv0)Getv0();
	459
	460	if (!v0)
	461	return -999.;
	462
	463	AliAODVertex *vtxPrimary = GetPrimaryVtx();
	464	Double_t posVtx[3] = {0.,0.,0.};
	465	vtxPrimary->GetXYZ(posVtx);
	466	return v0->CosPointingAngleXY(posVtx);
	467
	468	}
	469	//-----------------------------------------------------------------------------
	470	Double_t AliAODRecoCascadeHF::NormalizedV0DecayLength() const
	471	{
	472	//
	473	// Returns V0 normalized decay length wrt primary vertex
	474	//
	475
	476	AliAODv0 v0 = (AliAODv0)Getv0();
	477
	478	if (!v0)
	479	return -1.;
	480	//AliAODVertex *vtxPrimary = GetPrimaryVtx();
	481	//Double_t posVtx[3] = {0.,0.,0.};
	482	//vtxPrimary->GetXYZ(posVtx);
	483	//return v0->NormalizedDecayLength(posVtx);
	484	return v0->NormalizedDecayLength(GetPrimaryVtx());
	485
	486	}
	487	//-----------------------------------------------------------------------------
	488	Double_t AliAODRecoCascadeHF::NormalizedV0DecayLengthXY() const
	489	{
	490	//
	491	// Returns transverse V0 normalized decay length wrt primary vertex
	492	//
	493	AliAODv0 v0 = (AliAODv0)Getv0();
	494
	495	if (!v0)
	496	return -1.;
	497	//AliAODVertex *vtxPrimary = GetPrimaryVtx();
	498	//Double_t posVtx[3] = {0.,0.,0.};
	499	//vtxPrimary->GetXYZ(posVtx);
	500	//return v0->NormalizedDecayLengthXY(posVtx);
	501	return v0->NormalizedDecayLengthXY(GetPrimaryVtx());
	502
	503	}