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