PWG3/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 store as objects deriving from AliAODRecoDecay.
  19 // An example of usage can be found in the macro AliVertexingHFTest.C.
  20 // Can be used as a task of AliAnalysisManager by means of the interface
  21 // classes AliAnalysisTaskSEVertexingHF or AliAnalysisTaskVertexingHF.
  22 //
  23 //  Origin: E.Bruna, G.E.Bruno, A.Dainese, F.Prino, R.Romita
  24 //  Contact: andrea.dainese@lnl.infn.it
  25 //----------------------------------------------------------------------------
  26 #include <TTree.h>
  27 #include <TFile.h>
  28 #include <TDatabasePDG.h>
  29 #include <TString.h>
  30 #include "AliESDEvent.h"
  31 #include "AliVertexerTracks.h"
  32 #include "AliKFParticle.h"
  33 #include "AliESDVertex.h"
  34 #include "AliAODEvent.h"
  35 #include "AliAODRecoDecay.h"
  36 #include "AliAODRecoDecayHF.h"
  37 #include "AliAODRecoDecayHF2Prong.h"
  38 #include "AliAODRecoDecayHF3Prong.h"
  39 #include "AliAODRecoDecayHF4Prong.h"
  40 #include "AliAnalysisVertexingHF.h"
  41
  42 ClassImp(AliAnalysisVertexingHF)
  43
  44 //----------------------------------------------------------------------------
  45 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
  46 fBzkG(0.),
  47 fSecVtxWithKF(kFALSE),
  48 fUseTRef(kFALSE),
  49 fRecoPrimVtxSkippingTrks(kFALSE),
  50 fRmTrksFromPrimVtx(kFALSE),
  51 fV1(0x0),
  52 fDebug(0),
  53 fD0toKpi(kTRUE),
  54 fJPSItoEle(kTRUE),
  55 f3Prong(kTRUE),
  56 f4Prong(kTRUE),
  57 fITSrefit(kFALSE),
  58 fBothSPD(kTRUE),
  59 fMinITSCls(5),
  60 fMinPtCut(0.),
  61 fMind0rphiCut(0.)
  62 {
  63   // Default constructor
  64
  65   SetD0toKpiCuts();
  66   SetBtoJPSICuts();
  67   SetDplusCuts();
  68   SetDsCuts();
  69   SetLcCuts();
  70 }
  71 //--------------------------------------------------------------------------
  72 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
  73 TNamed(source),
  74 fBzkG(source.fBzkG),
  75 fSecVtxWithKF(source.fSecVtxWithKF),
  76 fUseTRef(source.fUseTRef),
  77 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
  78 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
  79 fV1(source.fV1),
  80 fDebug(source.fDebug),
  81 fD0toKpi(source.fD0toKpi),
  82 fJPSItoEle(source.fJPSItoEle),
  83 f3Prong(source.f3Prong),
  84 f4Prong(source.f4Prong),
  85 fITSrefit(source.fITSrefit),
  86 fBothSPD(source.fBothSPD),
  87 fMinITSCls(source.fMinITSCls),
  88 fMinPtCut(source.fMinPtCut),
  89 fMind0rphiCut(source.fMind0rphiCut)
  90 {
  91   //
  92   // Copy constructor
  93   //
  94   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
  95   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
  96   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
  97   for(Int_t i=0; i<1; i++)  fDsCuts[i]=source.fDsCuts[i];
  98   for(Int_t i=0; i<1; i++)  fLcCuts[i]=source.fLcCuts[i];
  99 }
 100 //--------------------------------------------------------------------------
 101 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
 102 {
 103   //
 104   // assignment operator
 105   //
 106   if(&source == this) return *this;
 107   fBzkG = source.fBzkG;
 108   fSecVtxWithKF = source.fSecVtxWithKF;
 109   fUseTRef = source.fUseTRef;
 110   fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
 111   fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
 112   fV1 = source.fV1;
 113   fDebug = source.fDebug;
 114   fD0toKpi = source.fD0toKpi;
 115   fJPSItoEle = source.fJPSItoEle;
 116   f3Prong = source.f3Prong;
 117   f4Prong = source.f4Prong;
 118   fITSrefit = source.fITSrefit;
 119   fBothSPD = source.fBothSPD;
 120   fMinITSCls = source.fMinITSCls;
 121   fMinPtCut = source.fMinPtCut;
 122   fMind0rphiCut = source.fMind0rphiCut;
 123
 124   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
 125   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
 126   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
 127   for(Int_t i=0; i<1; i++)  fDsCuts[i]=source.fDsCuts[i];
 128   for(Int_t i=0; i<1; i++)  fLcCuts[i]=source.fLcCuts[i];
 129
 130   return *this;
 131 }
 132 //----------------------------------------------------------------------------
 133 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
 134   // Destructor
 135   if(fV1) { delete fV1; fV1=0; }
 136 }
 137 //----------------------------------------------------------------------------
 138 void AliAnalysisVertexingHF::FindCandidatesESDtoAOD(AliESDEvent *esd,
 139                                         TClonesArray *aodVerticesHFTClArr,
 140                                         TClonesArray *aodD0toKpiTClArr,
 141                                         TClonesArray *aodJPSItoEleTClArr,
 142                                         TClonesArray *aodCharm3ProngTClArr,
 143                                         TClonesArray *aodCharm4ProngTClArr)
 144 {
 145   // Find heavy-flavour vertex candidates
 146   // Input:  ESD
 147   // Output: AOD (additional branches added)
 148
 149   if(!aodVerticesHFTClArr) {
 150     printf("ERROR: no aodVerticesHFTClArr");
 151     return;
 152   }
 153   if(fD0toKpi && !aodD0toKpiTClArr) {
 154     printf("ERROR: no aodD0toKpiTClArr");
 155     return;
 156   }
 157   if(fJPSItoEle && !aodJPSItoEleTClArr) {
 158     printf("ERROR: no aodJPSItoEleTClArr");
 159     return;
 160   }
 161   if(f3Prong && !aodCharm3ProngTClArr) {
 162     printf("ERROR: no aodCharm3ProngTClArr");
 163     return;
 164   }
 165   if(f4Prong && !aodCharm4ProngTClArr) {
 166     printf("ERROR: no aodCharm4ProngTClArr");
 167     return;
 168   }
 169
 170   // delete candidates from previous event and create references
 171   Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0;
 172   aodVerticesHFTClArr->Delete();
 173   iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
 174   TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
 175   if(fD0toKpi)   {
 176     aodD0toKpiTClArr->Delete();
 177     iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
 178   }
 179   if(fJPSItoEle) {
 180     aodJPSItoEleTClArr->Delete();
 181     iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
 182   }
 183   if(f3Prong) {
 184     aodCharm3ProngTClArr->Delete();
 185     i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
 186   }
 187   if(f4Prong) {
 188     aodCharm4ProngTClArr->Delete();
 189     i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
 190   }
 191   TClonesArray &aodD0toKpiRef     = *aodD0toKpiTClArr;
 192   TClonesArray &aodJPSItoEleRef   = *aodJPSItoEleTClArr;
 193   TClonesArray &aodCharm3ProngRef = *aodCharm3ProngTClArr;
 194   TClonesArray &aodCharm4ProngRef = *aodCharm4ProngTClArr;
 195
 196
 197   AliAODRecoDecayHF2Prong *io2Prong = 0;
 198   AliAODRecoDecayHF3Prong *io3Prong = 0;
 199   AliAODRecoDecayHF4Prong *io4Prong = 0;
 200
 201   Int_t    iTrkP1,iTrkP2,iTrkN1,iTrkN2,trkEntries;
 202   Int_t    nTrksP=0,nTrksN=0;
 203   Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2;
 204   Bool_t   okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
 205   AliESDtrack *postrack1 = 0;
 206   AliESDtrack *postrack2 = 0;
 207   AliESDtrack *negtrack1 = 0;
 208   AliESDtrack *negtrack2 = 0;
 209   Double_t dcaMax = fD0toKpiCuts[1];
 210   if(dcaMax < fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1];
 211   if(dcaMax < fDplusCuts[11])  dcaMax=fDplusCuts[11];
 212   if(fDebug) printf(" dca cut set to %f cm\n",dcaMax);
 213
 214   //------- SINGLE EVENT ANALYSIS ---------------------------------
 215   //
 216   Int_t ev = (Int_t)esd->GetEventNumberInFile();
 217   printf("--- Finding candidates in event %d\n",ev);
 218
 219
 220   // get Bz from ESD
 221   fBzkG = (Double_t)esd->GetMagneticField();
 222
 223   trkEntries = (Int_t)esd->GetNumberOfTracks();
 224   printf(" Number of tracks: %d\n",trkEntries);
 225
 226   if(trkEntries<2 || !esd->GetPrimaryVertex()) {
 227     return;
 228   }
 229
 230   // retrieve primary vertex from the AliESDEvent
 231   AliESDVertex copy(*(esd->GetPrimaryVertex()));
 232   SetPrimaryVertex(&copy);
 233
 234   // call function which applies sigle-track selection and
 235   // separates positives and negatives
 236   TObjArray trksP(trkEntries/2);
 237   TObjArray trksN(trkEntries/2);
 238   SelectTracks(esd,trksP,nTrksP,trksN,nTrksN);
 239
 240   printf(" Pos. tracks: %d    Neg. tracks: %d\n",nTrksP,nTrksN);
 241
 242   TObjArray *twoTrackArray1  = new TObjArray(2);
 243   TObjArray *twoTrackArray2  = new TObjArray(2);
 244   TObjArray *threeTrackArray = new TObjArray(3);
 245   TObjArray *fourTrackArray  = new TObjArray(4);
 246
 247   // LOOP ON  POSITIVE  TRACKS
 248   for(iTrkP1=0; iTrkP1<nTrksP; iTrkP1++) {
 249     if(fDebug && iTrkP1%1==0) printf("  Processing positive track number %d of %d\n",iTrkP1,nTrksP);
 250     // get track from tracks array
 251     postrack1 = (AliESDtrack*)trksP.UncheckedAt(iTrkP1);
 252
 253     // LOOP ON  NEGATIVE  TRACKS
 254     for(iTrkN1=0; iTrkN1<nTrksN; iTrkN1++) {
 255       if(fDebug && iTrkN1%1==0) printf("    Processing negative track number %d of %d\n",iTrkN1,nTrksN);
 256       // get track from tracks array
 257       negtrack1 = (AliESDtrack*)trksN.UncheckedAt(iTrkN1);
 258       // DCA between the two tracks
 259       dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 260       if(dcap1n1>dcaMax) { negtrack1=0; continue; }
 261       // Vertexing
 262       twoTrackArray1->AddAt(postrack1,0);
 263       twoTrackArray1->AddAt(negtrack1,1);
 264       AliESDVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1);
 265       if(!vertexp1n1) {
 266         twoTrackArray1->Clear();
 267         negtrack1=0;
 268         continue;
 269       }
 270       if(fD0toKpi || fJPSItoEle) {
 271         io2Prong = Make2Prong(twoTrackArray1,esd,vertexp1n1,dcap1n1,okD0,okJPSI);
 272         if(okD0 || okJPSI) {
 273           if(fUseTRef) {
 274             AliAODVertex *v = new(verticesHFRef[iVerticesHF++])
 275               AliAODVertex(*(io2Prong->GetOwnSecondaryVtx()));
 276             v->SetType(AliAODVertex::kUndef); // to be changed
 277             Double_t px[2]={io2Prong->PxProng(0),io2Prong->PxProng(1)};
 278             Double_t py[2]={io2Prong->PyProng(0),io2Prong->PyProng(1)};
 279             Double_t pz[2]={io2Prong->PzProng(0),io2Prong->PzProng(1)};
 280             Double_t d0[2]={io2Prong->Getd0Prong(0),io2Prong->Getd0Prong(1)};
 281             Double_t d0err[2]={io2Prong->Getd0errProng(0),io2Prong->Getd0errProng(1)};
 282             UShort_t id[2]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID()};
 283             if(okD0) {
 284               AliAODRecoDecayHF2Prong *rd=new(aodD0toKpiRef[iD0toKpi++])
 285                 AliAODRecoDecayHF2Prong(v,px,py,pz,d0,d0err,dcap1n1);
 286               if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) rd->SetOwnPrimaryVtx(io2Prong->GetOwnPrimaryVtx());
 287               rd->SetProngIDs(2,id);
 288               v->SetParent(rd);
 289             }
 290             if(okJPSI) {
 291               AliAODRecoDecayHF2Prong *rd=new(aodJPSItoEleRef[iJPSItoEle++])
 292                 AliAODRecoDecayHF2Prong(v,px,py,pz,d0,d0err,dcap1n1);
 293               if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) rd->SetOwnPrimaryVtx(io2Prong->GetOwnPrimaryVtx());
 294               rd->SetProngIDs(2,id);
 295               if(!okD0) v->SetParent(rd); // do something better here...
 296             }
 297             //printf("DCA: %f\n",rd->GetDCA());
 298           } else {
 299             if(okD0)   new(aodD0toKpiRef[iD0toKpi++]) AliAODRecoDecayHF2Prong(*io2Prong);
 300             if(okJPSI) new(aodJPSItoEleRef[iJPSItoEle++]) AliAODRecoDecayHF2Prong(*io2Prong);
 301           }
 302         }
 303         //delete io2Prong;
 304         io2Prong=NULL;
 305       }
 306
 307       twoTrackArray1->Clear();
 308       if(!f3Prong && !f4Prong)  {
 309         negtrack1=0;
 310         delete vertexp1n1;
 311         continue;
 312       }
 313
 314
 315       // 2nd LOOP  ON  POSITIVE  TRACKS
 316       for(iTrkP2=iTrkP1+1; iTrkP2<nTrksP; iTrkP2++) {
 317         // get track from tracks array
 318         postrack2 = (AliESDtrack*)trksP.UncheckedAt(iTrkP2);
 319         dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 320         if(dcap2n1>dcaMax) { postrack2=0; continue; }
 321         dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
 322         if(dcap1p2>dcaMax) { postrack2=0; continue; }
 323
 324         // Vertexing
 325         twoTrackArray2->AddAt(postrack2,0);
 326         twoTrackArray2->AddAt(negtrack1,1);
 327         AliESDVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2);
 328         if(!vertexp2n1) {
 329           twoTrackArray2->Clear();
 330           postrack2=0;
 331           continue;
 332         }
 333         if(f3Prong) {
 334           threeTrackArray->AddAt(postrack1,0);
 335           threeTrackArray->AddAt(negtrack1,1);
 336           threeTrackArray->AddAt(postrack2,2);
 337           io3Prong = Make3Prong(threeTrackArray,esd,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
 338           if(ok3Prong) {
 339             if(fUseTRef) {
 340               AliAODVertex *v = new(verticesHFRef[iVerticesHF++])
 341                 AliAODVertex(*(io3Prong->GetOwnSecondaryVtx()));
 342               v->SetType(AliAODVertex::kUndef);
 343               Double_t px[3]={io3Prong->PxProng(0),io3Prong->PxProng(1),io3Prong->PxProng(2)};
 344               Double_t py[3]={io3Prong->PyProng(0),io3Prong->PyProng(1),io3Prong->PyProng(2)};
 345               Double_t pz[3]={io3Prong->PzProng(0),io3Prong->PzProng(1),io3Prong->PzProng(2)};
 346               Double_t d0[3]={io3Prong->Getd0Prong(0),io3Prong->Getd0Prong(1),io3Prong->Getd0Prong(2)};
 347               Double_t d0err[3]={io3Prong->Getd0errProng(0),io3Prong->Getd0errProng(1),io3Prong->Getd0errProng(2)};
 348               Double_t dcas[3]={io3Prong->GetDCA(0),io3Prong->GetDCA(1),io3Prong->GetDCA(2)};
 349               UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID()};
 350               AliAODRecoDecayHF3Prong *rd=new(aodCharm3ProngRef[i3Prong++])
 351                 AliAODRecoDecayHF3Prong(v,px,py,pz,d0,d0err,dcas,io3Prong->GetSigmaVert(),io3Prong->GetDist12toPrim(),io3Prong->GetDist23toPrim(),io3Prong->GetCharge());
 352               if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) rd->SetOwnPrimaryVtx(io3Prong->GetOwnPrimaryVtx());
 353               rd->SetProngIDs(3,id);
 354               v->SetParent(rd);
 355             } else {
 356               new(aodD0toKpiRef[i3Prong++]) AliAODRecoDecayHF3Prong(*io3Prong);
 357             }
 358           }
 359           if(io3Prong) { /*delete io3Prong;*/ io3Prong=NULL; }
 360         }
 361         if(f4Prong) {
 362           // 3rd LOOP  ON  NEGATIVE  TRACKS (for 4 prong)
 363           for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) {
 364             // get track from tracks array
 365             negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2);
 366             dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 367             if(dcap1n2>dcaMax) { negtrack2=0; continue; }
 368             // Vertexing
 369             fourTrackArray->AddAt(postrack1,0);
 370             fourTrackArray->AddAt(negtrack1,1);
 371             fourTrackArray->AddAt(postrack2,2);
 372             fourTrackArray->AddAt(negtrack2,3);
 373             io4Prong = Make4Prong(fourTrackArray,esd,vertexp1n1,vertexp2n1,dcap1n1,dcap1n2,dcap2n1,ok4Prong);
 374             if(ok4Prong) {
 375               if(fUseTRef) {
 376                 AliAODVertex *v = new(verticesHFRef[iVerticesHF++])
 377                   AliAODVertex(*(io4Prong->GetOwnSecondaryVtx()));
 378                 v->SetType(AliAODVertex::kUndef);
 379                 Double_t px[4]={io4Prong->PxProng(0),io4Prong->PxProng(1),
 380                                 io4Prong->PxProng(2),io4Prong->PxProng(3)};
 381                 Double_t py[4]={io4Prong->PyProng(0),io4Prong->PyProng(1),
 382                                 io4Prong->PyProng(2),io4Prong->PyProng(3)};
 383                 Double_t pz[4]={io4Prong->PzProng(0),io4Prong->PzProng(1),
 384                                 io4Prong->PzProng(2),io4Prong->PzProng(3)};
 385                 Double_t d0[4]={io4Prong->Getd0Prong(0),io4Prong->Getd0Prong(1),
 386                                 io4Prong->Getd0Prong(2),io4Prong->Getd0Prong(3)};
 387                 Double_t d0err[4]={io4Prong->Getd0errProng(0),io4Prong->Getd0errProng(1),
 388                                    io4Prong->Getd0errProng(2),io4Prong->Getd0errProng(3)};
 389                 Double_t dcas[6]; io4Prong->GetDCAs(dcas);
 390                 UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
 391                 AliAODRecoDecayHF4Prong *rd=new(aodCharm4ProngRef[i4Prong++])
 392                   AliAODRecoDecayHF4Prong(v,px,py,pz,d0,d0err,dcas,io4Prong->GetDist12toPrim(),io4Prong->GetDist23toPrim(),io4Prong->GetDist14toPrim(),io4Prong->GetDist34toPrim(),io4Prong->GetCharge());
 393                 if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) rd->SetOwnPrimaryVtx(io4Prong->GetOwnPrimaryVtx());
 394                 rd->SetProngIDs(4,id);
 395                 v->SetParent(rd);
 396               } else {
 397                 new(aodD0toKpiRef[i4Prong++]) AliAODRecoDecayHF4Prong(*io4Prong);
 398               }
 399             }
 400             if(io4Prong) { /*delete io4Prong;*/ io4Prong=NULL; }
 401             fourTrackArray->Clear();
 402             negtrack2 = 0;
 403           } // end loop on negative tracks
 404         }
 405         postrack2 = 0;
 406         delete vertexp2n1;
 407       } // end 2nd loop on positive tracks
 408       twoTrackArray2->Clear();
 409
 410       // 2nd LOOP  ON  NEGATIVE  TRACKS
 411       for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) {
 412         // get track from tracks array
 413         negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2);
 414         dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 415         if(dcap1n2>dcaMax) { negtrack2=0; continue; }
 416         dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 417         if(dcan1n2>dcaMax) { negtrack2=0; continue; }
 418
 419         // Vertexing
 420         twoTrackArray2->AddAt(postrack1,0);
 421         twoTrackArray2->AddAt(negtrack2,1);
 422         AliESDVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2);
 423         if(!vertexp1n2) {
 424           twoTrackArray2->Clear();
 425           negtrack2=0;
 426           continue;
 427         }
 428         if(f3Prong) {
 429           threeTrackArray->AddAt(negtrack1,0);
 430           threeTrackArray->AddAt(postrack1,1);
 431           threeTrackArray->AddAt(negtrack2,2);
 432           io3Prong = Make3Prong(threeTrackArray,esd,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
 433           if(ok3Prong) {
 434             if(fUseTRef) {
 435               AliAODVertex *v = new(verticesHFRef[iVerticesHF++])
 436                 AliAODVertex(*(io3Prong->GetOwnSecondaryVtx()));
 437               v->SetType(AliAODVertex::kUndef);
 438               Double_t px[3]={io3Prong->PxProng(0),io3Prong->PxProng(1),io3Prong->PxProng(2)};
 439               Double_t py[3]={io3Prong->PyProng(0),io3Prong->PyProng(1),io3Prong->PyProng(2)};
 440               Double_t pz[3]={io3Prong->PzProng(0),io3Prong->PzProng(1),io3Prong->PzProng(2)};
 441               Double_t d0[3]={io3Prong->Getd0Prong(0),io3Prong->Getd0Prong(1),io3Prong->Getd0Prong(2)};
 442               Double_t d0err[3]={io3Prong->Getd0errProng(0),io3Prong->Getd0errProng(1),io3Prong->Getd0errProng(2)};
 443               Double_t dcas[3]={io3Prong->GetDCA(0),io3Prong->GetDCA(1),io3Prong->GetDCA(2)};
 444               UShort_t id[3]={(UShort_t)negtrack1->GetID(),(UShort_t)postrack1->GetID(),(UShort_t)negtrack2->GetID()};
 445               AliAODRecoDecayHF3Prong *rd=new(aodCharm3ProngRef[i3Prong++])
 446                 AliAODRecoDecayHF3Prong(v,px,py,pz,d0,d0err,dcas,io3Prong->GetSigmaVert(),io3Prong->GetDist12toPrim(),io3Prong->GetDist23toPrim(),io3Prong->GetCharge());
 447               if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) rd->SetOwnPrimaryVtx(io3Prong->GetOwnPrimaryVtx());
 448               rd->SetProngIDs(3,id);
 449               v->SetParent(rd);
 450             } else {
 451               new(aodD0toKpiRef[i3Prong++]) AliAODRecoDecayHF3Prong(*io3Prong);
 452             }
 453           }
 454           if(io3Prong) { /*delete io3Prong;*/ io3Prong=NULL; }
 455         }
 456         negtrack2 = 0;
 457         delete vertexp1n2;
 458       } // end 2nd loop on negative tracks
 459       twoTrackArray2->Clear();
 460
 461       negtrack1 = 0;
 462       delete vertexp1n1;
 463     } // end 1st loop on negative tracks
 464
 465     postrack1 = 0;
 466   }  // end 1st loop on positive tracks
 467
 468
 469   if(fD0toKpi) {
 470     printf(" D0->Kpi in event %d = %d;\n",
 471            (Int_t)esd->GetEventNumberInFile(),
 472            (Int_t)aodD0toKpiTClArr->GetEntriesFast());
 473   }
 474   if(fJPSItoEle) {
 475     printf(" JPSI->ee in event %d = %d;\n",
 476            (Int_t)esd->GetEventNumberInFile(),
 477            (Int_t)aodJPSItoEleTClArr->GetEntriesFast());
 478   }
 479   if(f3Prong) {
 480     printf(" Charm->3Prong in event %d = %d;\n",
 481            (Int_t)esd->GetEventNumberInFile(),
 482            (Int_t)aodCharm3ProngTClArr->GetEntriesFast());
 483   }
 484   if(f4Prong) {
 485     printf(" Charm->4Prong in event %d = %d;\n",
 486            (Int_t)esd->GetEventNumberInFile(),
 487            (Int_t)aodCharm4ProngTClArr->GetEntriesFast());
 488   }
 489
 490
 491   //printf("delete twoTr 1\n");
 492   twoTrackArray1->Delete(); delete twoTrackArray1;
 493   //printf("delete twoTr 2\n");
 494   twoTrackArray2->Delete(); delete twoTrackArray2;
 495   //printf("delete threeTr 1\n");
 496   threeTrackArray->Clear();
 497   threeTrackArray->Delete(); delete threeTrackArray;
 498   //printf("delete fourTr 1\n");
 499   fourTrackArray->Delete(); delete fourTrackArray;
 500
 501   //------- END SINGLE EVENT ANALYSIS --------------------------------
 502
 503   return;
 504 }
 505 //----------------------------------------------------------------------------
 506 void AliAnalysisVertexingHF::FindCandidates(AliESDEvent *esd,TTree *treeout[])
 507 {
 508   // Find heavy-flavour vertex candidates
 509   //
 510   // DEPRECATED: use FindCandidatesESDtoAOD!
 511
 512   fUseTRef=kFALSE; // cannot use TRefs outside AOD
 513
 514   AliAODRecoDecayHF2Prong *io2Prong = new AliAODRecoDecayHF2Prong();
 515   AliAODRecoDecayHF3Prong *io3Prong = new AliAODRecoDecayHF3Prong();
 516   AliAODRecoDecayHF4Prong *io4Prong = new AliAODRecoDecayHF4Prong();
 517   Int_t itree=0;
 518   Int_t itreeD0toKpi=-1,itreeJPSItoEle=-1,itree3Prong=-1,itree4Prong=-1;
 519   Int_t initEntriesD0toKpi=0,initEntriesJPSItoEle=0,initEntries3Prong=0,initEntries4Prong=0;
 520   if(fD0toKpi) {
 521     itreeD0toKpi=itree;
 522     treeout[itree]->SetBranchAddress("D0toKpi",&io2Prong);
 523     itree++;
 524     initEntriesD0toKpi = treeout[itreeD0toKpi]->GetEntries();
 525   }
 526   if(fJPSItoEle) {
 527     itreeJPSItoEle=itree;
 528     treeout[itree]->SetBranchAddress("JPSItoEle",&io2Prong);
 529     itree++;
 530     initEntriesJPSItoEle = treeout[itreeJPSItoEle]->GetEntries();
 531   }
 532   if(f3Prong) {
 533     itree3Prong=itree;
 534     treeout[itree]->SetBranchAddress("Charmto3Prong",&io3Prong);
 535     itree++;
 536     initEntries3Prong = treeout[itree3Prong]->GetEntries();
 537   }
 538   if(f4Prong) {
 539     itree4Prong=itree;
 540     treeout[itree]->SetBranchAddress("D0to4Prong",&io4Prong);
 541     itree++;
 542     initEntries4Prong = treeout[itree4Prong]->GetEntries();
 543   }
 544   delete io2Prong; io2Prong = NULL;
 545   delete io3Prong; io3Prong = NULL;
 546   delete io4Prong; io4Prong = NULL;
 547
 548   Int_t    iTrkP1,iTrkP2,iTrkN1,iTrkN2,trkEntries;
 549   Int_t    nTrksP=0,nTrksN=0;
 550   Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2;
 551   Bool_t   okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
 552   AliESDtrack *postrack1 = 0;
 553   AliESDtrack *postrack2 = 0;
 554   AliESDtrack *negtrack1 = 0;
 555   AliESDtrack *negtrack2 = 0;
 556   Double_t dcaMax = fD0toKpiCuts[1];
 557   if(dcaMax<fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1];
 558   if(dcaMax<fDplusCuts[11]) dcaMax=fDplusCuts[11];
 559   if(fDebug) printf(" dca cut set to %f cm\n",dcaMax);
 560
 561   Int_t ev = (Int_t)esd->GetEventNumberInFile();
 562   printf("--- Finding candidates in event %d\n",ev);
 563
 564   fBzkG = (Double_t)esd->GetMagneticField();
 565
 566   trkEntries = (Int_t)esd->GetNumberOfTracks();
 567   printf(" Number of tracks: %d\n",trkEntries);
 568   if(trkEntries<2) return;
 569
 570   // retrieve primary vertex from the AliESDEvent
 571   if(!esd->GetPrimaryVertex()) {
 572     printf(" No vertex in AliESD\n");
 573     return;
 574   }
 575   AliESDVertex copy(*(esd->GetPrimaryVertex()));
 576   SetPrimaryVertex(&copy);
 577
 578   // call function which applies sigle-track selection and
 579   // separetes positives and negatives
 580   TObjArray trksP(trkEntries/2);
 581   TObjArray trksN(trkEntries/2);
 582   SelectTracks(esd,trksP,nTrksP,
 583                    trksN,nTrksN);
 584
 585   printf(" Pos. tracks: %d    Neg. tracks: %d\n",nTrksP,nTrksN);
 586
 587   TObjArray *twoTrackArray1 = new TObjArray(2);
 588   TObjArray *twoTrackArray2 = new TObjArray(2);
 589   TObjArray *threeTrackArray = new TObjArray(3);
 590   TObjArray *fourTrackArray = new TObjArray(4);
 591
 592   // LOOP ON  POSITIVE  TRACKS
 593   for(iTrkP1=0; iTrkP1<nTrksP; iTrkP1++) {
 594     if(fDebug) if(iTrkP1%1==0) printf("  Processing positive track number %d of %d\n",iTrkP1,nTrksP);
 595     // get track from track array
 596     postrack1 = (AliESDtrack*)trksP.UncheckedAt(iTrkP1);
 597
 598     // LOOP ON  NEGATIVE  TRACKS
 599     for(iTrkN1=0; iTrkN1<nTrksN; iTrkN1++) {
 600       if(fDebug) if(iTrkN1%1==0) printf("    Processing negative track number %d of %d\n",iTrkN1,nTrksN);
 601       // get track from tracks array
 602       negtrack1 = (AliESDtrack*)trksN.UncheckedAt(iTrkN1);
 603       // DCA between the two tracks
 604       dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 605       if(dcap1n1>dcaMax) { negtrack1=0; continue; }
 606       // Vertexing
 607       twoTrackArray1->AddAt(postrack1,0);
 608       twoTrackArray1->AddAt(negtrack1,1);
 609       AliESDVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1);
 610       if(!vertexp1n1) {
 611         twoTrackArray1->Clear();
 612         negtrack1=0;
 613         continue;
 614       }
 615       if(fD0toKpi || fJPSItoEle) {
 616         io2Prong = Make2Prong(twoTrackArray1,esd,vertexp1n1,dcap1n1,okD0,okJPSI);
 617         if(okD0)   treeout[itreeD0toKpi]->Fill();
 618         if(okJPSI) treeout[itreeJPSItoEle]->Fill();
 619         delete io2Prong; io2Prong=NULL;
 620       }
 621
 622       twoTrackArray1->Clear();
 623       if(!f3Prong && !f4Prong)  {
 624         negtrack1=0;
 625         delete vertexp1n1;
 626         continue;
 627       }
 628
 629       // 2nd LOOP  ON  POSITIVE  TRACKS
 630       for(iTrkP2=iTrkP1+1; iTrkP2<nTrksP; iTrkP2++) {
 631         // get track from tracks array
 632         postrack2 = (AliESDtrack*)trksP.UncheckedAt(iTrkP2);
 633         dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 634         if(dcap2n1>dcaMax) { postrack2=0; continue; }
 635         dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
 636         if(dcap1p2>dcaMax) { postrack2=0; continue; }
 637
 638         // Vertexing
 639         twoTrackArray2->AddAt(postrack2,0);
 640         twoTrackArray2->AddAt(negtrack1,1);
 641         AliESDVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2);
 642         if(!vertexp2n1) {
 643           twoTrackArray2->Clear();
 644           postrack2=0;
 645           continue;
 646         }
 647         if(f3Prong) {
 648           threeTrackArray->AddAt(postrack1,0);
 649           threeTrackArray->AddAt(negtrack1,1);
 650           threeTrackArray->AddAt(postrack2,2);
 651           io3Prong = Make3Prong(threeTrackArray,esd,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
 652           if(ok3Prong) treeout[itree3Prong]->Fill();
 653           if(io3Prong) delete io3Prong; io3Prong=NULL;
 654         }
 655         if(f4Prong) {
 656           // 3rd LOOP  ON  NEGATIVE  TRACKS (for 4 prong)
 657           for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) {
 658             // get track from tracks array
 659             negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2);
 660             dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 661             if(dcap1n2>dcaMax) { negtrack2=0; continue; }
 662             // Vertexing
 663             fourTrackArray->AddAt(postrack1,0);
 664             fourTrackArray->AddAt(negtrack1,1);
 665             fourTrackArray->AddAt(postrack2,2);
 666             fourTrackArray->AddAt(negtrack2,3);
 667             io4Prong = Make4Prong(fourTrackArray,esd,vertexp1n1,vertexp2n1,dcap1n1,dcap1n2,dcap2n1,ok4Prong);
 668             if(ok4Prong) treeout[itree4Prong]->Fill();
 669             delete io4Prong; io4Prong=NULL;
 670             fourTrackArray->Clear();
 671             negtrack2 = 0;
 672           } // end loop on negative tracks
 673         }
 674         postrack2 = 0;
 675         delete vertexp2n1;
 676       } // end 2nd loop on positive tracks
 677       twoTrackArray2->Clear();
 678
 679       // 2nd LOOP  ON  NEGATIVE  TRACKS
 680       for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) {
 681         // get track from tracks array
 682         negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2);
 683         dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 684         if(dcap1n2>dcaMax) { negtrack2=0; continue; }
 685         dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 686         if(dcan1n2>dcaMax) { negtrack2=0; continue; }
 687
 688         // Vertexing
 689         twoTrackArray2->AddAt(postrack1,0);
 690         twoTrackArray2->AddAt(negtrack2,1);
 691         AliESDVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2);
 692         if(!vertexp1n2) {
 693           twoTrackArray2->Clear();
 694           negtrack2=0;
 695           continue;
 696         }
 697         if(f3Prong) {
 698           threeTrackArray->AddAt(negtrack1,0);
 699           threeTrackArray->AddAt(postrack1,1);
 700           threeTrackArray->AddAt(negtrack2,2);
 701           io3Prong = Make3Prong(threeTrackArray,esd,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
 702           if(ok3Prong) treeout[itree3Prong]->Fill();
 703           if(io3Prong) delete io3Prong; io3Prong=NULL;
 704           }
 705         negtrack2 = 0;
 706         delete vertexp1n2;
 707       } // end 2nd loop on negative tracks
 708       twoTrackArray2->Clear();
 709
 710       negtrack1 = 0;
 711       delete vertexp1n1;
 712     } // end 1st loop on negative tracks
 713
 714     postrack1 = 0;
 715   }  // end 1st loop on positive tracks
 716
 717
 718
 719   //printf("delete twoTr 1\n");
 720   twoTrackArray1->Delete(); delete twoTrackArray1;
 721   //printf("delete twoTr 2\n");
 722   twoTrackArray2->Delete(); delete twoTrackArray2;
 723   //printf("delete threeTr 1\n");
 724   threeTrackArray->Clear();
 725   threeTrackArray->Delete(); delete threeTrackArray;
 726   //printf("delete fourTr 1\n");
 727   fourTrackArray->Delete(); delete fourTrackArray;
 728
 729
 730   // create a copy of this class to be written to output file
 731   //AliAnalysisVertexingHF *copy = (AliAnalysisVertexingHF*)this->Clone("AnalysisVertexingHF");
 732
 733   // print statistics
 734   if(fD0toKpi) {
 735     printf(" D0->Kpi: event %d = %d; total = %d;\n",
 736            (Int_t)esd->GetEventNumberInFile(),
 737            (Int_t)treeout[itreeD0toKpi]->GetEntries()-initEntriesD0toKpi,
 738            (Int_t)treeout[itreeD0toKpi]->GetEntries());
 739   }
 740   if(fJPSItoEle) {
 741     printf(" JPSI->ee: event %d = %d; total = %d;\n",
 742            (Int_t)esd->GetEventNumberInFile(),
 743            (Int_t)treeout[itreeJPSItoEle]->GetEntries()-initEntriesJPSItoEle,
 744            (Int_t)treeout[itreeJPSItoEle]->GetEntries());
 745   }
 746   if(f3Prong) {
 747     printf(" Charm->3Prong: event %d = %d; total = %d;\n",
 748            (Int_t)esd->GetEventNumberInFile(),
 749            (Int_t)treeout[itree3Prong]->GetEntries()-initEntries3Prong,
 750            (Int_t)treeout[itree3Prong]->GetEntries());
 751   }
 752   if(f4Prong) {
 753     printf(" Charm->4Prong: event %d = %d; total = %d;\n",
 754            (Int_t)esd->GetEventNumberInFile(),
 755            (Int_t)treeout[itree4Prong]->GetEntries()-initEntries4Prong,
 756            (Int_t)treeout[itree4Prong]->GetEntries());
 757   }
 758
 759
 760   return;
 761 }
 762 //----------------------------------------------------------------------------
 763 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
 764                                    TObjArray *twoTrackArray1,AliESDEvent *esd,
 765                                    AliESDVertex *secVertexESD,Double_t dca,
 766                                    Bool_t &okD0,Bool_t &okJPSI) const
 767 {
 768   // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
 769   // reconstruction cuts
 770   // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
 771
 772   okD0=kFALSE; okJPSI=kFALSE;
 773
 774   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
 775
 776   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(0);
 777   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(1);
 778
 779   // propagate tracks to secondary vertex, to compute inv. mass
 780   postrack->RelateToVertex(secVertexESD,fBzkG,10.);
 781   negtrack->RelateToVertex(secVertexESD,fBzkG,10.);
 782
 783   Double_t momentum[3];
 784   postrack->GetPxPyPz(momentum);
 785   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 786   negtrack->GetPxPyPz(momentum);
 787   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 788
 789
 790   // invariant mass cut (try to improve coding here..)
 791   Bool_t okMassCut=kFALSE;
 792   if(!okMassCut && fD0toKpi) if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
 793   if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
 794   if(!okMassCut) {
 795     if(fDebug) printf(" candidate didn't pass mass cut\n");
 796     return 0x0;
 797   }
 798
 799
 800   AliESDVertex *primVertex = fV1;
 801   AliESDVertex *ownPrimVertex=0;
 802
 803   // primary vertex from *other* tracks in the event
 804   if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
 805     ownPrimVertex = OwnPrimaryVertex(2,twoTrackArray1,esd);
 806     if(!ownPrimVertex) {
 807       return 0x0;
 808     } else {
 809       if(ownPrimVertex->GetNContributors()<2) {
 810         delete ownPrimVertex;
 811         return 0x0;
 812       } else {
 813         primVertex = ownPrimVertex;
 814       }
 815     }
 816   }
 817
 818   Float_t d0z0[2],covd0z0[3];
 819   postrack->RelateToVertex(primVertex,fBzkG,10.);
 820   postrack->GetImpactParameters(d0z0,covd0z0);
 821   d0[0] = d0z0[0];
 822   d0err[0] = TMath::Sqrt(covd0z0[0]);
 823   negtrack->RelateToVertex(primVertex,fBzkG,10.);
 824   negtrack->GetImpactParameters(d0z0,covd0z0);
 825   d0[1] = d0z0[0];
 826   d0err[1] = TMath::Sqrt(covd0z0[0]);
 827
 828   // create the object AliAODRecoDecayHF2Prong
 829   Double_t pos[3],cov[6];
 830   secVertexESD->GetXYZ(pos); // position
 831   secVertexESD->GetCovMatrix(cov); //covariance matrix
 832   AliAODVertex *secVertexAOD = new AliAODVertex(pos,cov,secVertexESD->GetChi2toNDF());
 833   AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVertexAOD,px,py,pz,d0,d0err,dca);
 834   the2Prong->SetOwnSecondaryVtx(secVertexAOD);
 835   primVertex->GetXYZ(pos); // position
 836   primVertex->GetCovMatrix(cov); //covariance matrix
 837   AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF());
 838   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
 839
 840
 841   // select D0->Kpi
 842   Int_t checkD0,checkD0bar;
 843   if(fD0toKpi) okD0 = the2Prong->SelectD0(fD0toKpiCuts,checkD0,checkD0bar);
 844   //if(fDebug && fD0toKpi) printf("   %d\n",(Int_t)okD0);
 845   // select J/psi from B
 846   Int_t checkJPSI;
 847   if(fJPSItoEle) okJPSI = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI);
 848   //if(fDebug && fJPSItoEle) printf("   %d\n",(Int_t)okJPSI);
 849
 850
 851   if(okD0 || okJPSI) {
 852     // get PID info from ESD
 853     Double_t esdpid0[5];
 854     postrack->GetESDpid(esdpid0);
 855     Double_t esdpid1[5];
 856     negtrack->GetESDpid(esdpid1);
 857     Double_t esdpid[10];
 858     for(Int_t i=0;i<5;i++) {
 859       esdpid[i]   = esdpid0[i];
 860       esdpid[5+i] = esdpid1[i];
 861     }
 862     the2Prong->SetPID(2,esdpid);
 863   }
 864
 865   if(ownPrimVertex) delete ownPrimVertex;
 866
 867   return the2Prong;
 868 }
 869 //----------------------------------------------------------------------------
 870 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
 871                              TObjArray *threeTrackArray,AliESDEvent *esd,
 872                              AliESDVertex *vertexp1n1,AliESDVertex *vertexp2n1,
 873                              Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
 874                              Bool_t &ok3Prong) const
 875 {
 876   // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
 877   // reconstruction cuts
 878   // E.Bruna, F.Prino
 879
 880   ok3Prong=kFALSE;
 881   Double_t px[3],py[3],pz[3],d0[3],d0err[3];//d0z[3];
 882   Float_t d0z0[2],covd0z0[3];
 883
 884
 885   AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
 886   AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
 887   AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
 888
 889   AliESDVertex *primVertex = fV1;
 890
 891   postrack1->RelateToVertex(primVertex,fBzkG,10.);
 892   negtrack->RelateToVertex(primVertex,fBzkG,10.);
 893   postrack2->RelateToVertex(primVertex,fBzkG,10.);
 894
 895   Double_t momentum[3];
 896   postrack1->GetPxPyPz(momentum);
 897   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 898   negtrack->GetPxPyPz(momentum);
 899   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 900   postrack2->GetPxPyPz(momentum);
 901   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
 902
 903   postrack1->GetImpactParameters(d0z0,covd0z0);
 904   d0[0]=d0z0[0];
 905   d0err[0] = TMath::Sqrt(covd0z0[0]);
 906   negtrack->GetImpactParameters(d0z0,covd0z0);
 907   d0[1]=d0z0[0];
 908   d0err[1] = TMath::Sqrt(covd0z0[0]);
 909   postrack2->GetImpactParameters(d0z0,covd0z0);
 910   d0[2]=d0z0[0];
 911   d0err[2] = TMath::Sqrt(covd0z0[0]);
 912
 913
 914   // invariant mass cut for D+, Ds, Lc
 915   Bool_t okMassCut=kFALSE;
 916   if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
 917   if(!okMassCut) {
 918     if(fDebug) printf(" candidate didn't pass mass cut\n");
 919     return 0x0;
 920   }
 921
 922   //charge
 923   Short_t charge=(Short_t)(postrack1->GetSign()*postrack2->GetSign()*negtrack->GetSign());
 924
 925   AliESDVertex *ownPrimVertex = 0;
 926   // primary vertex from *other* tracks in the event
 927   if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
 928     ownPrimVertex = OwnPrimaryVertex(3,threeTrackArray,esd);
 929     if(!ownPrimVertex) {
 930       return 0x0;
 931     } else {
 932       if(ownPrimVertex->GetNContributors()<2) {
 933         delete ownPrimVertex;
 934         return 0x0;
 935       } else {
 936         primVertex = ownPrimVertex;
 937       }
 938     }
 939   }
 940
 941   // create the object AliAODRecoDecayHF3Prong
 942   AliESDVertex* secVert3Prong = ReconstructSecondaryVertex(threeTrackArray);
 943   if(!secVert3Prong) {
 944     if(ownPrimVertex) delete ownPrimVertex;
 945     return 0x0;
 946   }
 947   Double_t pos[3],cov[6],sigmavert;
 948   secVert3Prong->GetXYZ(pos); // position
 949   secVert3Prong->GetCovMatrix(cov); //covariance matrix
 950   sigmavert=secVert3Prong->GetDispersion();
 951
 952   AliAODVertex *secVert3PrAOD = new AliAODVertex(pos,cov,secVert3Prong->GetChi2toNDF());
 953   primVertex->GetXYZ(pos); // position
 954   primVertex->GetCovMatrix(cov); //covariance matrix
 955   AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF());
 956   Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
 957
 958   Double_t dist12=TMath::Sqrt((vertexp1n1->GetXv()-pos[0])*(vertexp1n1->GetXv()-pos[0])+(vertexp1n1->GetYv()-pos[1])*(vertexp1n1->GetYv()-pos[1])+(vertexp1n1->GetZv()-pos[2])*(vertexp1n1->GetZv()-pos[2]));
 959   Double_t dist23=TMath::Sqrt((vertexp2n1->GetXv()-pos[0])*(vertexp2n1->GetXv()-pos[0])+(vertexp2n1->GetYv()-pos[1])*(vertexp2n1->GetYv()-pos[1])+(vertexp2n1->GetZv()-pos[2])*(vertexp2n1->GetZv()-pos[2]));
 960
 961   AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert3PrAOD,px,py,pz,d0,d0err,dca,sigmavert,dist12,dist23,charge);
 962   the3Prong->SetOwnSecondaryVtx(secVert3PrAOD);
 963   the3Prong->SetOwnPrimaryVtx(primVertexAOD);
 964
 965
 966   // select D+->Kpipi, Ds->KKpi, Lc->pKpi
 967   if(f3Prong) {
 968     ok3Prong = kFALSE;
 969     Int_t ok1,ok2;
 970     if(the3Prong->SelectDplus(fDplusCuts))   ok3Prong = kTRUE;
 971     if(the3Prong->SelectDs(fDsCuts,ok1,ok2)) ok3Prong = kTRUE;
 972     if(the3Prong->SelectLc(fLcCuts,ok1,ok2)) ok3Prong = kTRUE;
 973   }
 974   //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
 975   if(ok3Prong) {
 976     // get PID info from ESD
 977     Double_t esdpid0[5];
 978     postrack1->GetESDpid(esdpid0);
 979     Double_t esdpid1[5];
 980     negtrack->GetESDpid(esdpid1);
 981     Double_t esdpid2[5];
 982     postrack2->GetESDpid(esdpid2);
 983
 984
 985     Double_t esdpid[15];
 986     for(Int_t i=0;i<5;i++) {
 987       esdpid[i]   = esdpid0[i];
 988       esdpid[5+i] = esdpid1[i];
 989       esdpid[10+i] = esdpid2[i];
 990     }
 991     the3Prong->SetPID(3,esdpid);
 992   }
 993
 994   if(ownPrimVertex) delete ownPrimVertex;
 995
 996   return the3Prong;
 997 }
 998 //----------------------------------------------------------------------------
 999 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1000                              TObjArray *fourTrackArray,AliESDEvent *esd,
1001                              AliESDVertex *vertexp1n1,AliESDVertex *vertexp2n1,
1002                              Double_t dcap1n1,Double_t dcap1n2,Double_t dcap2n1,
1003                              Bool_t &ok4Prong) const
1004 {
1005   // Make 4Prong candidates and check if they pass D0toKpipipi
1006   // reconstruction cuts
1007   // G.E.Bruno, R.Romita
1008
1009   ok4Prong=kFALSE;
1010
1011   Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1012   //Float_t d0z0[2],covd0z0[3];
1013
1014   px[0]=dcap1n1*dcap1n2*dcap2n1; // TO BE CHANGED (done just to removed compilation warning about dca... not used)
1015
1016   //charge
1017   Short_t charge=0;
1018
1019   AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1020   AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1021   AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1022   AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1023
1024   AliESDVertex *primVertex = fV1;
1025
1026   postrack1->RelateToVertex(primVertex,fBzkG,10.);
1027   negtrack1->RelateToVertex(primVertex,fBzkG,10.);
1028   postrack2->RelateToVertex(primVertex,fBzkG,10.);
1029   negtrack2->RelateToVertex(primVertex,fBzkG,10.);
1030
1031   Double_t momentum[3];
1032   postrack1->GetPxPyPz(momentum);
1033   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1034   negtrack1->GetPxPyPz(momentum);
1035   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1036   postrack2->GetPxPyPz(momentum);
1037   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1038   negtrack2->GetPxPyPz(momentum);
1039   px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1040
1041   // invariant mass cut for rho or D0 (try to improve coding here..)
1042   //Bool_t okMassCut=kFALSE;
1043   //if(!okMassCut) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
1044   //if(!okMassCut) {
1045   //  if(fDebug) printf(" candidate didn't pass mass cut\n");
1046   //  return 0x0;
1047   //}
1048
1049   AliESDVertex *ownPrimVertex = 0;
1050   // primary vertex from *other* tracks in the event
1051   if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
1052     ownPrimVertex = OwnPrimaryVertex(4,fourTrackArray,esd);
1053     if(!ownPrimVertex) {
1054       return 0x0;
1055     } else {
1056       if(ownPrimVertex->GetNContributors()<2) {
1057         delete ownPrimVertex;
1058         return 0x0;
1059       } else {
1060         primVertex = ownPrimVertex;
1061       }
1062     }
1063   }
1064
1065   // create the object AliAODRecoDecayHF4Prong
1066   AliESDVertex* secVert4Prong = ReconstructSecondaryVertex(fourTrackArray);
1067   if(!secVert4Prong) {
1068     if(ownPrimVertex) delete ownPrimVertex;
1069     return 0x0;
1070   }
1071   Double_t pos[3],cov[6],sigmavert;
1072   secVert4Prong->GetXYZ(pos); // position
1073   secVert4Prong->GetCovMatrix(cov); //covariance matrix
1074   sigmavert=secVert4Prong->GetDispersion();
1075
1076   AliAODVertex *secVert4PrAOD = new AliAODVertex(pos,cov,secVert4Prong->GetChi2toNDF());
1077   primVertex->GetXYZ(pos); // position
1078   primVertex->GetCovMatrix(cov); //covariance matrix
1079   AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF());
1080   //Double_t dca[6]={dcap1n1,dcap2n1,dcap1p2,0.,0.,0.}; //
1081   Double_t dca[6]={0.,0.,0.,0.,0.,0.}; //  modify it
1082
1083   Double_t dist12=TMath::Sqrt((vertexp1n1->GetXv()-pos[0])*(vertexp1n1->GetXv()-pos[0])+(vertexp1n1->GetYv()-pos[1])*(vertexp1n1->GetYv()-pos[1])+(vertexp1n1->GetZv()-pos[2])*(vertexp1n1->GetZv()-pos[2]));
1084   Double_t dist23=TMath::Sqrt((vertexp2n1->GetXv()-pos[0])*(vertexp2n1->GetXv()-pos[0])+(vertexp2n1->GetYv()-pos[1])*(vertexp2n1->GetYv()-pos[1])+(vertexp2n1->GetZv()-pos[2])*(vertexp2n1->GetZv()-pos[2]));
1085   Double_t dist14=0.; // to be implemented
1086   Double_t dist34=0.; // to be implemented
1087
1088   //AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert4PrAOD,px,py,pz,d0,d0err,dca,sigmavert,dist12,dist23,charge);
1089   AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert4PrAOD,px,py,pz,d0,d0err,dca,dist12,dist23,dist14,dist34,charge);
1090   the4Prong->SetOwnPrimaryVtx(primVertexAOD);
1091   the4Prong->SetOwnSecondaryVtx(secVert4PrAOD);
1092
1093
1094   // use the following two lines once AliAODRecoDecayHF4Prong::SelectD0 is available
1095   // select D0->Kpipipi
1096   //Int_t checkD0,checkD0bar;
1097   // ok4Prong=the4Prong->SelectD0(fD04pCuts,checkD0,checkD0bar);
1098   ok4Prong=kFALSE;  //for the time being ...
1099
1100
1101   // get PID info from ESD
1102   Double_t esdpid0[5];
1103   postrack1->GetESDpid(esdpid0);
1104   Double_t esdpid1[5];
1105   negtrack1->GetESDpid(esdpid1);
1106   Double_t esdpid2[5];
1107   postrack2->GetESDpid(esdpid2);
1108   Double_t esdpid3[5];
1109   negtrack2->GetESDpid(esdpid3);
1110
1111   Double_t esdpid[20];
1112   for(Int_t i=0;i<5;i++) {
1113     esdpid[i]   = esdpid0[i];
1114     esdpid[5+i] = esdpid1[i];
1115     esdpid[10+i] = esdpid2[i];
1116     esdpid[15+i] = esdpid3[i];
1117   }
1118   the4Prong->SetPID(4,esdpid);
1119
1120   if(ownPrimVertex) delete ownPrimVertex;
1121
1122   return the4Prong;
1123 }
1124 //-----------------------------------------------------------------------------
1125 AliESDVertex* AliAnalysisVertexingHF::OwnPrimaryVertex(Int_t ntrks,
1126                                                        TObjArray *trkArray,
1127                                                        AliESDEvent *esd) const
1128 {
1129   // Returns primary vertex specific to this candidate
1130
1131   AliVertexerTracks *vertexer1 = new AliVertexerTracks(esd->GetMagneticField());
1132   AliESDVertex *ownPrimVertex = 0;
1133
1134   // recalculating the vertex
1135   if(fRecoPrimVtxSkippingTrks) {
1136     if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1137       Float_t diamondcovxy[3];
1138       esd->GetDiamondCovXY(diamondcovxy);
1139       Double_t pos[3]={esd->GetDiamondX(),esd->GetDiamondY(),0.};
1140       Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.};
1141       AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1142       vertexer1->SetVtxStart(diamond);
1143       delete diamond; diamond=NULL;
1144       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1145         vertexer1->SetOnlyFitter();
1146     }
1147     Int_t skipped[10];
1148     AliESDtrack *t = 0;
1149     for(Int_t i=0; i<ntrks; i++) {
1150       t = (AliESDtrack*)trkArray->UncheckedAt(i);
1151       skipped[i] = (Int_t)t->GetID();
1152     }
1153     vertexer1->SetSkipTracks(ntrks,skipped);
1154     ownPrimVertex = (AliESDVertex*)vertexer1->FindPrimaryVertex(esd);
1155   }
1156
1157   // removing the prongs tracks
1158   if(fRmTrksFromPrimVtx) {
1159     TObjArray rmArray(ntrks);
1160     UShort_t *rmId = new UShort_t[ntrks];
1161     AliESDtrack *esdTrack = 0;
1162     AliESDtrack *t = 0;
1163     for(Int_t i=0; i<ntrks; i++) {
1164       t = (AliESDtrack*)trkArray->UncheckedAt(i);
1165       esdTrack = new AliESDtrack(*t);
1166       rmArray.AddLast(esdTrack);
1167       rmId[i]=(UShort_t)esdTrack->GetID();
1168       delete esdTrack;
1169     }
1170     Float_t diamondxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
1171     ownPrimVertex = vertexer1->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1172     delete [] rmId; rmId=NULL;
1173     rmArray.Delete();
1174   }
1175
1176   delete vertexer1; vertexer1=NULL;
1177
1178   return ownPrimVertex;
1179 }
1180 //-----------------------------------------------------------------------------
1181 void AliAnalysisVertexingHF::PrintStatus() const {
1182   // Print parameters being used
1183
1184   printf("Preselections:\n");
1185   printf("    fITSrefit   = %d\n",(Int_t)fITSrefit);
1186   printf("    fBothSPD   = %d\n",(Int_t)fBothSPD);
1187   printf("    fMinITSCls   = %d\n",fMinITSCls);
1188   printf("    fMinPtCut   = %f GeV/c\n",fMinPtCut);
1189   printf("    fMind0rphiCut   = %f cm\n",fMind0rphiCut);
1190   if(fSecVtxWithKF) {
1191     printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1192   } else {
1193     printf("Secondary vertex with AliVertexerTracks\n");
1194   }
1195   if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1196   if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1197   if(fD0toKpi) {
1198     printf("Reconstruct D0->Kpi candidates with cuts:\n");
1199     printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
1200     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
1201     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
1202     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
1203     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
1204     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
1205     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
1206     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
1207     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
1208   }
1209   if(fJPSItoEle) {
1210     printf("Reconstruct J/psi from B candidates with cuts:\n");
1211     printf("    |M-MJPSI| [GeV]    < %f\n",fBtoJPSICuts[0]);
1212     printf("    dca    [cm]  < %f\n",fBtoJPSICuts[1]);
1213     printf("    cosThetaStar     < %f\n",fBtoJPSICuts[2]);
1214     printf("    pTP     [GeV/c]    > %f\n",fBtoJPSICuts[3]);
1215     printf("    pTN    [GeV/c]    > %f\n",fBtoJPSICuts[4]);
1216     printf("    |d0P|  [cm]  < %f\n",fBtoJPSICuts[5]);
1217     printf("    |d0N| [cm]  < %f\n",fBtoJPSICuts[6]);
1218     printf("    d0d0  [cm^2] < %f\n",fBtoJPSICuts[7]);
1219     printf("    cosThetaPoint    > %f\n",fBtoJPSICuts[8]);
1220   }
1221   if(f3Prong) {
1222     printf("Reconstruct 3 prong candidates.\n");
1223     printf("  D+->Kpipi cuts:\n");
1224     printf("    |M-MD+| [GeV]    < %f\n",fDplusCuts[0]);
1225     printf("    pTK     [GeV/c]    > %f\n",fDplusCuts[1]);
1226     printf("    pTPi    [GeV/c]    > %f\n",fDplusCuts[2]);
1227     printf("    |d0K|  [cm]  > %f\n",fDplusCuts[3]);
1228     printf("    |d0Pi| [cm]  > %f\n",fDplusCuts[4]);
1229     printf("    dist12    [cm]  < %f\n",fDplusCuts[5]);
1230     printf("    sigmavert [cm]   < %f\n",fDplusCuts[6]);
1231     printf("    dist prim-sec [cm] > %f\n",fDplusCuts[7]);
1232     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDplusCuts[8]);
1233     printf("    cosThetaPoint    > %f\n",fDplusCuts[9]);
1234     printf("    Sum d0^2 [cm^2]  > %f\n",fDplusCuts[10]);
1235     printf("    dca cut [cm]  < %f\n",fDplusCuts[11]);
1236     printf("  Ds->KKpi cuts:\n");
1237     printf("    |M-MDs| [GeV]    < %f\n",fDsCuts[0]);
1238     printf("  Lc->pKpi cuts:\n");
1239     printf("    |M-MD+| [GeV]    < %f\n",fLcCuts[0]);
1240   }
1241
1242   return;
1243 }
1244 //-----------------------------------------------------------------------------
1245 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1246                                              Int_t nprongs,
1247                                              Double_t *px,
1248                                              Double_t *py,
1249                                              Double_t *pz) const {
1250   // Check invariant mass cut
1251
1252   Short_t dummycharge=0;
1253   Double_t *dummyd0 = new Double_t[nprongs];
1254   for(Int_t ip=0;ip<nprongs;ip++) dummyd0[ip]=0.;
1255   AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
1256   delete [] dummyd0;
1257
1258   UInt_t pdg2[2],pdg3[3];
1259   Double_t mPDG,minv;
1260
1261   Bool_t retval=kFALSE;
1262   switch (decay)
1263     {
1264     case 0:                  // D0->Kpi
1265       pdg2[0]=211; pdg2[1]=321;
1266       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1267       minv = rd->InvMass(nprongs,pdg2);
1268       if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1269       pdg2[0]=321; pdg2[1]=211;
1270       minv = rd->InvMass(nprongs,pdg2);
1271       if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1272       break;
1273     case 1:                  // JPSI->ee
1274       pdg2[0]=11; pdg2[1]=11;
1275       mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
1276       minv = rd->InvMass(nprongs,pdg2);
1277       if(TMath::Abs(minv-mPDG)<fBtoJPSICuts[0]) retval=kTRUE;
1278       break;
1279     case 2:                  // D+->Kpipi
1280       pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
1281       mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
1282       minv = rd->InvMass(nprongs,pdg3);
1283       if(TMath::Abs(minv-mPDG)<fDplusCuts[0]) retval=kTRUE;
1284                             // Ds+->KKpi
1285       pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
1286       mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1287       minv = rd->InvMass(nprongs,pdg3);
1288       if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1289       pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
1290       minv = rd->InvMass(nprongs,pdg3);
1291       if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1292                             // Lc->pKpi
1293       pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
1294       mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1295       minv = rd->InvMass(nprongs,pdg3);
1296       if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1297       pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
1298       minv = rd->InvMass(nprongs,pdg3);
1299       if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1300       break;
1301     default:
1302       break;
1303     }
1304
1305   delete rd;
1306
1307   return retval;
1308 }
1309 //-----------------------------------------------------------------------------
1310 void AliAnalysisVertexingHF::SelectTracks(AliESDEvent *esd,
1311                                           TObjArray &trksP,Int_t &nTrksP,
1312                                           TObjArray &trksN,Int_t &nTrksN) const
1313 {
1314   // Fill two TObjArrays with positive and negative tracks and
1315   // apply single-track preselection
1316
1317   nTrksP=0,nTrksN=0;
1318
1319   Int_t entries = (Int_t)esd->GetNumberOfTracks();
1320
1321   // transfer ITS tracks from ESD to arrays and to a tree
1322   for(Int_t i=0; i<entries; i++) {
1323
1324     AliESDtrack *esdtrack = esd->GetTrack(i);
1325     UInt_t status = esdtrack->GetStatus();
1326
1327     // require refit in ITS
1328     if(fITSrefit && !(status&AliESDtrack::kITSrefit)) {
1329       if(fDebug) printf("track %d is not kITSrefit\n",i);
1330       continue;
1331     }
1332
1333     // require minimum # of ITS points
1334     if(esdtrack->GetNcls(0)<fMinITSCls)  {
1335       if(fDebug) printf("track %d has %d ITS cls\n",i,esdtrack->GetNcls(0));
1336       continue;
1337     }
1338     // require points on the 2 pixel layers
1339     if(fBothSPD)
1340       if(!TESTBIT(esdtrack->GetITSClusterMap(),0) ||
1341          !TESTBIT(esdtrack->GetITSClusterMap(),1)) continue;
1342
1343     // single track selection
1344     if(!SingleTrkCuts(*esdtrack)) continue;
1345
1346     if(esdtrack->GetSign()<0) { // negative track
1347       trksN.AddLast(esdtrack);
1348       nTrksN++;
1349     } else {                 // positive track
1350       trksP.AddLast(esdtrack);
1351       nTrksP++;
1352     }
1353
1354   } // loop on ESD tracks
1355
1356   return;
1357 }
1358 //-----------------------------------------------------------------------------
1359 void AliAnalysisVertexingHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
1360                                    Double_t cut2,Double_t cut3,Double_t cut4,
1361                                    Double_t cut5,Double_t cut6,
1362                                    Double_t cut7,Double_t cut8)
1363 {
1364   // Set the cuts for D0 selection
1365   fD0toKpiCuts[0] = cut0;
1366   fD0toKpiCuts[1] = cut1;
1367   fD0toKpiCuts[2] = cut2;
1368   fD0toKpiCuts[3] = cut3;
1369   fD0toKpiCuts[4] = cut4;
1370   fD0toKpiCuts[5] = cut5;
1371   fD0toKpiCuts[6] = cut6;
1372   fD0toKpiCuts[7] = cut7;
1373   fD0toKpiCuts[8] = cut8;
1374
1375   return;
1376 }
1377 //-----------------------------------------------------------------------------
1378 void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9])
1379 {
1380   // Set the cuts for D0 selection
1381
1382   for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
1383
1384   return;
1385 }
1386 //-----------------------------------------------------------------------------
1387 void AliAnalysisVertexingHF::SetBtoJPSICuts(Double_t cut0,Double_t cut1,
1388                                    Double_t cut2,Double_t cut3,Double_t cut4,
1389                                    Double_t cut5,Double_t cut6,
1390                                    Double_t cut7,Double_t cut8)
1391 {
1392   // Set the cuts for J/psi from B selection
1393   fBtoJPSICuts[0] = cut0;
1394   fBtoJPSICuts[1] = cut1;
1395   fBtoJPSICuts[2] = cut2;
1396   fBtoJPSICuts[3] = cut3;
1397   fBtoJPSICuts[4] = cut4;
1398   fBtoJPSICuts[5] = cut5;
1399   fBtoJPSICuts[6] = cut6;
1400   fBtoJPSICuts[7] = cut7;
1401   fBtoJPSICuts[8] = cut8;
1402
1403   return;
1404 }
1405 //-----------------------------------------------------------------------------
1406 void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9])
1407 {
1408   // Set the cuts for J/psi from B selection
1409
1410   for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i];
1411
1412   return;
1413 }
1414 //-----------------------------------------------------------------------------
1415 void AliAnalysisVertexingHF::SetDplusCuts(Double_t cut0,Double_t cut1,
1416                                    Double_t cut2,Double_t cut3,Double_t cut4,
1417                                    Double_t cut5,Double_t cut6,
1418                                    Double_t cut7,Double_t cut8,
1419                                    Double_t cut9,Double_t cut10,Double_t cut11)
1420 {
1421   // Set the cuts for Dplus->Kpipi selection
1422   fDplusCuts[0] = cut0;
1423   fDplusCuts[1] = cut1;
1424   fDplusCuts[2] = cut2;
1425   fDplusCuts[3] = cut3;
1426   fDplusCuts[4] = cut4;
1427   fDplusCuts[5] = cut5;
1428   fDplusCuts[6] = cut6;
1429   fDplusCuts[7] = cut7;
1430   fDplusCuts[8] = cut8;
1431   fDplusCuts[9] = cut9;
1432   fDplusCuts[10] = cut10;
1433   fDplusCuts[11] = cut11;
1434
1435   return;
1436 }
1437 //-----------------------------------------------------------------------------
1438 void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12])
1439 {
1440   // Set the cuts for Dplus->Kpipi selection
1441
1442   for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i];
1443
1444   return;
1445 }
1446 //-----------------------------------------------------------------------------
1447 void AliAnalysisVertexingHF::SetDsCuts(Double_t cut0)
1448 {
1449   // Set the cuts for Ds->KKpi selection
1450   fDsCuts[0] = cut0;
1451
1452   return;
1453 }
1454 //-----------------------------------------------------------------------------
1455 void AliAnalysisVertexingHF::SetDsCuts(const Double_t cuts[1])
1456 {
1457   // Set the cuts for Ds->KKpi selection
1458
1459   for(Int_t i=0; i<1; i++) fDsCuts[i] = cuts[i];
1460
1461   return;
1462 }
1463 //-----------------------------------------------------------------------------
1464 void AliAnalysisVertexingHF::SetLcCuts(Double_t cut0)
1465 {
1466   // Set the cuts for Lc->pKpi selection
1467   fLcCuts[0] = cut0;
1468
1469   return;
1470 }
1471 //-----------------------------------------------------------------------------
1472 void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[1])
1473 {
1474   // Set the cuts for Lc->pKpi selection
1475
1476   for(Int_t i=0; i<1; i++) fLcCuts[i] = cuts[i];
1477
1478   return;
1479 }
1480 //-----------------------------------------------------------------------------
1481 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack& trk) const
1482 {
1483   // Check if track passes some kinematical cuts
1484
1485   if(TMath::Abs(trk.Pt()) < fMinPtCut) {
1486     //printf("pt %f\n",1./trk.GetParameter()[4]);
1487     return kFALSE;
1488   }
1489   trk.RelateToVertex(fV1,fBzkG,10.);
1490   Float_t d0z0[2],covd0z0[3];
1491   trk.GetImpactParameters(d0z0,covd0z0);
1492   if(TMath::Abs(d0z0[0]) < fMind0rphiCut) {
1493     printf("d0rphi %f\n",TMath::Abs(d0z0[0]));
1494     return kFALSE;
1495   }
1496
1497   return kTRUE;
1498 }
1499 //-----------------------------------------------------------------------------
1500 AliESDVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray) const
1501 {
1502   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1503
1504   AliESDVertex *vertex = 0;
1505
1506   if(!fSecVtxWithKF) { // AliVertexerTracks
1507
1508     AliVertexerTracks *vertexer2 = new AliVertexerTracks(fBzkG);
1509     vertexer2->SetDebug(0);
1510     vertexer2->SetVtxStart(fV1);
1511     vertex = (AliESDVertex*)vertexer2->VertexForSelectedESDTracks(trkArray);
1512     delete vertexer2;
1513
1514     if(vertex->GetNContributors()!=trkArray->GetEntriesFast()) {
1515       if(fDebug) printf("vertexing failed\n");
1516       delete vertex; vertex=0;
1517     }
1518
1519   } else { // Kalman Filter vertexer (AliKFParticle)
1520
1521     AliKFParticle::SetField(fBzkG);
1522
1523     AliKFParticle vertexKF;
1524
1525     Int_t nTrks = trkArray->GetEntriesFast();
1526     for(Int_t i=0; i<nTrks; i++) {
1527       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1528       AliKFParticle daughterKF(*esdTrack,211);
1529       vertexKF.AddDaughter(daughterKF);
1530     }
1531     vertex = new AliESDVertex();
1532     vertexKF.CopyToESDVertex(*vertex);
1533
1534   }
1535
1536   return vertex;
1537 }
1538 //-----------------------------------------------------------------------------
1539
1540
1541