]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
D0->Kpipipi analysis included (Rossella)
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Feb 2009 00:35:02 +0000 (00:35 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Feb 2009 00:35:02 +0000 (00:35 +0000)
PWG3/vertexingHF/AliAODRecoDecayHF4Prong.cxx
PWG3/vertexingHF/AliAODRecoDecayHF4Prong.h
PWG3/vertexingHF/AliAnalysisVertexingHF.cxx
PWG3/vertexingHF/AliAnalysisVertexingHF.h

index 2185a8cc28404bfb91eac62bf8cbe27fcc6b3760..da325a68580b4ca5036b0c51047bc47367e68598 100644 (file)
@@ -30,11 +30,9 @@ ClassImp(AliAODRecoDecayHF4Prong)
 //--------------------------------------------------------------------------
 AliAODRecoDecayHF4Prong::AliAODRecoDecayHF4Prong() :
   AliAODRecoDecayHF(), 
-  //fSigmaVert(0),
   fDist12toPrim(0),
-  fDist23toPrim(0),
-  fDist14toPrim(0),
-  fDist34toPrim(0)
+  fDist3toPrim(0),
+  fDist4toPrim(0)
 {
   //
   // Default Constructor
@@ -45,15 +43,14 @@ AliAODRecoDecayHF4Prong::AliAODRecoDecayHF4Prong(AliAODVertex *vtx2,
                                                 Double_t *px,Double_t *py,Double_t *pz,
                                                 Double_t *d0,Double_t *d0err,
                                                 Double_t *dca, //Double_t sigvert,
-                                                Double_t dist12,Double_t dist23,
-                                                Double_t dist14,Double_t dist34,
+                                                Double_t dist12,Double_t dist3,
+                                                Double_t dist4,
                                                 Short_t charge) :
   AliAODRecoDecayHF(vtx2,4,charge,px,py,pz,d0,d0err),
   // fSigmaVert(sigvert),
   fDist12toPrim(dist12),
-  fDist23toPrim(dist23),
-  fDist14toPrim(dist14),
-  fDist34toPrim(dist34)
+  fDist3toPrim(dist3),
+  fDist4toPrim(dist4)
 {
   //
   // Constructor with AliAODVertex for decay vertex
@@ -64,15 +61,14 @@ AliAODRecoDecayHF4Prong::AliAODRecoDecayHF4Prong(AliAODVertex *vtx2,
 AliAODRecoDecayHF4Prong::AliAODRecoDecayHF4Prong(AliAODVertex *vtx2,
                                                 Double_t *d0,Double_t *d0err,
                                                 Double_t *dca, //Double_t sigvert,
-                                                Double_t dist12,Double_t dist23, 
-                                                Double_t dist14,Double_t dist34, 
+                                                Double_t dist12,Double_t dist3, 
+                                                Double_t dist4, 
                                                 Short_t charge) :
   AliAODRecoDecayHF(vtx2,4,charge,d0,d0err),
   //fSigmaVert(sigvert),
   fDist12toPrim(dist12),
-  fDist23toPrim(dist23),
-  fDist14toPrim(dist14),
-  fDist34toPrim(dist34)
+  fDist3toPrim(dist3),
+  fDist4toPrim(dist4)
 {
   //
   // Constructor with AliAODVertex for decay vertex and without prongs momenta
@@ -84,9 +80,8 @@ AliAODRecoDecayHF4Prong::AliAODRecoDecayHF4Prong(const AliAODRecoDecayHF4Prong &
   AliAODRecoDecayHF(source),
   //fSigmaVert(source.fSigmaVert),
   fDist12toPrim(source.fDist12toPrim),
-  fDist23toPrim(source.fDist23toPrim),
-  fDist14toPrim(source.fDist14toPrim),
-  fDist34toPrim(source.fDist34toPrim)
+  fDist3toPrim(source.fDist3toPrim),
+  fDist4toPrim(source.fDist4toPrim)
 {
   //
   // Copy constructor
@@ -103,49 +98,168 @@ AliAODRecoDecayHF4Prong &AliAODRecoDecayHF4Prong::operator=(const AliAODRecoDeca
   AliAODRecoDecayHF::operator=(source);
 
   fDist12toPrim= source.fDist12toPrim;
-  fDist23toPrim= source.fDist23toPrim;
-  fDist14toPrim= source.fDist14toPrim;
-  fDist34toPrim= source.fDist34toPrim;
+  fDist3toPrim= source.fDist3toPrim;
+  fDist4toPrim= source.fDist4toPrim;
   //fSigmaVert= source.fSigmaVert;
 
   return *this;
 }
 //--------------------------------------------------------------------------
-Bool_t AliAODRecoDecayHF4Prong::SelectD0(const Double_t *cuts,Int_t &okD0,Int_t &okD0bar)
-  const {
+void AliAODRecoDecayHF4Prong::InvMassD0(Double_t mD0[2]) const {
+  //
+  // the masses for D0 hypothesis
+  // 
+  UInt_t pdg[4];
+  pdg[0]=211; pdg[1]=321; pdg[2]=211; pdg[3]=211;
+  mD0[0]=InvMass(4,pdg);
+  pdg[1]=211; pdg[3]=321;
+  mD0[1]=InvMass(4,pdg);
+
+  return;
+}
+//--------------------------------------------------------------------------
+void AliAODRecoDecayHF4Prong::InvMassD0bar(Double_t mD0bar[2]) const {
+  //
+  // the two masses for D0bar hypothesis
+  //
+  UInt_t pdg[4];
+  pdg[0]=321; pdg[1]=211; pdg[2]=211; pdg[3]=211;
+  mD0bar[0]=InvMass(4,pdg);
+  pdg[0]=211; pdg[2]=321;
+  mD0bar[1]=InvMass(4,pdg);
+  
+  return;
+}
+//--------------------------------------------------------------------------
+Bool_t AliAODRecoDecayHF4Prong::SelectD0(const Double_t *cuts,Int_t &okD0,Int_t &okD0bar) const {
 //
 // This function compares the D0 with a set of cuts:
 // 
-// to be implemented 
-//
-// cuts[1] = ... 
-// cuts[2] = ...
-// cuts[3] = ...
-// cuts[4] = ...
+// cuts[0] = invariant mass cut 
+// cuts[1] = DCA between opposite sign tracks 
+// cuts[2] = Distance between primary and two tracks vertex fDist12toPrim
+// cuts[3] = Distance between primary and three tracks vertex fDist3toPrim
+// cuts[4] = Distance between primary and two tracks vertex fDist4toPrim
+// cuts[5] = Cosinus of the pointing angle
+// cuts[6] = Transverse momentum of the D0 candidate
+// cuts[7] = Mass Pi+Pi- = mass of the rho0
+// cuts[8] = PID cut (one K in the quadruplet)
 //
 // If candidate D0 does not pass the cuts return kFALSE
 //
 
   okD0=0; okD0bar=0;
-  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
-  Double_t mD0=InvMassD0();
-  if(TMath::Abs(mD0-mD0PDG)>cuts[0]) return kFALSE;
+//  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
+//  Double_t mD0=InvMassD0();
+//  if(TMath::Abs(mD0-mD0PDG)>cuts[0]) return kFALSE;
 
   //single track
   //if(TMath::Abs(PtProng(2)) < cuts[2] || TMath::Abs(Getd0Prong(2))<cuts[4])return kFALSE;//Pion2
-
-  //DCA
-  //for(Int_t i=0;i<3;i++) if(cuts[11]>0 && GetDCA(i)>cuts[11])return kFALSE;
-
+  
+  //DCA btw opposite sign tracks
+  if(cuts[1]>0.){
+    if(GetDCA(0)>cuts[1]) return kFALSE;
+    if(GetDCA(1)>cuts[1]) return kFALSE;
+    if(GetDCA(2)>cuts[1]) return kFALSE;
+    if(GetDCA(3)>cuts[1]) return kFALSE;
+    if(GetDCA(4)>cuts[1]) return kFALSE;
+    if(GetDCA(5)>cuts[1]) return kFALSE;
+  }
   //2track cuts
-  //if(fDist12toPrim<cuts[5] || fDist23toPrim<cuts[5])return kFALSE;
-  //if(Getd0Prong(0)*Getd0Prong(1)<0 && Getd0Prong(2)*Getd0Prong(1)<0)return kFALSE;
-
-  //sec vert
-  //if(fSigmaVert>cuts[6])return kFALSE;
+  if(cuts[2]>0.){
+    if(fDist12toPrim>10.)return kFALSE;
+    if(fDist12toPrim<cuts[2])return kFALSE;
+  }
+  //3track cuts
+  if(cuts[3]>0.){
+    if(fDist3toPrim<cuts[3])return kFALSE;
+  }
+  //4track cuts
+  if(cuts[4]>0.){
+    if(fDist4toPrim<cuts[4])return kFALSE;
+  }
+  if(cuts[5]>-1.1){
+    if(CosPointingAngle()<cuts[5])return kFALSE;
+  }
+  if(cuts[6]>0.){
+    if(Pt()<cuts[6])return kFALSE;
+  }
+  if(cuts[7]>0.){
+    Double_t massD0[2];
+    Double_t massD0bar[2];
+    Bool_t good=CutRhoMass(massD0,massD0bar,cuts[0],cuts[7]);
+    if(!good) return kFALSE;
+  }
 
-  //if(DecayLength()<cuts[7])return kFALSE;
-  
-  return kFALSE; // for the time being 
-  // return kTRUE; // after implementation 
+  return kTRUE;
+}
+//----------------------------------------------------------------------------
+Bool_t AliAODRecoDecayHF4Prong::CutRhoMass(Double_t massD0[2],Double_t massD0bar[2],Double_t cutMass,Double_t cutRho) const {
+  //
+  // cut on the mass of rho->pipi
+  //
+  Bool_t isGood=kFALSE;
+  Int_t nprongs=4;
+  for(Int_t i=0;i<2;i++){massD0[i]=0.;massD0bar[i]=0.;}
+  Bool_t isRho=kFALSE;
+  Bool_t isTrue=kFALSE;
+  Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
+  Double_t mPDGRho=TDatabasePDG::Instance()->GetParticle(113)->Mass();
+  Double_t minv01=InvMassRho(0,1);
+  if(TMath::Abs(minv01-mPDGRho)<cutRho) isRho=kTRUE;
+  if(isRho){
+    UInt_t pdg1[4]={211,211,321,211};
+    Double_t mass1=InvMass(nprongs,pdg1);
+    if(TMath::Abs(mass1-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
+    if(isTrue) massD0bar[1]=mass1;
+    isTrue=kFALSE;
+    UInt_t pdg2[4]={211,211,211,321};
+    Double_t mass2=InvMass(4,pdg2);
+    if(TMath::Abs(mass2-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
+    if(isTrue) massD0[1]=mass2;
+    isTrue=kFALSE;
+  }
+  Double_t minv03=InvMassRho(0,3);
+  if(TMath::Abs(minv03-mPDGRho)<cutRho) isRho=kTRUE;
+  if(isRho){
+    UInt_t pdg1[4]={211,211,321,211};
+    Double_t mass1=InvMass(4,pdg1);
+    if(TMath::Abs(mass1-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
+    if(isTrue) massD0bar[1]=mass1;
+    isTrue=kFALSE;
+    UInt_t pdg2[4]={211,321,211,211};
+    Double_t mass2=InvMass(4,pdg2);
+    if(TMath::Abs(mass2-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
+    if(isTrue) massD0[0]=mass2;
+    isTrue=kFALSE;
+  }
+  Double_t minv12=InvMassRho(1,2);
+  if(TMath::Abs(minv12-mPDGRho)<cutRho) isRho=kTRUE;
+  if(isRho){
+    UInt_t pdg1[4]={321,211,211,211};
+    Double_t mass1=InvMass(4,pdg1);
+    if(TMath::Abs(mass1-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
+    if(isTrue) massD0bar[0]=mass1;
+    isTrue=kFALSE;
+    UInt_t pdg2[4]={211,211,211,321};
+    Double_t mass2=InvMass(4,pdg2);
+    if(TMath::Abs(mass2-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
+    if(isTrue) massD0[1]=mass2;
+    isTrue=kFALSE;
+  }
+  Double_t minv23=InvMassRho(2,3);
+  if(TMath::Abs(minv23-mPDGRho)<cutRho) isRho=kTRUE;
+  if(isRho){
+    UInt_t pdg1[4]={321,211,211,211};
+    Double_t mass1=InvMass(4,pdg1);
+    if(TMath::Abs(mass1-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
+    if(isTrue) massD0bar[0]=mass1;
+    isTrue=kFALSE;
+    UInt_t pdg2[4]={211,321,211,211};
+    Double_t mass2=InvMass(4,pdg2);
+    if(TMath::Abs(mass2-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
+    if(isTrue) massD0[0]=mass2;
+    isTrue=kFALSE;
+  }
+  return isGood;
 }
index cdee57a3e3afbaaca8ebc9c08e4c13274210cd55..8fdc9a6f55588bbd65fb983b70b9114c34329897 100644 (file)
@@ -22,14 +22,14 @@ class AliAODRecoDecayHF4Prong : public AliAODRecoDecayHF {
                           Double_t *px,Double_t *py,Double_t *pz,
                           Double_t *d0,Double_t *d0err,
                           Double_t *dca, //Double_t sigvert,
-                          Double_t dist12,Double_t dist23, 
-                          Double_t dist14,Double_t dist34, 
+                          Double_t dist12,Double_t dist3, 
+                          Double_t dist4, 
                            Short_t charge);
    AliAODRecoDecayHF4Prong(AliAODVertex *vtx2,
                           Double_t *d0,Double_t *d0err,
                           Double_t *dca, //Double_t sigvert,
-                          Double_t dist12,Double_t dist23, 
-                          Double_t dist14,Double_t dist34, 
+                          Double_t dist12,Double_t dist3, 
+                          Double_t dist4, 
                           Short_t charge);
 
   AliAODRecoDecayHF4Prong(const AliAODRecoDecayHF4Prong& source);
@@ -38,11 +38,11 @@ class AliAODRecoDecayHF4Prong : public AliAODRecoDecayHF {
   virtual ~AliAODRecoDecayHF4Prong() {}  
  
   void GetDCAs(Double_t dca[6]) const 
-    {for(Int_t i=0;i<6;i++) dca[i]=GetDCA(i);} // convention:fDCA[0]=p0p1,fDCA[1]=p0p2,fDCA[2]=p0p3,fDCA[3]=p1p2...
+    {for(Int_t i=0;i<6;i++) dca[i]=GetDCA(i);} 
+  // convention:fDCA[0]=p0p1,fDCA[1]=p0p2,fDCA[2]=p0p3,fDCA[3]=p1p2...
   Double_t GetDist12toPrim() const {return fDist12toPrim;}
-  Double_t GetDist23toPrim() const {return fDist12toPrim;}
-  Double_t GetDist14toPrim() const {return fDist12toPrim;}
-  Double_t GetDist34toPrim() const {return fDist12toPrim;}
+  Double_t GetDist3toPrim() const {return fDist3toPrim;}
+  Double_t GetDist4toPrim() const {return fDist4toPrim;}
 
   // D0->pi+K- pipi and D0bar->K+pi- pipi (in the order given) 
   Double_t ED0() const {return E(421);}
@@ -50,26 +50,24 @@ class AliAODRecoDecayHF4Prong : public AliAODRecoDecayHF {
   Double_t CtD0() const {return Ct(421);} 
   Double_t CtD0(Double_t point[3]) const {return AliAODRecoDecay::Ct(421,point);}
   Double_t CtD0(AliAODVertex *vtx1) const {return AliAODRecoDecay::Ct(421,vtx1);}
-  Double_t InvMassRho() const {return InvMass2Prongs(2,3,211,211);} 
-  Double_t InvMassD0() const {UInt_t pdg[4]={211,321,211,211};return InvMass(4,pdg);}
-  Double_t InvMassD0bar() const {UInt_t pdg[4]={321,211,211,211};return InvMass(4,pdg);}
-  void InvMassD0(Double_t &mD0,Double_t &mD0bar) const
-    {mD0=InvMassD0();mD0bar=InvMassD0bar();return;}
+  Double_t InvMassRho(Int_t i,Int_t j) const {return InvMass2Prongs(i,j,211,211);} 
+  Bool_t CutRhoMass(Double_t massD0[2],Double_t massD0bar[2],Double_t CutMass,Double_t CutRho) const; 
 
+  void InvMassD0(Double_t mD0[2]) const;
+  void InvMassD0bar(Double_t mD0bar[2]) const;
 
   Bool_t   SelectD0(const Double_t* cuts,Int_t &okD0,Int_t &okD0bar) const;
 
-
  private:
 
   //Double_t fSigmaVert; // track dispersion around the secondary vertex
   Double_t fDist12toPrim; //distance prim vert - 2 opposite sign track vertex 
-  Double_t fDist23toPrim; //distance prim vert - 2 opposite sign track vertex 
-  Double_t fDist14toPrim; //distance prim vert - 2 opposite sign track vertex 
-  Double_t fDist34toPrim; //distance prim vert - 2 opposite sign track vertex 
+  Double_t fDist3toPrim; //distance prim vert - 3 track vertex 
+  Double_t fDist4toPrim; //distance prim vert - 4 track vertex 
   //Double_t fDist123toPrim;  //distance 
 
-  ClassDef(AliAODRecoDecayHF4Prong,1)  // base class for AOD reconstructed 
+  ClassDef(AliAODRecoDecayHF4Prong,2)  // base class for AOD reconstructed 
                                        // heavy-flavour 3-prong decays
 };
 
index 16058c6f7ed0f10f11d0cac39d09dadef9aeb8b0..b03a636b9e0451389638fad538bb132de7437ce5 100644 (file)
@@ -77,7 +77,7 @@ fTrackFilterSoftPi(0x0)
   SetDsCuts();
   SetLcCuts();
   SetDstarCuts();
-
+  SetD0to4ProngsCuts();
 }
 //--------------------------------------------------------------------------
 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) : 
@@ -107,6 +107,7 @@ fTrackFilterSoftPi(source.fTrackFilterSoftPi)
   for(Int_t i=0; i<13; i++) fDsCuts[i]=source.fDsCuts[i];
   for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
   for(Int_t i=0; i<5; i++)  fDstarCuts[i]=source.fDstarCuts[i];
+  for(Int_t i=0; i<9; i++)  fD0to4ProngsCuts[i]=source.fD0to4ProngsCuts[i];
 }
 //--------------------------------------------------------------------------
 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
@@ -136,6 +137,7 @@ AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVerte
   for(Int_t i=0; i<13; i++) fDsCuts[i]=source.fDsCuts[i];
   for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
   for(Int_t i=0; i<5; i++)  fDstarCuts[i]=source.fDstarCuts[i];
+  for(Int_t i=0; i<9; i++)  fD0to4ProngsCuts[i]=source.fD0to4ProngsCuts[i];
 
   return *this;
 }
@@ -236,7 +238,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
   AliAODRecoCascadeHF     *ioCascade = 0;
 
   Int_t    iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries;
-  Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcaCasc;
+  Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaCasc;
   Bool_t   okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
   Bool_t   okDstar=kFALSE,okD0fromDstar=kFALSE;
   AliESDtrack *postrack1 = 0;
@@ -247,6 +249,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
   Double_t dcaMax = fD0toKpiCuts[1];
   if(dcaMax < fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1];
   if(dcaMax < fDplusCuts[11])  dcaMax=fDplusCuts[11];
+  if(dcaMax < fD0to4ProngsCuts[1])  dcaMax=fD0to4ProngsCuts[1];
   AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
 
 
@@ -516,6 +519,16 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
        // 4 prong candidates
        if(f4Prong) {
 
+         // back to primary vertex
+         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+         postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
+         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+         // Vertexing for these 3 (can be taken from above?)
+          threeTrackArray->AddAt(postrack1,0);
+          threeTrackArray->AddAt(negtrack1,1);
+          threeTrackArray->AddAt(postrack2,2);
+          AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
+
          // 3rd LOOP  ON  NEGATIVE  TRACKS (for 4 prong) 
          for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
 
@@ -537,6 +550,8 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
            negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
            dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
            if(dcap1n2>dcaMax) { negtrack2=0; continue; }
+            dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
+            if(dcap2n2>dcaMax) { negtrack2=0; continue; }
 
            // Vertexing
            fourTrackArray->AddAt(postrack1,0);
@@ -545,7 +560,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
            fourTrackArray->AddAt(negtrack2,3);
 
            AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
-           io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp2n1,dcap1n1,dcap1n2,dcap2n1,ok4Prong);
+           io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
            if(ok4Prong) {
              AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
              if(fInputAOD) AddDaughterRefs(v4Prong,event,fourTrackArray);
@@ -560,6 +575,9 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
 
          } // end loop on negative tracks
 
+          threeTrackArray->Clear();
+         delete vertexp1n1p2;
+
        }
 
        postrack2 = 0;
@@ -944,10 +962,12 @@ AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
 //----------------------------------------------------------------------------
 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
                              TObjArray *fourTrackArray,AliVEvent *event,
-                            AliAODVertex *secVert,
-                            AliAODVertex *vertexp1n1,AliAODVertex *vertexp2n1,
-                            Double_t dcap1n1,Double_t dcap1n2,Double_t dcap2n1,
-                            Bool_t &ok4Prong) const
+                             AliAODVertex *secVert,
+                             AliAODVertex *vertexp1n1,
+                             AliAODVertex *vertexp1n1p2,
+                             Double_t dcap1n1,Double_t dcap1n2,
+                            Double_t dcap2n1,Double_t dcap2n2,
+                             Bool_t &ok4Prong) const
 {
   // Make 4Prong candidates and check if they pass D0toKpipipi
   // reconstruction cuts
@@ -958,8 +978,6 @@ AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
 
   Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
 
-  px[0]=dcap1n1*dcap1n2*dcap2n1; // TO BE CHANGED (done just to removed compilation warning about dca... not used)
-
   AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
   AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
   AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
@@ -986,51 +1004,66 @@ AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
   //if(!okMassCut) {
   //  if(fDebug) printf(" candidate didn't pass mass cut\n");
   //  return 0x0;
-  //}
+  // }
 
   // primary vertex to be used by this candidate
   AliAODVertex *primVertexAOD  = PrimaryVertex(fourTrackArray,event);
   if(!primVertexAOD) return 0x0;
 
-  /*
-    Double_t d0z0[2],covd0z0[3];
-    postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
-    negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
-    postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
-    negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
-  */
+  Double_t d0z0[2],covd0z0[3];
+  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
+  d0[0]=d0z0[0];
+  d0err[0] = TMath::Sqrt(covd0z0[0]);
+  negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
+  d0[1]=d0z0[0];
+  d0err[1] = TMath::Sqrt(covd0z0[0]);
+  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
+  d0[2]=d0z0[0];
+  d0err[2] = TMath::Sqrt(covd0z0[0]);
+  negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
+  d0[3]=d0z0[0];
+  d0err[3] = TMath::Sqrt(covd0z0[0]);
+
 
   // create the object AliAODRecoDecayHF4Prong
   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
-  Double_t dca[6]={0.,0.,0.,0.,0.,0.}; //  modify it
+  Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
   Double_t dist12=TMath::Sqrt((vertexp1n1->GetX()-pos[0])*(vertexp1n1->GetX()-pos[0])+(vertexp1n1->GetY()-pos[1])*(vertexp1n1->GetY()-pos[1])+(vertexp1n1->GetZ()-pos[2])*(vertexp1n1->GetZ()-pos[2]));
-  Double_t dist23=TMath::Sqrt((vertexp2n1->GetX()-pos[0])*(vertexp2n1->GetX()-pos[0])+(vertexp2n1->GetY()-pos[1])*(vertexp2n1->GetY()-pos[1])+(vertexp2n1->GetZ()-pos[2])*(vertexp2n1->GetZ()-pos[2]));
-  Double_t dist14=0.; // to be implemented
-  Double_t dist34=0.; // to be implemented
+  Double_t dist3=TMath::Sqrt((vertexp1n1p2->GetX()-pos[0])*(vertexp1n1p2->GetX()-pos[0])+(vertexp1n1p2->GetY()-pos[1])*(vertexp1n1p2->GetY()-pos[1])+(vertexp1n1p2->GetZ()-pos[2])*(vertexp1n1p2->GetZ()-pos[2]));
+  Double_t dist4=TMath::Sqrt((secVert->GetX()-pos[0])*(secVert->GetX()-pos[0])+(secVert->GetY()-pos[1])*(secVert->GetY()-pos[1])+(secVert->GetZ()-pos[2])*(secVert->GetZ()-pos[2]));
   Short_t charge=0;
-  AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist23,dist14,dist34,charge);
+  AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
   the4Prong->SetOwnPrimaryVtx(primVertexAOD);
   UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
   the4Prong->SetProngIDs(4,id);
 
+  // invariant mass cut for rho or D0 (try to improve coding here..)
+  Bool_t okMassCut=kFALSE;
+  if(!okMassCut && fD0to4ProngsCuts[8]==0.){      //no PID, to be implemented with PID
+    if(SelectInvMass(4,4,px,py,pz)) okMassCut=kTRUE;
+  }
+  if(!okMassCut) {
+    //if(fDebug) printf(" candidate didn't pass mass cut\n");
+    //printf(" candidate didn't pass mass cut\n");
+    return 0x0;
+  }
+
 
-  // use the following two lines once AliAODRecoDecayHF4Prong::SelectD0 is available
-  // select D0->Kpipipi
-  //Int_t checkD0,checkD0bar;   
-  // ok4Prong=the4Prong->SelectD0(fD04pCuts,checkD0,checkD0bar); 
-  ok4Prong=kFALSE;  //for the time being ...
+  Int_t checkD0,checkD0bar;
+  ok4Prong=the4Prong->SelectD0(fD0to4ProngsCuts,checkD0,checkD0bar);
 
   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) the4Prong->UnsetOwnPrimaryVtx();
 
   // get PID info from ESD
-  Double_t esdpid0[5];
-  postrack1->GetESDpid(esdpid0);
-  Double_t esdpid1[5];
-  negtrack1->GetESDpid(esdpid1);
-  Double_t esdpid2[5];
-  postrack2->GetESDpid(esdpid2);
-  Double_t esdpid3[5];
-  negtrack2->GetESDpid(esdpid3);
+  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
+  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
+  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
+  if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
+  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
+  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
+  Double_t esdpid3[5]={0.,0.,0.,0.,0.};
+  if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
 
   Double_t esdpid[20];
   for(Int_t i=0;i<5;i++) {
@@ -1040,7 +1073,7 @@ AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
     esdpid[15+i] = esdpid3[i];
   }
   the4Prong->SetPID(4,esdpid);
-
+  
   return the4Prong;
 }
 //-----------------------------------------------------------------------------
@@ -1307,7 +1340,7 @@ Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
   AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
   delete [] dummyd0;
 
-  UInt_t pdg2[2],pdg3[3];
+  UInt_t pdg2[2],pdg3[3],pdg4[4];
   Double_t mPDG,minv;
 
   Bool_t retval=kFALSE;
@@ -1356,6 +1389,25 @@ Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
       minv = rd->InvMass(nprongs,pdg2);
       if(TMath::Abs(minv-mPDG)<fDstarCuts[0]) retval=kTRUE;
       break;
+    case 4:                 // D0->K3pi without PID
+      pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
+      mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
+      minv = rd->InvMass(nprongs,pdg4);
+      if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0])
+       retval=kTRUE;
+      pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
+      mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
+      minv = rd->InvMass(nprongs,pdg4);
+      if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
+      pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
+      mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
+      minv = rd->InvMass(nprongs,pdg4);
+      if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
+      pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
+      mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
+      minv = rd->InvMass(nprongs,pdg4);
+      if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
+      break;
     default:
       printf("SelectInvMass(): wrong decay selection\n");
       break;
@@ -1659,6 +1711,35 @@ void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[12])
   return;
 }
 //-----------------------------------------------------------------------------
+void AliAnalysisVertexingHF::SetD0to4ProngsCuts(Double_t cut0,Double_t cut1,
+                                   Double_t cut2,Double_t cut3,Double_t cut4,
+                                   Double_t cut5,Double_t cut6,
+                                   Double_t cut7,Double_t cut8)
+{
+  // Set the cuts for D0->Kpipipi selection
+
+  fD0to4ProngsCuts[0] = cut0;
+  fD0to4ProngsCuts[1] = cut1;
+  fD0to4ProngsCuts[2] = cut2;
+  fD0to4ProngsCuts[3] = cut3;
+  fD0to4ProngsCuts[4] = cut4;
+  fD0to4ProngsCuts[5] = cut5;
+  fD0to4ProngsCuts[6] = cut6;
+  fD0to4ProngsCuts[7] = cut7;
+  fD0to4ProngsCuts[8] = cut8;
+
+  return;
+}
+//-----------------------------------------------------------------------------
+void AliAnalysisVertexingHF::SetD0to4ProngsCuts(const Double_t cuts[9])
+{
+  // Set the cuts for D0->Kpipipi selection
+
+  for(Int_t i=0; i<9; i++) fD0to4ProngsCuts[i] = cuts[i];
+
+  return;
+}
+//-----------------------------------------------------------------------------
 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
                                             Bool_t &okDisplaced,Bool_t &okSoftPi) const 
 {
index b5c2b64e46c5bff1251d46a58b64f438d38559f1..aded400859d82573169b9b2e134fe0300193c460 100644 (file)
@@ -110,6 +110,11 @@ class AliAnalysisVertexingHF : public TNamed {
                       Double_t cut2=-1., Double_t cut3=1000.,
                       Double_t cut4=1.6); 
   void SetDstarCuts(const Double_t cuts[5]); 
+  void SetD0to4ProngsCuts(Double_t cut0=1000.,Double_t cut1=100000.,
+                      Double_t cut2=0.,Double_t cut3=0.,Double_t cut4=0.,
+                      Double_t cut5=-1.1,Double_t cut6=0.,
+                      Double_t cut7=0.,Double_t cut8=0.);
+  void SetD0to4ProngsCuts(const Double_t cuts[9]);
   const Double_t *GetD0toKpiCuts() const {return fD0toKpiCuts;}
   const Double_t *GetD0fromDstarCuts() const {return fD0fromDstarCuts;}
   const Double_t *GetBtoJPSICuts() const {return fBtoJPSICuts;}
@@ -117,6 +122,7 @@ class AliAnalysisVertexingHF : public TNamed {
   const Double_t *GetDsCuts() const {return fDsCuts;}
   const Double_t *GetLcCuts() const {return fLcCuts;}
   const Double_t *GetDstarCuts() const {return fDstarCuts;}
+  const Double_t *GetD0to4ProngsCuts() const {return fD0to4ProngsCuts;}
 
   //
  private:
@@ -234,6 +240,18 @@ class AliAnalysisVertexingHF : public TNamed {
                         // 3 = PtMax of pi_s [GeV/c]
                         // 4 = theta, angle between the track of pi_s and D0 decay plane [rad]
 
+  Double_t fD0to4ProngsCuts[9]; // cuts on D0->K3pi candidates
+                        // (to be passed to AliAODRecoDecayHF4Prong::SelectD0())
+                        // 0 = inv. mass half width of D0 [GeV]
+                        // 1 = DCA between opposite sign tracks 
+                        // 2 = Distance between primary and two tracks vertex fDist12toPrim
+                        // 3 = Distance between primary and three tracks vertex fDist3toPrim
+                        // 4 = Distance between primary and two tracks vertex fDist4toPrim
+                        // 5 = Cosinus of the pointing angle
+                        // 6 = Transverse momentum of the D0 candidate
+                        // 7 = Mass Pi+Pi- = mass of the rho0
+                        // 8 = PID cut (one K in the quadruplet)
+
   //
   void AddDaughterRefs(AliAODVertex *v,AliVEvent *event,
                       TObjArray *trkArray) const;
@@ -248,12 +266,12 @@ class AliAnalysisVertexingHF : public TNamed {
                                      Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
                                      Bool_t &ok3Prong) const;
   AliAODRecoDecayHF4Prong* Make4Prong(TObjArray *fourTrackArray,AliVEvent *event,
-                                     AliAODVertex *secVert,
-                                     AliAODVertex *vertexp1n1,
-                                     AliAODVertex *vertexp2n1,
-                                     Double_t dcap1n1,Double_t dcap1n2,
-                                     Double_t dcap2n1,
-                                     Bool_t &ok4Prong) const;
+                                      AliAODVertex *secVert,
+                                      AliAODVertex *vertexp1n1,
+                                      AliAODVertex *vertexp1n1p2,
+                                      Double_t dcap1n1,Double_t dcap1n2,
+                                      Double_t dcap2n1,Double_t dcap2n2,
+                                      Bool_t &ok4Prong) const;
   AliAODRecoCascadeHF* MakeCascade(TObjArray *twoTrackArray,AliVEvent *event,
                                   AliAODVertex *secVert,
                                   AliAODRecoDecayHF2Prong *rd2Prong,
@@ -270,7 +288,7 @@ class AliAnalysisVertexingHF : public TNamed {
   void   SetPrimaryVertex(AliESDVertex *v1) { fV1 = v1; }
   Bool_t SingleTrkCuts(AliESDtrack *trk,Bool_t &okDisplaced,Bool_t &okSoftPi) const;
   //
-  ClassDef(AliAnalysisVertexingHF,7)  // Reconstruction of HF decay candidates
+  ClassDef(AliAnalysisVertexingHF,8);  // Reconstruction of HF decay candidates
 };