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