PWG3/vertexingHF/AliAnalysisVertexingHF.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, 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
  16 //----------------------------------------------------------------------------
  17 //    Implementation of the heavy-flavour vertexing analysis class
  18 // Candidates are stored in the AOD as objects deriving from AliAODRecoDecay.
  19 // To be used as a task of AliAnalysisManager by means of the interface
  20 // class AliAnalysisTaskSEVertexingHF.
  21 // An example of usage in the macro AliAnalysisTaskSEVertexingHFTest.C.
  22 //
  23 //  Contact: andrea.dainese@lnl.infn.it
  24 //  Contributors: E.Bruna, G.E.Bruno, A.Dainese, C.Di Gliglio,
  25 //                F.Prino, R.Romita, X.M.Zhang
  26 //----------------------------------------------------------------------------
  27 #include <TFile.h>
  28 #include <TDatabasePDG.h>
  29 #include <TString.h>
  30 #include "AliLog.h"
  31 #include "AliVEvent.h"
  32 #include "AliVVertex.h"
  33 #include "AliVTrack.h"
  34 #include "AliVertexerTracks.h"
  35 #include "AliKFVertex.h"
  36 #include "AliESDEvent.h"
  37 #include "AliESDVertex.h"
  38 #include "AliExternalTrackParam.h"
  39 #include "AliNeutralTrackParam.h"
  40 #include "AliESDtrack.h"
  41 #include "AliESDtrackCuts.h"
  42 #include "AliAODEvent.h"
  43 #include "AliAODRecoDecay.h"
  44 #include "AliAODRecoDecayHF.h"
  45 #include "AliAODRecoDecayHF2Prong.h"
  46 #include "AliAODRecoDecayHF3Prong.h"
  47 #include "AliAODRecoDecayHF4Prong.h"
  48 #include "AliAODRecoCascadeHF.h"
  49 #include "AliAnalysisFilter.h"
  50 #include "AliAnalysisVertexingHF.h"
  51
  52 ClassImp(AliAnalysisVertexingHF)
  53
  54 //----------------------------------------------------------------------------
  55 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
  56 fInputAOD(kFALSE),
  57 fBzkG(0.),
  58 fSecVtxWithKF(kFALSE),
  59 fRecoPrimVtxSkippingTrks(kFALSE),
  60 fRmTrksFromPrimVtx(kFALSE),
  61 fV1(0x0),
  62 fD0toKpi(kTRUE),
  63 fJPSItoEle(kTRUE),
  64 f3Prong(kTRUE),
  65 f4Prong(kTRUE),
  66 fDstar(kTRUE),
  67 fLikeSign(kFALSE),
  68 fTrackFilter(0x0),
  69 fTrackFilterSoftPi(0x0),
  70 fFindVertexForDstar(kTRUE)
  71 {
  72   // Default constructor
  73
  74   SetD0toKpiCuts();
  75   SetBtoJPSICuts();
  76   SetDplusCuts();
  77   SetDsCuts();
  78   SetLcCuts();
  79   SetDstarCuts();
  80   SetD0to4ProngsCuts();
  81 }
  82 //--------------------------------------------------------------------------
  83 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
  84 TNamed(source),
  85 fInputAOD(source.fInputAOD),
  86 fBzkG(source.fBzkG),
  87 fSecVtxWithKF(source.fSecVtxWithKF),
  88 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
  89 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
  90 fV1(source.fV1),
  91 fD0toKpi(source.fD0toKpi),
  92 fJPSItoEle(source.fJPSItoEle),
  93 f3Prong(source.f3Prong),
  94 f4Prong(source.f4Prong),
  95 fDstar(source.fDstar),
  96 fLikeSign(source.fLikeSign),
  97 fTrackFilter(source.fTrackFilter),
  98 fTrackFilterSoftPi(source.fTrackFilterSoftPi),
  99 fFindVertexForDstar(source.fFindVertexForDstar)
 100 {
 101   //
 102   // Copy constructor
 103   //
 104   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
 105   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
 106   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
 107   for(Int_t i=0; i<13; i++) fDsCuts[i]=source.fDsCuts[i];
 108   for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
 109   for(Int_t i=0; i<5; i++)  fDstarCuts[i]=source.fDstarCuts[i];
 110   for(Int_t i=0; i<9; i++)  fD0to4ProngsCuts[i]=source.fD0to4ProngsCuts[i];
 111 }
 112 //--------------------------------------------------------------------------
 113 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
 114 {
 115   //
 116   // assignment operator
 117   //
 118   if(&source == this) return *this;
 119   fInputAOD = source.fInputAOD;
 120   fBzkG = source.fBzkG;
 121   fSecVtxWithKF = source.fSecVtxWithKF;
 122   fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
 123   fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
 124   fV1 = source.fV1;
 125   fD0toKpi = source.fD0toKpi;
 126   fJPSItoEle = source.fJPSItoEle;
 127   f3Prong = source.f3Prong;
 128   f4Prong = source.f4Prong;
 129   fDstar = source.fDstar;
 130   fLikeSign = source.fLikeSign;
 131   fTrackFilter = source.fTrackFilter;
 132   fTrackFilterSoftPi = source.fTrackFilterSoftPi;
 133   fFindVertexForDstar = source.fFindVertexForDstar;
 134
 135   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
 136   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
 137   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
 138   for(Int_t i=0; i<13; i++) fDsCuts[i]=source.fDsCuts[i];
 139   for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
 140   for(Int_t i=0; i<5; i++)  fDstarCuts[i]=source.fDstarCuts[i];
 141   for(Int_t i=0; i<9; i++)  fD0to4ProngsCuts[i]=source.fD0to4ProngsCuts[i];
 142
 143   return *this;
 144 }
 145 //----------------------------------------------------------------------------
 146 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
 147   // Destructor
 148   if(fV1) { delete fV1; fV1=0; }
 149   if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
 150   if(fTrackFilterSoftPi) { delete fTrackFilterSoftPi; fTrackFilterSoftPi=0; }
 151 }
 152 //----------------------------------------------------------------------------
 153 void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
 154                                             TClonesArray *aodVerticesHFTClArr,
 155                                             TClonesArray *aodD0toKpiTClArr,
 156                                             TClonesArray *aodJPSItoEleTClArr,
 157                                             TClonesArray *aodCharm3ProngTClArr,
 158                                             TClonesArray *aodCharm4ProngTClArr,
 159                                             TClonesArray *aodDstarTClArr,
 160                                             TClonesArray *aodLikeSign2ProngTClArr,
 161                                             TClonesArray *aodLikeSign3ProngTClArr)
 162 {
 163   // Find heavy-flavour vertex candidates
 164   // Input:  ESD or AOD
 165   // Output: AOD (additional branches added)
 166
 167
 168   TString evtype = event->IsA()->GetName();
 169   fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
 170
 171   if(fInputAOD) {
 172     AliDebug(2,"Creating HF candidates from AOD");
 173   } else {
 174     AliDebug(2,"Creating HF candidates from ESD");
 175   }
 176
 177   if(!aodVerticesHFTClArr) {
 178     printf("ERROR: no aodVerticesHFTClArr");
 179     return;
 180   }
 181   if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
 182     printf("ERROR: no aodD0toKpiTClArr");
 183     return;
 184   }
 185   if(fJPSItoEle && !aodJPSItoEleTClArr) {
 186     printf("ERROR: no aodJPSItoEleTClArr");
 187     return;
 188   }
 189   if(f3Prong && !aodCharm3ProngTClArr) {
 190     printf("ERROR: no aodCharm3ProngTClArr");
 191     return;
 192   }
 193   if(f4Prong && !aodCharm4ProngTClArr) {
 194     printf("ERROR: no aodCharm4ProngTClArr");
 195     return;
 196   }
 197   if(fDstar && !aodDstarTClArr) {
 198     printf("ERROR: no aodDstarTClArr");
 199     return;
 200   }
 201   if(fLikeSign && !aodLikeSign2ProngTClArr) {
 202     printf("ERROR: no aodLikeSign2ProngTClArr");
 203     return;
 204   }
 205   if(fLikeSign && f3Prong && !aodLikeSign3ProngTClArr) {
 206     printf("ERROR: no aodLikeSign2ProngTClArr");
 207     return;
 208   }
 209
 210   // delete candidates from previous event and create references
 211   Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0,iLikeSign2Prong=0,iLikeSign3Prong=0;
 212   aodVerticesHFTClArr->Delete();
 213   iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
 214   TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
 215   if(fD0toKpi || fDstar)   {
 216     aodD0toKpiTClArr->Delete();
 217     iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
 218   }
 219   if(fJPSItoEle) {
 220     aodJPSItoEleTClArr->Delete();
 221     iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
 222   }
 223   if(f3Prong) {
 224     aodCharm3ProngTClArr->Delete();
 225     i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
 226   }
 227   if(f4Prong) {
 228     aodCharm4ProngTClArr->Delete();
 229     i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
 230   }
 231   if(fDstar) {
 232     aodDstarTClArr->Delete();
 233     iDstar = aodDstarTClArr->GetEntriesFast();
 234   }
 235   if(fLikeSign) {
 236     aodLikeSign2ProngTClArr->Delete();
 237     iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast();
 238   }
 239   if(fLikeSign && f3Prong) {
 240     aodLikeSign3ProngTClArr->Delete();
 241     iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast();
 242   }
 243
 244   TClonesArray &aodD0toKpiRef        = *aodD0toKpiTClArr;
 245   TClonesArray &aodJPSItoEleRef      = *aodJPSItoEleTClArr;
 246   TClonesArray &aodCharm3ProngRef    = *aodCharm3ProngTClArr;
 247   TClonesArray &aodCharm4ProngRef    = *aodCharm4ProngTClArr;
 248   TClonesArray &aodDstarRef          = *aodDstarTClArr;
 249   TClonesArray &aodLikeSign2ProngRef = *aodLikeSign2ProngTClArr;
 250   TClonesArray &aodLikeSign3ProngRef = *aodLikeSign3ProngTClArr;
 251
 252
 253   AliAODRecoDecayHF2Prong *io2Prong  = 0;
 254   AliAODRecoDecayHF3Prong *io3Prong  = 0;
 255   AliAODRecoDecayHF4Prong *io4Prong  = 0;
 256   AliAODRecoCascadeHF     *ioCascade = 0;
 257
 258   Int_t    iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries;
 259   Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaCasc;
 260   Bool_t   okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
 261   Bool_t   okDstar=kFALSE,okD0fromDstar=kFALSE;
 262   AliESDtrack *postrack1 = 0;
 263   AliESDtrack *postrack2 = 0;
 264   AliESDtrack *negtrack1 = 0;
 265   AliESDtrack *negtrack2 = 0;
 266   AliESDtrack *trackPi   = 0;
 267   Double_t dcaMax = fD0toKpiCuts[1];
 268   if(dcaMax < fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1];
 269   if(dcaMax < fDplusCuts[11])  dcaMax=fDplusCuts[11];
 270   if(dcaMax < fD0to4ProngsCuts[1])  dcaMax=fD0to4ProngsCuts[1];
 271   AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
 272
 273
 274   // get Bz
 275   fBzkG = (Double_t)event->GetMagneticField();
 276
 277   trkEntries = (Int_t)event->GetNumberOfTracks();
 278   AliDebug(1,Form(" Number of tracks: %d",trkEntries));
 279
 280   if(trkEntries<2) {
 281     AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
 282     return;
 283   }
 284
 285   const AliVVertex *primary = event->GetPrimaryVertex();
 286   if(!primary) {
 287     AliDebug(1," No primary vertex from tracks");
 288     return;
 289   }
 290   TString primTitle = primary->GetTitle();
 291   if(!primTitle.Contains("VertexerTracks") ||
 292      primary->GetNContributors()<=0) {
 293     AliDebug(1," No primary vertex from tracks");
 294     return;
 295   }
 296
 297   // call function that applies sigle-track selection,
 298   // for displaced tracks and soft pions (both charges) for D*,
 299   // and retrieves primary vertex
 300   TObjArray seleTrksArray(trkEntries);
 301   UChar_t  *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi
 302   Int_t     nSeleTrks=0;
 303   SelectTracksAndCopyVertex(event,seleTrksArray,nSeleTrks,seleFlags);
 304
 305   AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
 306
 307   TObjArray *twoTrackArray1    = new TObjArray(2);
 308   TObjArray *twoTrackArray2    = new TObjArray(2);
 309   TObjArray *twoTrackArrayCasc = new TObjArray(2);
 310   TObjArray *threeTrackArray   = new TObjArray(3);
 311   TObjArray *fourTrackArray    = new TObjArray(4);
 312
 313   Double_t dispersion;
 314   Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
 315
 316   AliAODRecoDecayHF   *rd = 0;
 317   AliAODRecoCascadeHF *rc = 0;
 318
 319   // LOOP ON  POSITIVE  TRACKS
 320   for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
 321
 322     if(iTrkP1%1==0) AliDebug(1,Form("  1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));
 323
 324     // get track from tracks array
 325     postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
 326
 327     if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
 328
 329     if(postrack1->Charge()<0 && !fLikeSign) continue;
 330
 331     // LOOP ON  NEGATIVE  TRACKS
 332     for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
 333
 334       if(iTrkN1%1==0) AliDebug(1,Form("    1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));
 335
 336       if(iTrkN1==iTrkP1) continue;
 337
 338       // get track from tracks array
 339       negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
 340
 341       if(negtrack1->Charge()>0 && !fLikeSign) continue;
 342
 343       if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
 344
 345       if(postrack1->Charge()==negtrack1->Charge()) { // like-sign
 346         isLikeSign2Prong=kTRUE;
 347         if(!fLikeSign)    continue;
 348         if(iTrkN1<iTrkP1) continue; // this is needed to avoid double-counting of like-sign
 349       } else { // unlike-sign
 350         isLikeSign2Prong=kFALSE;
 351         if(postrack1->Charge()<0 || negtrack1->Charge()>0) continue;  // this is needed to avoid double-counting of unlike-sign
 352       }
 353
 354       // back to primary vertex
 355       postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 356       negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 357       //printf("********** %d %d\n",postrack1->GetID(),negtrack1->GetID());
 358
 359       // DCA between the two tracks
 360       dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 361       if(dcap1n1>dcaMax) { negtrack1=0; continue; }
 362
 363       // Vertexing
 364       twoTrackArray1->AddAt(postrack1,0);
 365       twoTrackArray1->AddAt(negtrack1,1);
 366       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
 367       if(!vertexp1n1) {
 368         twoTrackArray1->Clear();
 369         negtrack1=0;
 370         continue;
 371       }
 372
 373       // 2 prong candidate
 374       if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) {
 375
 376         io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
 377
 378         if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI) || (isLikeSign2Prong && (okD0 || okJPSI))) {
 379           // add the vertex and the decay to the AOD
 380           AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
 381           if(fInputAOD) AddDaughterRefs(v2Prong,event,twoTrackArray1);
 382           if(!isLikeSign2Prong) {
 383             if(okD0) {
 384               rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
 385               rd->SetSecondaryVtx(v2Prong);
 386               v2Prong->SetParent(rd);
 387             }
 388             if(okJPSI) {
 389               rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
 390               rd->SetSecondaryVtx(v2Prong);
 391               if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
 392             }
 393           } else { // isLikeSign2Prong
 394             rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
 395             rd->SetSecondaryVtx(v2Prong);
 396             v2Prong->SetParent(rd);
 397           }
 398         }
 399         // D* candidates
 400         if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
 401           // write references in io2Prong
 402           if(fInputAOD) {
 403             AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
 404           } else {
 405             vertexp1n1->AddDaughter(postrack1);
 406             vertexp1n1->AddDaughter(negtrack1);
 407           }
 408           io2Prong->SetSecondaryVtx(vertexp1n1);
 409           //printf("--->  %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
 410           // create a track from the D0
 411           AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
 412
 413           // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
 414           for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
 415
 416             if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
 417
 418             if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
 419
 420             if(iTrkSoftPi%1==0) AliDebug(1,Form("    1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));
 421
 422             trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
 423             if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
 424
 425             // get track from tracks array
 426             trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
 427             trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
 428
 429             twoTrackArrayCasc->AddAt(trackPi,0);
 430             twoTrackArrayCasc->AddAt(trackD0,1);
 431
 432             AliAODVertex *vertexCasc = 0;
 433
 434             if(fFindVertexForDstar) {
 435               // DCA between the two tracks
 436               dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
 437               // Vertexing
 438               vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
 439             } else {
 440               // assume Dstar decays at the primary vertex
 441               Double_t pos[3],cov[6],chi2perNDF;
 442               fV1->GetXYZ(pos);
 443               fV1->GetCovMatrix(cov);
 444               chi2perNDF = fV1->GetChi2toNDF();
 445               vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
 446               dcaCasc = 0.;
 447             }
 448             if(!vertexCasc) {
 449               twoTrackArrayCasc->Clear();
 450               trackPi=0;
 451               continue;
 452             }
 453
 454             ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
 455             if(okDstar) {
 456               // add the D0 to the AOD (if not already done)
 457               if(!okD0) {
 458                 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
 459                 if(fInputAOD) AddDaughterRefs(v2Prong,event,twoTrackArray1);
 460                 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
 461                 rd->SetSecondaryVtx(v2Prong);
 462                 v2Prong->SetParent(rd);
 463                 okD0=kTRUE; // this is done to add it only once
 464               }
 465               // add the vertex and the cascade to the AOD
 466               AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
 467               if(fInputAOD) {
 468                 AddDaughterRefs(vCasc,event,twoTrackArrayCasc); // add the pion
 469               } else {
 470                 vCasc->AddDaughter(rd); // just to fill ref #0
 471               }
 472               vCasc->AddDaughter(rd); // add the D0 (in ref #1)
 473               rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
 474               rc->SetSecondaryVtx(vCasc);
 475               vCasc->SetParent(rc);
 476             }
 477             twoTrackArrayCasc->Clear();
 478             trackPi=0;
 479             if(ioCascade) {delete ioCascade; ioCascade=NULL;}
 480             delete vertexCasc; vertexCasc=NULL;
 481           } // end loop on soft pi tracks
 482
 483           if(trackD0) {delete trackD0; trackD0=NULL;}
 484
 485         }
 486         if(io2Prong) {delete io2Prong; io2Prong=NULL;}
 487       }
 488
 489       twoTrackArray1->Clear();
 490       if( (!f3Prong && !f4Prong) ||
 491           (isLikeSign2Prong && !f3Prong) ) {
 492         negtrack1=0;
 493         delete vertexp1n1;
 494         continue;
 495       }
 496
 497
 498       // 2nd LOOP  ON  POSITIVE  TRACKS
 499       for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
 500
 501         if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
 502
 503         if(iTrkP2%1==0) AliDebug(1,Form("    2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));
 504
 505         // get track from tracks array
 506         postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
 507
 508         if(postrack2->Charge()<0) continue;
 509
 510         if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
 511
 512         if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
 513           if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
 514             isLikeSign3Prong=kTRUE;
 515           } else { // not ok
 516             continue;
 517           }
 518         } else { // normal triplet (+-+)
 519           isLikeSign3Prong=kFALSE;
 520         }
 521
 522         // back to primary vertex
 523         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 524         postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 525         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 526         //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
 527
 528         dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 529         if(dcap2n1>dcaMax) { postrack2=0; continue; }
 530         dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
 531         if(dcap1p2>dcaMax) { postrack2=0; continue; }
 532
 533         // Vertexing
 534         twoTrackArray2->AddAt(postrack2,0);
 535         twoTrackArray2->AddAt(negtrack1,1);
 536         AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
 537         if(!vertexp2n1) {
 538           twoTrackArray2->Clear();
 539           postrack2=0;
 540           continue;
 541         }
 542
 543         // 3 prong candidates
 544         if(f3Prong) {
 545           if(postrack2->Charge()>0) {
 546             threeTrackArray->AddAt(postrack1,0);
 547             threeTrackArray->AddAt(negtrack1,1);
 548             threeTrackArray->AddAt(postrack2,2);
 549           } else {
 550             threeTrackArray->AddAt(negtrack1,0);
 551             threeTrackArray->AddAt(postrack1,1);
 552             threeTrackArray->AddAt(postrack2,2);
 553           }
 554
 555           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
 556           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
 557           if(ok3Prong) {
 558             AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
 559             if(fInputAOD) AddDaughterRefs(v3Prong,event,threeTrackArray);
 560             if(!isLikeSign3Prong) {
 561               rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
 562               rd->SetSecondaryVtx(v3Prong);
 563               v3Prong->SetParent(rd);
 564             } else { // isLikeSign3Prong
 565               rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
 566               rd->SetSecondaryVtx(v3Prong);
 567               v3Prong->SetParent(rd);
 568             }
 569           }
 570           if(io3Prong) {delete io3Prong; io3Prong=NULL;}
 571           if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
 572         }
 573
 574         // 4 prong candidates
 575         if(f4Prong
 576            // don't make 4 prong with like-sign pairs and triplets
 577            && !isLikeSign2Prong && !isLikeSign3Prong
 578            // track-to-track dca cuts already now
 579            && dcap1n1 < fD0to4ProngsCuts[1]
 580            && dcap2n1 < fD0to4ProngsCuts[1]) {
 581
 582           // back to primary vertex
 583           postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 584           postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 585           negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 586           // Vertexing for these 3 (can be taken from above?)
 587           threeTrackArray->AddAt(postrack1,0);
 588           threeTrackArray->AddAt(negtrack1,1);
 589           threeTrackArray->AddAt(postrack2,2);
 590           AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
 591
 592           // 3rd LOOP  ON  NEGATIVE  TRACKS (for 4 prong)
 593           for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
 594
 595             if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
 596
 597             if(iTrkN2%1==0) AliDebug(1,Form("    3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
 598
 599             // get track from tracks array
 600             negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
 601
 602             if(negtrack2->Charge()>0) continue;
 603
 604             if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
 605
 606             // back to primary vertex
 607             postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 608             postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 609             negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 610             negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 611             dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 612             if(dcap1n2 > fD0to4ProngsCuts[1]) { negtrack2=0; continue; }
 613             dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 614             if(dcap2n2 > fD0to4ProngsCuts[1]) { negtrack2=0; continue; }
 615
 616             // Vertexing
 617             fourTrackArray->AddAt(postrack1,0);
 618             fourTrackArray->AddAt(negtrack1,1);
 619             fourTrackArray->AddAt(postrack2,2);
 620             fourTrackArray->AddAt(negtrack2,3);
 621
 622             AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
 623             io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
 624             if(ok4Prong) {
 625               AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
 626               if(fInputAOD) AddDaughterRefs(v4Prong,event,fourTrackArray);
 627               rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
 628               rd->SetSecondaryVtx(v4Prong);
 629               v4Prong->SetParent(rd);
 630             }
 631
 632             if(io4Prong) {delete io4Prong; io4Prong=NULL;}
 633             if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;}
 634             fourTrackArray->Clear();
 635             negtrack2 = 0;
 636
 637           } // end loop on negative tracks
 638
 639           threeTrackArray->Clear();
 640           delete vertexp1n1p2;
 641
 642         }
 643
 644         postrack2 = 0;
 645         delete vertexp2n1;
 646
 647       } // end 2nd loop on positive tracks
 648
 649       twoTrackArray2->Clear();
 650
 651       // 2nd LOOP  ON  NEGATIVE  TRACKS (for 3 prong -+-)
 652       for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
 653
 654         if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
 655
 656         if(iTrkN2%1==0) AliDebug(1,Form("    2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
 657
 658         // get track from tracks array
 659         negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
 660
 661         if(negtrack2->Charge()>0) continue;
 662
 663         if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
 664
 665         if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
 666           if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
 667             isLikeSign3Prong=kTRUE;
 668           } else { // not ok
 669             continue;
 670           }
 671         } else { // normal triplet (-+-)
 672           isLikeSign3Prong=kFALSE;
 673         }
 674
 675         // back to primary vertex
 676         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 677         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 678         negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 679         //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
 680
 681         dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 682         if(dcap1n2>dcaMax) { negtrack2=0; continue; }
 683         dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 684         if(dcan1n2>dcaMax) { negtrack2=0; continue; }
 685
 686         // Vertexing
 687         twoTrackArray2->AddAt(postrack1,0);
 688         twoTrackArray2->AddAt(negtrack2,1);
 689
 690         AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
 691         if(!vertexp1n2) {
 692           twoTrackArray2->Clear();
 693           negtrack2=0;
 694           continue;
 695         }
 696
 697         if(f3Prong) {
 698           threeTrackArray->AddAt(negtrack1,0);
 699           threeTrackArray->AddAt(postrack1,1);
 700           threeTrackArray->AddAt(negtrack2,2);
 701           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
 702           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
 703           if(ok3Prong) {
 704             AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
 705             if(fInputAOD) AddDaughterRefs(v3Prong,event,threeTrackArray);
 706             if(!isLikeSign3Prong) {
 707               rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
 708               rd->SetSecondaryVtx(v3Prong);
 709               v3Prong->SetParent(rd);
 710             } else { // isLikeSign3Prong
 711               rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
 712               rd->SetSecondaryVtx(v3Prong);
 713               v3Prong->SetParent(rd);
 714             }
 715           }
 716           if(io3Prong) {delete io3Prong; io3Prong=NULL;}
 717           if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
 718         }
 719
 720         negtrack2 = 0;
 721         delete vertexp1n2;
 722
 723       } // end 2nd loop on negative tracks
 724
 725       twoTrackArray2->Clear();
 726
 727       negtrack1 = 0;
 728       delete vertexp1n1;
 729     } // end 1st loop on negative tracks
 730
 731     postrack1 = 0;
 732   }  // end 1st loop on positive tracks
 733
 734
 735   AliDebug(1,Form(" Total HF vertices in event = %d;",
 736                   (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
 737   if(fD0toKpi) {
 738     AliDebug(1,Form(" D0->Kpi in event = %d;",
 739                     (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
 740   }
 741   if(fJPSItoEle) {
 742     AliDebug(1,Form(" JPSI->ee in event = %d;",
 743                     (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
 744   }
 745   if(f3Prong) {
 746     AliDebug(1,Form(" Charm->3Prong in event = %d;",
 747                     (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
 748   }
 749   if(f4Prong) {
 750     AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
 751                     (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
 752   }
 753   if(fDstar) {
 754     AliDebug(1,Form(" D*->D0pi in event = %d;\n",
 755                     (Int_t)aodDstarTClArr->GetEntriesFast()));
 756   }
 757   if(fLikeSign) {
 758     AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
 759                     (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
 760   }
 761   if(fLikeSign && f3Prong) {
 762     AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
 763                     (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
 764   }
 765
 766
 767   twoTrackArray1->Delete();  delete twoTrackArray1;
 768   twoTrackArray2->Delete();  delete twoTrackArray2;
 769   twoTrackArrayCasc->Delete();  delete twoTrackArrayCasc;
 770   threeTrackArray->Clear();
 771   threeTrackArray->Delete(); delete threeTrackArray;
 772   fourTrackArray->Delete();  delete fourTrackArray;
 773   delete [] seleFlags; seleFlags=NULL;
 774
 775   if(fInputAOD) {
 776     seleTrksArray.Delete();
 777   }
 778
 779   return;
 780 }
 781 //----------------------------------------------------------------------------
 782 void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,AliVEvent *event,
 783                                              TObjArray *trkArray) const
 784 {
 785   // Add the AOD tracks as daughters of the vertex (TRef)
 786
 787   Int_t nTrks = trkArray->GetEntriesFast();
 788
 789   AliExternalTrackParam *track = 0;
 790   AliAODTrack *aodTrack = 0;
 791   Int_t id;
 792
 793   for(Int_t i=0; i<nTrks; i++) {
 794     track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
 795     id = (Int_t)track->GetID();
 796     //printf("---> %d\n",id);
 797     if(id<0) continue; // this track is a AliAODRecoDecay
 798     aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
 799     v->AddDaughter(aodTrack);
 800   }
 801
 802   return;
 803 }
 804 //----------------------------------------------------------------------------
 805 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
 806                                    TObjArray *twoTrackArray,AliVEvent *event,
 807                                    AliAODVertex *secVert,
 808                                    AliAODRecoDecayHF2Prong *rd2Prong,
 809                                    Double_t dca,
 810                                    Bool_t &okDstar) const
 811 {
 812   // Make the cascade as a 2Prong decay and check if it passes Dstar
 813   // reconstruction cuts
 814
 815   okDstar = kFALSE;
 816
 817   Bool_t dummy1,dummy2,dummy3;
 818
 819   // We use Make2Prong to construct the AliAODRecoCascadeHF
 820   // (which inherits from AliAODRecoDecayHF2Prong)
 821   AliAODRecoCascadeHF *theCascade =
 822     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
 823                                      dummy1,dummy2,dummy3);
 824   if(!theCascade) return 0x0;
 825
 826   // charge
 827   AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
 828   theCascade->SetCharge(trackPi->Charge());
 829
 830   //--- selection cuts
 831   //
 832   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
 833   tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
 834   tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
 835   AliAODVertex *primVertexAOD=0;
 836   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
 837     // take event primary vertex
 838     primVertexAOD = PrimaryVertex();
 839     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
 840     rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
 841   }
 842   // select D*->D0pi
 843   if(fDstar) {
 844     Bool_t testD0=kTRUE;
 845     okDstar = tmpCascade->SelectDstar(fDstarCuts,fD0fromDstarCuts,testD0);
 846   }
 847   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
 848   tmpCascade->UnsetOwnPrimaryVtx();
 849   delete tmpCascade; tmpCascade=NULL;
 850   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
 851     rd2Prong->UnsetOwnPrimaryVtx();
 852   }
 853   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
 854   //---
 855
 856   return theCascade;
 857 }
 858 //-----------------------------------------------------------------------------
 859 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
 860                                    TObjArray *twoTrackArray,AliVEvent *event,
 861                                    AliAODVertex *secVert,Double_t dca,
 862                                    Bool_t &okD0,Bool_t &okJPSI,
 863                                    Bool_t &okD0fromDstar) const
 864 {
 865   // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
 866   // reconstruction cuts
 867   // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
 868
 869   okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
 870
 871   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
 872
 873   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
 874   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
 875
 876   // propagate tracks to secondary vertex, to compute inv. mass
 877   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
 878   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
 879
 880   Double_t momentum[3];
 881   postrack->GetPxPyPz(momentum);
 882   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 883   negtrack->GetPxPyPz(momentum);
 884   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 885
 886
 887   // invariant mass cut (try to improve coding here..)
 888   Bool_t okMassCut=kFALSE;
 889   if(!okMassCut && fD0toKpi)   if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
 890   if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
 891   if(!okMassCut && fDstar)     if(SelectInvMass(3,2,px,py,pz)) okMassCut=kTRUE;
 892   if(!okMassCut) {
 893     AliDebug(2," candidate didn't pass mass cut");
 894     return 0x0;
 895   }
 896
 897   // primary vertex to be used by this candidate
 898   AliAODVertex *primVertexAOD  = PrimaryVertex(twoTrackArray,event);
 899   if(!primVertexAOD) return 0x0;
 900
 901   Double_t d0z0[2],covd0z0[3];
 902   postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 903   d0[0] = d0z0[0];
 904   d0err[0] = TMath::Sqrt(covd0z0[0]);
 905   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 906   d0[1] = d0z0[0];
 907   d0err[1] = TMath::Sqrt(covd0z0[0]);
 908
 909   // create the object AliAODRecoDecayHF2Prong
 910   AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
 911   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
 912   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
 913   the2Prong->SetProngIDs(2,id);
 914   delete primVertexAOD; primVertexAOD=NULL;
 915
 916
 917   // select D0->Kpi
 918   Int_t checkD0,checkD0bar;
 919   if(fD0toKpi)   okD0          = the2Prong->SelectD0(fD0toKpiCuts,checkD0,checkD0bar);
 920   //if(fDebug && fD0toKpi) printf("   %d\n",(Int_t)okD0);
 921   // select J/psi from B
 922   Int_t checkJPSI;
 923   if(fJPSItoEle) okJPSI        = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI);
 924   //if(fDebug && fJPSItoEle) printf("   %d\n",(Int_t)okJPSI);
 925   // select D0->Kpi from Dstar
 926   if(fDstar)     okD0fromDstar = the2Prong->SelectD0(fD0fromDstarCuts,checkD0,checkD0bar);
 927   //if(fDebug && fDstar) printf("   %d\n",(Int_t)okD0fromDstar);
 928
 929   // remove the primary vertex (was used only for selection)
 930   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
 931     the2Prong->UnsetOwnPrimaryVtx();
 932   }
 933
 934   // get PID info from ESD
 935   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
 936   if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
 937   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
 938   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
 939   Double_t esdpid[10];
 940   for(Int_t i=0;i<5;i++) {
 941     esdpid[i]   = esdpid0[i];
 942     esdpid[5+i] = esdpid1[i];
 943   }
 944   the2Prong->SetPID(2,esdpid);
 945
 946   return the2Prong;
 947 }
 948 //----------------------------------------------------------------------------
 949 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
 950                              TObjArray *threeTrackArray,AliVEvent *event,
 951                              AliAODVertex *secVert,Double_t dispersion,
 952                              AliAODVertex *vertexp1n1,AliAODVertex *vertexp2n1,
 953                              Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
 954                              Bool_t &ok3Prong) const
 955 {
 956   // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
 957   // reconstruction cuts
 958   // E.Bruna, F.Prino
 959
 960
 961   ok3Prong=kFALSE;
 962   if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
 963
 964   Double_t px[3],py[3],pz[3],d0[3],d0err[3];
 965   Double_t momentum[3];
 966
 967
 968   AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
 969   AliESDtrack *negtrack  = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
 970   AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
 971
 972   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
 973   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
 974   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
 975   postrack1->GetPxPyPz(momentum);
 976   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 977   negtrack->GetPxPyPz(momentum);
 978   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 979   postrack2->GetPxPyPz(momentum);
 980   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
 981
 982   // invariant mass cut for D+, Ds, Lc
 983   Bool_t okMassCut=kFALSE;
 984   if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
 985   if(!okMassCut) {
 986     AliDebug(2," candidate didn't pass mass cut");
 987     return 0x0;
 988   }
 989
 990   // primary vertex to be used by this candidate
 991   AliAODVertex *primVertexAOD  = PrimaryVertex(threeTrackArray,event);
 992   if(!primVertexAOD) return 0x0;
 993
 994   Double_t d0z0[2],covd0z0[3];
 995   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 996   d0[0]=d0z0[0];
 997   d0err[0] = TMath::Sqrt(covd0z0[0]);
 998   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 999   d0[1]=d0z0[0];
1000   d0err[1] = TMath::Sqrt(covd0z0[0]);
1001   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1002   d0[2]=d0z0[0];
1003   d0err[2] = TMath::Sqrt(covd0z0[0]);
1004
1005
1006   // create the object AliAODRecoDecayHF3Prong
1007   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1008   Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1009   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]));
1010   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]));
1011   Short_t charge=(Short_t)(postrack1->Charge()*postrack2->Charge()*negtrack->Charge());
1012   AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1013   the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1014   UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1015   the3Prong->SetProngIDs(3,id);
1016
1017   delete primVertexAOD; primVertexAOD=NULL;
1018
1019   // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1020   if(f3Prong) {
1021     ok3Prong = kFALSE;
1022     Int_t ok1,ok2;
1023     if(the3Prong->SelectDplus(fDplusCuts))   ok3Prong = kTRUE;
1024     if(the3Prong->SelectDs(fDsCuts,ok1,ok2)) ok3Prong = kTRUE;
1025     if(the3Prong->SelectLc(fLcCuts,ok1,ok2)) ok3Prong = kTRUE;
1026   }
1027   //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1028
1029   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1030     the3Prong->UnsetOwnPrimaryVtx();
1031   }
1032
1033   // get PID info from ESD
1034   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1035   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1036   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1037   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1038   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1039   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1040
1041   Double_t esdpid[15];
1042   for(Int_t i=0;i<5;i++) {
1043     esdpid[i]    = esdpid0[i];
1044     esdpid[5+i]  = esdpid1[i];
1045     esdpid[10+i] = esdpid2[i];
1046   }
1047   the3Prong->SetPID(3,esdpid);
1048
1049   return the3Prong;
1050 }
1051 //----------------------------------------------------------------------------
1052 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1053                              TObjArray *fourTrackArray,AliVEvent *event,
1054                              AliAODVertex *secVert,
1055                              AliAODVertex *vertexp1n1,
1056                              AliAODVertex *vertexp1n1p2,
1057                              Double_t dcap1n1,Double_t dcap1n2,
1058                              Double_t dcap2n1,Double_t dcap2n2,
1059                              Bool_t &ok4Prong) const
1060 {
1061   // Make 4Prong candidates and check if they pass D0toKpipipi
1062   // reconstruction cuts
1063   // G.E.Bruno, R.Romita
1064
1065   ok4Prong=kFALSE;
1066   if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
1067
1068   Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1069
1070   AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1071   AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1072   AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1073   AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1074
1075   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1076   negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1077   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1078   negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1079
1080   Double_t momentum[3];
1081   postrack1->GetPxPyPz(momentum);
1082   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1083   negtrack1->GetPxPyPz(momentum);
1084   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1085   postrack2->GetPxPyPz(momentum);
1086   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1087   negtrack2->GetPxPyPz(momentum);
1088   px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1089
1090   // invariant mass cut for rho or D0 (try to improve coding here..)
1091   Bool_t okMassCut=kFALSE;
1092   if(!okMassCut && fD0to4ProngsCuts[8]==0.){      //no PID, to be implemented with PID
1093     if(SelectInvMass(4,4,px,py,pz)) okMassCut=kTRUE;
1094   }
1095   if(!okMassCut) {
1096     //if(fDebug) printf(" candidate didn't pass mass cut\n");
1097     //printf(" candidate didn't pass mass cut\n");
1098     return 0x0;
1099   }
1100
1101   // primary vertex to be used by this candidate
1102   AliAODVertex *primVertexAOD  = PrimaryVertex(fourTrackArray,event);
1103   if(!primVertexAOD) return 0x0;
1104
1105   Double_t d0z0[2],covd0z0[3];
1106   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1107   d0[0]=d0z0[0];
1108   d0err[0] = TMath::Sqrt(covd0z0[0]);
1109   negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1110   d0[1]=d0z0[0];
1111   d0err[1] = TMath::Sqrt(covd0z0[0]);
1112   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1113   d0[2]=d0z0[0];
1114   d0err[2] = TMath::Sqrt(covd0z0[0]);
1115   negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1116   d0[3]=d0z0[0];
1117   d0err[3] = TMath::Sqrt(covd0z0[0]);
1118
1119
1120   // create the object AliAODRecoDecayHF4Prong
1121   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1122   Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
1123   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]));
1124   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]));
1125   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]));
1126   Short_t charge=0;
1127   AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
1128   the4Prong->SetOwnPrimaryVtx(primVertexAOD);
1129   UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
1130   the4Prong->SetProngIDs(4,id);
1131
1132   delete primVertexAOD; primVertexAOD=NULL;
1133
1134   Int_t checkD0,checkD0bar;
1135   ok4Prong=the4Prong->SelectD0(fD0to4ProngsCuts,checkD0,checkD0bar);
1136
1137   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1138     the4Prong->UnsetOwnPrimaryVtx();
1139   }
1140
1141
1142   // get PID info from ESD
1143   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1144   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1145   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1146   if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
1147   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1148   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1149   Double_t esdpid3[5]={0.,0.,0.,0.,0.};
1150   if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
1151
1152   Double_t esdpid[20];
1153   for(Int_t i=0;i<5;i++) {
1154     esdpid[i]    = esdpid0[i];
1155     esdpid[5+i]  = esdpid1[i];
1156     esdpid[10+i] = esdpid2[i];
1157     esdpid[15+i] = esdpid3[i];
1158   }
1159   the4Prong->SetPID(4,esdpid);
1160
1161   return the4Prong;
1162 }
1163 //-----------------------------------------------------------------------------
1164 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(TObjArray *trkArray,
1165                                                     AliVEvent *event) const
1166 {
1167   // Returns primary vertex to be used for this candidate
1168
1169   AliESDVertex *vertexESD = 0;
1170   AliAODVertex *vertexAOD = 0;
1171
1172
1173   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1174     // primary vertex from the input event
1175
1176     vertexESD = new AliESDVertex(*fV1);
1177
1178   } else {
1179     // primary vertex specific to this candidate
1180
1181     Int_t nTrks = trkArray->GetEntriesFast();
1182     AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1183
1184     if(fRecoPrimVtxSkippingTrks) {
1185       // recalculating the vertex
1186
1187       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1188         Float_t diamondcovxy[3];
1189         event->GetDiamondCovXY(diamondcovxy);
1190         Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1191         Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.};
1192         AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1193         vertexer->SetVtxStart(diamond);
1194         delete diamond; diamond=NULL;
1195         if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1196           vertexer->SetOnlyFitter();
1197       }
1198       Int_t skipped[1000];
1199       Int_t nTrksToSkip=0,id;
1200       AliExternalTrackParam *t = 0;
1201       for(Int_t i=0; i<nTrks; i++) {
1202         t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1203         id = (Int_t)t->GetID();
1204         if(id<0) continue;
1205         skipped[nTrksToSkip++] = id;
1206       }
1207       // TEMPORARY FIX
1208       // For AOD, skip also tracks without covariance matrix
1209       if(fInputAOD) {
1210         Double_t covtest[21];
1211         for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1212           AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1213           if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1214             id = (Int_t)t->GetID();
1215             if(id<0) continue;
1216             skipped[nTrksToSkip++] = id;
1217           }
1218         }
1219       }
1220       //
1221       vertexer->SetSkipTracks(nTrksToSkip,skipped);
1222       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1223
1224     } else if(fRmTrksFromPrimVtx) {
1225       // removing the prongs tracks
1226
1227       TObjArray rmArray(nTrks);
1228       UShort_t *rmId = new UShort_t[nTrks];
1229       AliESDtrack *esdTrack = 0;
1230       AliESDtrack *t = 0;
1231       for(Int_t i=0; i<nTrks; i++) {
1232         t = (AliESDtrack*)trkArray->UncheckedAt(i);
1233         esdTrack = new AliESDtrack(*t);
1234         rmArray.AddLast(esdTrack);
1235         if(esdTrack->GetID()>=0) {
1236           rmId[i]=(UShort_t)esdTrack->GetID();
1237         } else {
1238           rmId[i]=9999;
1239         }
1240       }
1241       Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
1242       vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1243       delete [] rmId; rmId=NULL;
1244       rmArray.Delete();
1245
1246     }
1247
1248     if(!vertexESD) return vertexAOD;
1249     if(vertexESD->GetNContributors()<=0) {
1250       AliDebug(2,"vertexing failed");
1251       delete vertexESD; vertexESD=NULL;
1252       return vertexAOD;
1253     }
1254
1255     delete vertexer; vertexer=NULL;
1256
1257   }
1258
1259   // convert to AliAODVertex
1260   Double_t pos[3],cov[6],chi2perNDF;
1261   vertexESD->GetXYZ(pos); // position
1262   vertexESD->GetCovMatrix(cov); //covariance matrix
1263   chi2perNDF = vertexESD->GetChi2toNDF();
1264   delete vertexESD; vertexESD=NULL;
1265
1266   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1267
1268   return vertexAOD;
1269 }
1270 //-----------------------------------------------------------------------------
1271 void AliAnalysisVertexingHF::PrintStatus() const {
1272   // Print parameters being used
1273
1274   //printf("Preselections:\n");
1275   //   fTrackFilter->Dump();
1276   if(fSecVtxWithKF) {
1277     printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1278   } else {
1279     printf("Secondary vertex with AliVertexerTracks\n");
1280   }
1281   if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1282   if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1283   if(fD0toKpi) {
1284     printf("Reconstruct D0->Kpi candidates with cuts:\n");
1285     printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
1286     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
1287     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
1288     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
1289     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
1290     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
1291     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
1292     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
1293     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
1294   }
1295   if(fDstar) {
1296     if(fFindVertexForDstar) {
1297       printf("    Reconstruct a secondary vertex for the D*\n");
1298     } else {
1299       printf("    Assume the D* comes from the primary vertex\n");
1300     }
1301     printf("    |M-MD*| [GeV]    < %f\n",fDstarCuts[0]);
1302     printf("    |M_Kpipi-M_Kpi-(MD*-MD0)| [GeV]  < %f\n",fDstarCuts[1]);
1303     printf("    pTpisoft [GeV/c]    > %f\n",fDstarCuts[2]);
1304     printf("    pTpisoft [GeV/c]    < %f\n",fDstarCuts[3]);
1305     printf("    Theta(pisoft,D0plane) < %f\n",fDstarCuts[4]);
1306     printf("Reconstruct D*->D0pi candidates with cuts:\n");
1307     printf("   D0 from D* cuts:\n");
1308     printf("    |M-MD0| [GeV]    < %f\n",fD0fromDstarCuts[0]);
1309     printf("    dca    [cm]  < %f\n",fD0fromDstarCuts[1]);
1310     printf("    cosThetaStar     < %f\n",fD0fromDstarCuts[2]);
1311     printf("    pTK     [GeV/c]    > %f\n",fD0fromDstarCuts[3]);
1312     printf("    pTpi    [GeV/c]    > %f\n",fD0fromDstarCuts[4]);
1313     printf("    |d0K|  [cm]  < %f\n",fD0fromDstarCuts[5]);
1314     printf("    |d0pi| [cm]  < %f\n",fD0fromDstarCuts[6]);
1315     printf("    d0d0  [cm^2] < %f\n",fD0fromDstarCuts[7]);
1316     printf("    cosThetaPoint    > %f\n",fD0fromDstarCuts[8]);
1317   }
1318   if(fJPSItoEle) {
1319     printf("Reconstruct J/psi from B candidates with cuts:\n");
1320     printf("    |M-MJPSI| [GeV]    < %f\n",fBtoJPSICuts[0]);
1321     printf("    dca    [cm]  < %f\n",fBtoJPSICuts[1]);
1322     printf("    cosThetaStar     < %f\n",fBtoJPSICuts[2]);
1323     printf("    pTP     [GeV/c]    > %f\n",fBtoJPSICuts[3]);
1324     printf("    pTN    [GeV/c]    > %f\n",fBtoJPSICuts[4]);
1325     printf("    |d0P|  [cm]  < %f\n",fBtoJPSICuts[5]);
1326     printf("    |d0N| [cm]  < %f\n",fBtoJPSICuts[6]);
1327     printf("    d0d0  [cm^2] < %f\n",fBtoJPSICuts[7]);
1328     printf("    cosThetaPoint    > %f\n",fBtoJPSICuts[8]);
1329   }
1330   if(f3Prong) {
1331     printf("Reconstruct 3 prong candidates.\n");
1332     printf("  D+->Kpipi cuts:\n");
1333     printf("    |M-MD+| [GeV]    < %f\n",fDplusCuts[0]);
1334     printf("    pTK     [GeV/c]    > %f\n",fDplusCuts[1]);
1335     printf("    pTPi    [GeV/c]    > %f\n",fDplusCuts[2]);
1336     printf("    |d0K|  [cm]  > %f\n",fDplusCuts[3]);
1337     printf("    |d0Pi| [cm]  > %f\n",fDplusCuts[4]);
1338     printf("    dist12    [cm]  < %f\n",fDplusCuts[5]);
1339     printf("    sigmavert [cm]   < %f\n",fDplusCuts[6]);
1340     printf("    dist prim-sec [cm] > %f\n",fDplusCuts[7]);
1341     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDplusCuts[8]);
1342     printf("    cosThetaPoint    > %f\n",fDplusCuts[9]);
1343     printf("    Sum d0^2 [cm^2]  > %f\n",fDplusCuts[10]);
1344     printf("    dca cut [cm]  < %f\n",fDplusCuts[11]);
1345     printf("  Ds->KKpi cuts:\n");
1346     printf("    |M-MDs| [GeV]    < %f\n",fDsCuts[0]);
1347     printf("    pTK     [GeV/c]    > %f\n",fDsCuts[1]);
1348     printf("    pTPi    [GeV/c]    > %f\n",fDsCuts[2]);
1349     printf("    |d0K|  [cm]  > %f\n",fDsCuts[3]);
1350     printf("    |d0Pi| [cm]  > %f\n",fDsCuts[4]);
1351     printf("    dist12    [cm]  < %f\n",fDsCuts[5]);
1352     printf("    sigmavert [cm]   < %f\n",fDsCuts[6]);
1353     printf("    dist prim-sec [cm] > %f\n",fDsCuts[7]);
1354     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDsCuts[8]);
1355     printf("    cosThetaPoint    > %f\n",fDsCuts[9]);
1356     printf("    Sum d0^2 [cm^2]  > %f\n",fDsCuts[10]);
1357     printf("    dca cut [cm]  < %f\n",fDsCuts[11]);
1358     printf("    Inv. Mass  phi/K0* [GeV]  < %f\n",fDsCuts[12]);
1359     printf("  Lc->pKpi cuts:\n");
1360     printf("    |M-MLc| [GeV]    < %f\n",fLcCuts[0]);
1361     printf("    pTP     [GeV/c]    > %f\n",fLcCuts[1]);
1362     printf("    pTPi and pTK [GeV/c]    > %f\n",fLcCuts[2]);
1363     printf("    |d0P|  [cm]  > %f\n",fLcCuts[3]);
1364     printf("    |d0Pi| and |d0K| [cm]  > %f\n",fLcCuts[4]);
1365     printf("    dist12    [cm]  < %f\n",fLcCuts[5]);
1366     printf("    sigmavert [cm]   < %f\n",fLcCuts[6]);
1367     printf("    dist prim-sec [cm] > %f\n",fLcCuts[7]);
1368     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fLcCuts[8]);
1369     printf("    cosThetaPoint    > %f\n",fLcCuts[9]);
1370     printf("    Sum d0^2 [cm^2]  > %f\n",fLcCuts[10]);
1371     printf("    dca cut [cm]  < %f\n",fLcCuts[11]);
1372   }
1373
1374   return;
1375 }
1376 //-----------------------------------------------------------------------------
1377 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
1378                                                                  Double_t &dispersion,Bool_t useTRefArray) const
1379 {
1380   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1381
1382   AliESDVertex *vertexESD = 0;
1383   AliAODVertex *vertexAOD = 0;
1384
1385   if(!fSecVtxWithKF) { // AliVertexerTracks
1386
1387     AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1388     vertexer->SetVtxStart(fV1);
1389     vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1390     delete vertexer; vertexer=NULL;
1391
1392     if(!vertexESD) return vertexAOD;
1393
1394     if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
1395       AliDebug(2,"vertexing failed");
1396       delete vertexESD; vertexESD=NULL;
1397       return vertexAOD;
1398     }
1399
1400   } else { // Kalman Filter vertexer (AliKFParticle)
1401
1402     AliKFParticle::SetField(fBzkG);
1403
1404     AliKFVertex vertexKF;
1405
1406     Int_t nTrks = trkArray->GetEntriesFast();
1407     for(Int_t i=0; i<nTrks; i++) {
1408       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1409       AliKFParticle daughterKF(*esdTrack,211);
1410       vertexKF.AddDaughter(daughterKF);
1411     }
1412     vertexESD = new AliESDVertex(vertexKF.Parameters(),
1413                                  vertexKF.CovarianceMatrix(),
1414                                  vertexKF.GetChi2(),
1415                                  vertexKF.GetNContributors());
1416
1417   }
1418
1419   // convert to AliAODVertex
1420   Double_t pos[3],cov[6],chi2perNDF;
1421   vertexESD->GetXYZ(pos); // position
1422   vertexESD->GetCovMatrix(cov); //covariance matrix
1423   chi2perNDF = vertexESD->GetChi2toNDF();
1424   dispersion = vertexESD->GetDispersion();
1425   delete vertexESD; vertexESD=NULL;
1426
1427   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
1428   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
1429
1430   return vertexAOD;
1431 }
1432 //-----------------------------------------------------------------------------
1433 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1434                                              Int_t nprongs,
1435                                              Double_t *px,
1436                                              Double_t *py,
1437                                              Double_t *pz) const {
1438   // Check invariant mass cut
1439
1440   Short_t dummycharge=0;
1441   Double_t *dummyd0 = new Double_t[nprongs];
1442   for(Int_t ip=0;ip<nprongs;ip++) dummyd0[ip]=0.;
1443   AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
1444   delete [] dummyd0; dummyd0=NULL;
1445
1446   UInt_t pdg2[2],pdg3[3],pdg4[4];
1447   Double_t mPDG,minv;
1448
1449   Bool_t retval=kFALSE;
1450   switch (decay)
1451     {
1452     case 0:                  // D0->Kpi
1453       pdg2[0]=211; pdg2[1]=321;
1454       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1455       minv = rd->InvMass(nprongs,pdg2);
1456       if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1457       pdg2[0]=321; pdg2[1]=211;
1458       minv = rd->InvMass(nprongs,pdg2);
1459       if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1460       break;
1461     case 1:                  // JPSI->ee
1462       pdg2[0]=11; pdg2[1]=11;
1463       mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
1464       minv = rd->InvMass(nprongs,pdg2);
1465       if(TMath::Abs(minv-mPDG)<fBtoJPSICuts[0]) retval=kTRUE;
1466       break;
1467     case 2:                  // D+->Kpipi
1468       pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
1469       mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
1470       minv = rd->InvMass(nprongs,pdg3);
1471       if(TMath::Abs(minv-mPDG)<fDplusCuts[0]) retval=kTRUE;
1472                             // Ds+->KKpi
1473       pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
1474       mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1475       minv = rd->InvMass(nprongs,pdg3);
1476       if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1477       pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
1478       minv = rd->InvMass(nprongs,pdg3);
1479       if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1480                             // Lc->pKpi
1481       pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
1482       mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1483       minv = rd->InvMass(nprongs,pdg3);
1484       if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1485       pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
1486       minv = rd->InvMass(nprongs,pdg3);
1487       if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1488       break;
1489     case 3:                  // D*->D0pi
1490       pdg2[0]=421; pdg2[1]=211;
1491       mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
1492       minv = rd->InvMass(nprongs,pdg2);
1493       if(TMath::Abs(minv-mPDG)<fDstarCuts[0]) retval=kTRUE;
1494       break;
1495     case 4:                 // D0->K3pi without PID
1496       pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
1497       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1498       minv = rd->InvMass(nprongs,pdg4);
1499       if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0])
1500         retval=kTRUE;
1501       pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
1502       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1503       minv = rd->InvMass(nprongs,pdg4);
1504       if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1505       pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
1506       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1507       minv = rd->InvMass(nprongs,pdg4);
1508       if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1509       pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
1510       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1511       minv = rd->InvMass(nprongs,pdg4);
1512       if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1513       break;
1514     default:
1515       printf("SelectInvMass(): wrong decay selection\n");
1516       break;
1517     }
1518
1519   delete rd; rd=NULL;
1520
1521   return retval;
1522 }
1523 //-----------------------------------------------------------------------------
1524 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(AliVEvent *event,
1525                                    TObjArray &seleTrksArray,Int_t &nSeleTrks,
1526                                                        UChar_t *seleFlags)
1527 {
1528   // Apply single-track preselection.
1529   // Fill a TObjArray with selected tracks (for displaced vertices or
1530   // soft pion from D*). Selection flag stored in seleFlags.
1531   // Create the AliESDVertex object (convert from AliAODVertex if necessary)
1532   // In case of AOD input, also fill fAODMap for track index<->ID
1533
1534   const AliVVertex *vprimary = event->GetPrimaryVertex();
1535
1536   if(fV1) { delete fV1; fV1=NULL; }
1537
1538   Int_t nindices=0;
1539   UShort_t *indices = 0;
1540   Double_t pos[3],cov[6];
1541
1542   if(!fInputAOD) { // ESD
1543     fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
1544   } else {         // AOD
1545     vprimary->GetXYZ(pos);
1546     vprimary->GetCovarianceMatrix(cov);
1547     fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1548     indices = new UShort_t[event->GetNumberOfTracks()];
1549     for(Int_t j=0; j<100000; j++) fAODMap[j]=-1;
1550   }
1551
1552
1553   Int_t entries = (Int_t)event->GetNumberOfTracks();
1554   Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
1555   nSeleTrks=0;
1556
1557   // transfer ITS tracks from event to arrays
1558   for(Int_t i=0; i<entries; i++) {
1559     AliVTrack *track = (AliVTrack*)event->GetTrack(i);
1560
1561     // TEMPORARY: check that the cov matrix is there
1562     Double_t covtest[21];
1563     if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1564     //
1565
1566     if(fInputAOD) {
1567       AliAODTrack *aodt = (AliAODTrack*)track;
1568       if(aodt->GetUsedForPrimVtxFit()) {
1569         indices[nindices]=aodt->GetID(); nindices++;
1570       }
1571       fAODMap[(Int_t)aodt->GetID()] = i;
1572     }
1573
1574     AliESDtrack *esdt = 0;
1575
1576     if(!fInputAOD) {
1577       esdt = (AliESDtrack*)track;
1578     } else {
1579       esdt = new AliESDtrack(track);
1580     }
1581
1582     // single track selection
1583     okDisplaced=kFALSE; okSoftPi=kFALSE;
1584     if(SingleTrkCuts(esdt,okDisplaced,okSoftPi)) {
1585       seleTrksArray.AddLast(esdt);
1586       seleFlags[nSeleTrks]=0;
1587       if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
1588       if(okSoftPi)    SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
1589       nSeleTrks++;
1590     } else {
1591       if(fInputAOD) delete esdt;
1592       esdt = NULL;
1593       continue;
1594     }
1595
1596   } // end loop on tracks
1597
1598   // primary vertex from AOD
1599   if(fInputAOD) {
1600     delete fV1; fV1=NULL;
1601     vprimary->GetXYZ(pos);
1602     vprimary->GetCovarianceMatrix(cov);
1603     Double_t chi2toNDF = vprimary->GetChi2perNDF();
1604     Int_t ncontr=nindices;
1605     if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
1606     Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
1607     fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
1608     fV1->SetTitle(vprimary->GetTitle());
1609     fV1->SetIndices(nindices,indices);
1610   }
1611   if(indices) { delete [] indices; indices=NULL; }
1612
1613
1614   return;
1615 }
1616 //-----------------------------------------------------------------------------
1617 void AliAnalysisVertexingHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
1618                                    Double_t cut2,Double_t cut3,Double_t cut4,
1619                                    Double_t cut5,Double_t cut6,
1620                                    Double_t cut7,Double_t cut8)
1621 {
1622   // Set the cuts for D0 selection
1623   fD0toKpiCuts[0] = cut0;
1624   fD0toKpiCuts[1] = cut1;
1625   fD0toKpiCuts[2] = cut2;
1626   fD0toKpiCuts[3] = cut3;
1627   fD0toKpiCuts[4] = cut4;
1628   fD0toKpiCuts[5] = cut5;
1629   fD0toKpiCuts[6] = cut6;
1630   fD0toKpiCuts[7] = cut7;
1631   fD0toKpiCuts[8] = cut8;
1632
1633   return;
1634 }
1635 //-----------------------------------------------------------------------------
1636 void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9])
1637 {
1638   // Set the cuts for D0 selection
1639
1640   for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
1641
1642   return;
1643 }
1644 //-----------------------------------------------------------------------------
1645 void AliAnalysisVertexingHF::SetD0fromDstarCuts(Double_t cut0,Double_t cut1,
1646                                    Double_t cut2,Double_t cut3,Double_t cut4,
1647                                    Double_t cut5,Double_t cut6,
1648                                    Double_t cut7,Double_t cut8)
1649 {
1650   // Set the cuts for D0 from D* selection
1651   fD0fromDstarCuts[0] = cut0;
1652   fD0fromDstarCuts[1] = cut1;
1653   fD0fromDstarCuts[2] = cut2;
1654   fD0fromDstarCuts[3] = cut3;
1655   fD0fromDstarCuts[4] = cut4;
1656   fD0fromDstarCuts[5] = cut5;
1657   fD0fromDstarCuts[6] = cut6;
1658   fD0fromDstarCuts[7] = cut7;
1659   fD0fromDstarCuts[8] = cut8;
1660
1661   return;
1662 }
1663 //-----------------------------------------------------------------------------
1664 void AliAnalysisVertexingHF::SetD0fromDstarCuts(const Double_t cuts[9])
1665 {
1666   // Set the cuts for D0 from D* selection
1667
1668   for(Int_t i=0; i<9; i++) fD0fromDstarCuts[i] = cuts[i];
1669
1670   return;
1671 }
1672 //-----------------------------------------------------------------------------
1673 void AliAnalysisVertexingHF::SetDstarCuts(Double_t cut0,Double_t cut1,
1674                                           Double_t cut2,Double_t cut3,
1675                                           Double_t cut4)
1676 {
1677   // Set the cuts for D* selection
1678   fDstarCuts[0] = cut0;
1679   fDstarCuts[1] = cut1;
1680   fDstarCuts[2] = cut2;
1681   fDstarCuts[3] = cut3;
1682   fDstarCuts[4] = cut4;
1683
1684   return;
1685 }
1686 //-----------------------------------------------------------------------------
1687 void AliAnalysisVertexingHF::SetDstarCuts(const Double_t cuts[5])
1688 {
1689   // Set the cuts for D* selection
1690
1691   for(Int_t i=0; i<5; i++) fDstarCuts[i] = cuts[i];
1692
1693   return;
1694 }
1695 //-----------------------------------------------------------------------------
1696 void AliAnalysisVertexingHF::SetBtoJPSICuts(Double_t cut0,Double_t cut1,
1697                                    Double_t cut2,Double_t cut3,Double_t cut4,
1698                                    Double_t cut5,Double_t cut6,
1699                                    Double_t cut7,Double_t cut8)
1700 {
1701   // Set the cuts for J/psi from B selection
1702   fBtoJPSICuts[0] = cut0;
1703   fBtoJPSICuts[1] = cut1;
1704   fBtoJPSICuts[2] = cut2;
1705   fBtoJPSICuts[3] = cut3;
1706   fBtoJPSICuts[4] = cut4;
1707   fBtoJPSICuts[5] = cut5;
1708   fBtoJPSICuts[6] = cut6;
1709   fBtoJPSICuts[7] = cut7;
1710   fBtoJPSICuts[8] = cut8;
1711
1712   return;
1713 }
1714 //-----------------------------------------------------------------------------
1715 void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9])
1716 {
1717   // Set the cuts for J/psi from B selection
1718
1719   for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i];
1720
1721   return;
1722 }
1723 //-----------------------------------------------------------------------------
1724 void AliAnalysisVertexingHF::SetDplusCuts(Double_t cut0,Double_t cut1,
1725                                    Double_t cut2,Double_t cut3,Double_t cut4,
1726                                    Double_t cut5,Double_t cut6,
1727                                    Double_t cut7,Double_t cut8,
1728                                    Double_t cut9,Double_t cut10,Double_t cut11)
1729 {
1730   // Set the cuts for Dplus->Kpipi selection
1731   fDplusCuts[0] = cut0;
1732   fDplusCuts[1] = cut1;
1733   fDplusCuts[2] = cut2;
1734   fDplusCuts[3] = cut3;
1735   fDplusCuts[4] = cut4;
1736   fDplusCuts[5] = cut5;
1737   fDplusCuts[6] = cut6;
1738   fDplusCuts[7] = cut7;
1739   fDplusCuts[8] = cut8;
1740   fDplusCuts[9] = cut9;
1741   fDplusCuts[10] = cut10;
1742   fDplusCuts[11] = cut11;
1743
1744   return;
1745 }
1746 //-----------------------------------------------------------------------------
1747 void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12])
1748 {
1749   // Set the cuts for Dplus->Kpipi selection
1750
1751   for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i];
1752
1753   return;
1754 }
1755 //-----------------------------------------------------------------------------
1756 void AliAnalysisVertexingHF::SetDsCuts(Double_t cut0,Double_t cut1,
1757                                    Double_t cut2,Double_t cut3,Double_t cut4,
1758                                    Double_t cut5,Double_t cut6,
1759                                    Double_t cut7,Double_t cut8,
1760                                    Double_t cut9,Double_t cut10,
1761                                    Double_t cut11,Double_t cut12)
1762 {
1763   // Set the cuts for Ds->KKpi selection
1764   fDsCuts[0] = cut0;
1765   fDsCuts[1] = cut1;
1766   fDsCuts[2] = cut2;
1767   fDsCuts[3] = cut3;
1768   fDsCuts[4] = cut4;
1769   fDsCuts[5] = cut5;
1770   fDsCuts[6] = cut6;
1771   fDsCuts[7] = cut7;
1772   fDsCuts[8] = cut8;
1773   fDsCuts[9] = cut9;
1774   fDsCuts[10] = cut10;
1775   fDsCuts[11] = cut11;
1776   fDsCuts[12] = cut12;
1777
1778   return;
1779 }
1780 //-----------------------------------------------------------------------------
1781 void AliAnalysisVertexingHF::SetDsCuts(const Double_t cuts[13])
1782 {
1783   // Set the cuts for Ds->KKpi selection
1784
1785   for(Int_t i=0; i<13; i++) fDsCuts[i] = cuts[i];
1786
1787   return;
1788 }
1789 //-----------------------------------------------------------------------------
1790 void AliAnalysisVertexingHF::SetLcCuts(Double_t cut0,Double_t cut1,
1791                                    Double_t cut2,Double_t cut3,Double_t cut4,
1792                                    Double_t cut5,Double_t cut6,
1793                                    Double_t cut7,Double_t cut8,
1794                                    Double_t cut9,Double_t cut10,Double_t cut11)
1795 {
1796   // Set the cuts for Lc->pKpi selection
1797   fLcCuts[0] = cut0;
1798   fLcCuts[1] = cut1;
1799   fLcCuts[2] = cut2;
1800   fLcCuts[3] = cut3;
1801   fLcCuts[4] = cut4;
1802   fLcCuts[5] = cut5;
1803   fLcCuts[6] = cut6;
1804   fLcCuts[7] = cut7;
1805   fLcCuts[8] = cut8;
1806   fLcCuts[9] = cut9;
1807   fLcCuts[10] = cut10;
1808   fLcCuts[11] = cut11;
1809
1810   return;
1811 }
1812 //-----------------------------------------------------------------------------
1813 void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[12])
1814 {
1815   // Set the cuts for Lc->pKpi selection
1816
1817   for(Int_t i=0; i<12; i++) fLcCuts[i] = cuts[i];
1818
1819   return;
1820 }
1821 //-----------------------------------------------------------------------------
1822 void AliAnalysisVertexingHF::SetD0to4ProngsCuts(Double_t cut0,Double_t cut1,
1823                                    Double_t cut2,Double_t cut3,Double_t cut4,
1824                                    Double_t cut5,Double_t cut6,
1825                                    Double_t cut7,Double_t cut8)
1826 {
1827   // Set the cuts for D0->Kpipipi selection
1828
1829   fD0to4ProngsCuts[0] = cut0;
1830   fD0to4ProngsCuts[1] = cut1;
1831   fD0to4ProngsCuts[2] = cut2;
1832   fD0to4ProngsCuts[3] = cut3;
1833   fD0to4ProngsCuts[4] = cut4;
1834   fD0to4ProngsCuts[5] = cut5;
1835   fD0to4ProngsCuts[6] = cut6;
1836   fD0to4ProngsCuts[7] = cut7;
1837   fD0to4ProngsCuts[8] = cut8;
1838
1839   return;
1840 }
1841 //-----------------------------------------------------------------------------
1842 void AliAnalysisVertexingHF::SetD0to4ProngsCuts(const Double_t cuts[9])
1843 {
1844   // Set the cuts for D0->Kpipipi selection
1845
1846   for(Int_t i=0; i<9; i++) fD0to4ProngsCuts[i] = cuts[i];
1847
1848   return;
1849 }
1850 //-----------------------------------------------------------------------------
1851 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
1852                                              Bool_t &okDisplaced,Bool_t &okSoftPi) const
1853 {
1854   // Check if track passes some kinematical cuts
1855
1856   // this is needed to store the impact parameters
1857   trk->RelateToVertex(fV1,fBzkG,kVeryBig);
1858
1859   UInt_t selectInfo;
1860   //
1861   // Track selection, displaced tracks
1862   selectInfo = 0;
1863   if(fTrackFilter) {
1864     selectInfo = fTrackFilter->IsSelected(trk);
1865   }
1866   if(selectInfo) okDisplaced=kTRUE;
1867   // Track selection, soft pions
1868   selectInfo = 0;
1869   if(fDstar && fTrackFilterSoftPi) {
1870     selectInfo = fTrackFilterSoftPi->IsSelected(trk);
1871   }
1872   if(selectInfo) okSoftPi=kTRUE;
1873
1874   if(okDisplaced || okSoftPi) return kTRUE;
1875
1876   return kFALSE;
1877 }
1878 //-----------------------------------------------------------------------------