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 // class AliAnalysisTaskVertexingHF.
  22 //
  23 //  Origin: E.Bruna, G.E.Bruno, A.Dainese, F.Prino, R.Romita
  24 //----------------------------------------------------------------------------
  25 #include <TTree.h>
  26 #include <TFile.h>
  27 #include <TDatabasePDG.h>
  28 #include <TString.h>
  29 #include "AliESDEvent.h"
  30 #include "AliVertexerTracks.h"
  31 #include "AliKFParticle.h"
  32 #include "AliESDVertex.h"
  33 #include "AliESDv0.h"
  34 #include "AliAODRecoDecay.h"
  35 #include "AliAODRecoDecayHF.h"
  36 #include "AliAODRecoDecayHF2Prong.h"
  37 #include "AliAODRecoDecayHF3Prong.h"
  38 #include "AliAODRecoDecayHF4Prong.h"
  39 #include "AliAnalysisVertexingHF.h"
  40
  41 ClassImp(AliAnalysisVertexingHF)
  42
  43 //----------------------------------------------------------------------------
  44 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
  45 fBzkG(0.),
  46 fSecVtxWithKF(kFALSE),
  47 fRecoPrimVtxSkippingTrks(kFALSE),
  48 fRmTrksFromPrimVtx(kFALSE),
  49 fV1(0x0),
  50 fDebug(0),
  51 fD0toKpi(kTRUE),
  52 fJPSItoEle(kTRUE),
  53 f3Prong(kTRUE),
  54 f4Prong(kTRUE),
  55 fITSrefit(kFALSE),
  56 fBothSPD(kTRUE),
  57 fMinITSCls(5),
  58 fMinPtCut(0.),
  59 fMind0rphiCut(0.)
  60 {
  61   // Default constructor
  62
  63   SetD0toKpiCuts();
  64   SetBtoJPSICuts();
  65   SetDplusCuts();
  66   SetDsCuts();
  67   SetLcCuts();
  68 }
  69 //--------------------------------------------------------------------------
  70 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
  71 TNamed(source),
  72 fBzkG(source.fBzkG),
  73 fSecVtxWithKF(source.fSecVtxWithKF),
  74 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
  75 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
  76 fV1(source.fV1),
  77 fDebug(source.fDebug),
  78 fD0toKpi(source.fD0toKpi),
  79 fJPSItoEle(source.fJPSItoEle),
  80 f3Prong(source.f3Prong),
  81 f4Prong(source.f4Prong),
  82 fITSrefit(source.fITSrefit),
  83 fBothSPD(source.fBothSPD),
  84 fMinITSCls(source.fMinITSCls),
  85 fMinPtCut(source.fMinPtCut),
  86 fMind0rphiCut(source.fMind0rphiCut)
  87 {
  88   //
  89   // Copy constructor
  90   //
  91   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
  92   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
  93   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
  94   for(Int_t i=0; i<1; i++)  fDsCuts[i]=source.fDsCuts[i];
  95   for(Int_t i=0; i<1; i++)  fLcCuts[i]=source.fLcCuts[i];
  96 }
  97 //--------------------------------------------------------------------------
  98 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
  99 {
 100   //
 101   // assignment operator
 102   //
 103   if(&source == this) return *this;
 104   fBzkG = source.fBzkG;
 105   fSecVtxWithKF = source.fSecVtxWithKF;
 106   fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
 107   fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
 108   fV1 = source.fV1;
 109   fDebug = source.fDebug;
 110   fD0toKpi = source.fD0toKpi;
 111   fJPSItoEle = source.fJPSItoEle;
 112   f3Prong = source.f3Prong;
 113   f4Prong = source.f4Prong;
 114   fITSrefit = source.fITSrefit;
 115   fBothSPD = source.fBothSPD;
 116   fMinITSCls = source.fMinITSCls;
 117   fMinPtCut = source.fMinPtCut;
 118   fMind0rphiCut = source.fMind0rphiCut;
 119
 120   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
 121   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
 122   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
 123   for(Int_t i=0; i<1; i++)  fDsCuts[i]=source.fDsCuts[i];
 124   for(Int_t i=0; i<1; i++)  fLcCuts[i]=source.fLcCuts[i];
 125
 126   return *this;
 127 }
 128 //----------------------------------------------------------------------------
 129 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
 130   // Destructor
 131   if(fV1) { delete fV1; fV1=0; }
 132 }
 133 //----------------------------------------------------------------------------
 134 void AliAnalysisVertexingHF::FindCandidates(AliESDEvent *esd,TTree treeout[])
 135 {
 136   // Find heavy-flavour vertex candidates
 137
 138
 139   AliAODRecoDecayHF2Prong *io2Prong = new AliAODRecoDecayHF2Prong();
 140   AliAODRecoDecayHF3Prong *io3Prong = new AliAODRecoDecayHF3Prong();
 141   AliAODRecoDecayHF4Prong *io4Prong = new AliAODRecoDecayHF4Prong();
 142   Int_t itree=0;
 143   Int_t itreeD0toKpi=-1,itreeJPSItoEle=-1,itree3Prong=-1,itree4Prong=-1;
 144   Int_t initEntriesD0toKpi=0,initEntriesJPSItoEle=0,initEntries3Prong=0,initEntries4Prong=0;
 145   if(fD0toKpi) {
 146     itreeD0toKpi=itree;
 147     treeout[itree].SetBranchAddress("D0toKpi",&io2Prong);
 148     itree++;
 149     initEntriesD0toKpi = treeout[itreeD0toKpi].GetEntries();
 150   }
 151   if(fJPSItoEle) {
 152     itreeJPSItoEle=itree;
 153     treeout[itree].SetBranchAddress("JPSItoEle",&io2Prong);
 154     itree++;
 155     initEntriesJPSItoEle = treeout[itreeJPSItoEle].GetEntries();
 156   }
 157   if(f3Prong) {
 158     itree3Prong=itree;
 159     treeout[itree].SetBranchAddress("Charmto3Prong",&io3Prong);
 160     itree++;
 161     initEntries3Prong = treeout[itree3Prong].GetEntries();
 162   }
 163   if(f4Prong) {
 164     itree4Prong=itree;
 165     treeout[itree].SetBranchAddress("D0to4Prong",&io4Prong);
 166     itree++;
 167     initEntries4Prong = treeout[itree4Prong].GetEntries();
 168   }
 169   delete io2Prong; io2Prong = NULL;
 170   delete io3Prong; io3Prong = NULL;
 171   delete io4Prong; io4Prong = NULL;
 172
 173   Int_t    iTrkP1,iTrkP2,iTrkN1,iTrkN2,trkEntries;
 174   Int_t    nTrksP=0,nTrksN=0;
 175   Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2;
 176   Bool_t   okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
 177   AliESDtrack *postrack1 = 0;
 178   AliESDtrack *postrack2 = 0;
 179   AliESDtrack *negtrack1 = 0;
 180   AliESDtrack *negtrack2 = 0;
 181   Double_t dcaMax = fD0toKpiCuts[1];
 182   if(dcaMax<fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1];
 183   if(dcaMax<fDplusCuts[11]) dcaMax=fDplusCuts[11];
 184   if(fDebug) printf(" dca cut set to %f cm\n",dcaMax);
 185
 186   Int_t ev = (Int_t)esd->GetEventNumberInFile();
 187   printf("--- Finding candidates in event %d\n",ev);
 188
 189   fBzkG = (Double_t)esd->GetMagneticField();
 190
 191   trkEntries = (Int_t)esd->GetNumberOfTracks();
 192   printf(" Number of tracks: %d\n",trkEntries);
 193   if(trkEntries<2) return;
 194
 195   // retrieve primary vertex from the AliESDEvent
 196   if(!esd->GetPrimaryVertex()) {
 197     printf(" No vertex in AliESD\n");
 198     return;
 199   }
 200   AliESDVertex copy(*(esd->GetPrimaryVertex()));
 201   SetPrimaryVertex(&copy);
 202
 203   // call function which applies sigle-track selection and
 204   // separetes positives and negatives
 205   TObjArray trksP(trkEntries/2); // will become TClonesArray
 206   TObjArray trksN(trkEntries/2); // will become TClonesArray
 207   SelectTracks(esd,trksP,nTrksP,
 208                    trksN,nTrksN);
 209
 210   printf(" Pos. tracks: %d    Neg. tracks: %d\n",nTrksP,nTrksN);
 211
 212   TObjArray *twoTrackArray1 = new TObjArray(2);
 213   TObjArray *twoTrackArray2 = new TObjArray(2);
 214   TObjArray *threeTrackArray = new TObjArray(3);
 215   TObjArray *fourTrackArray = new TObjArray(4);
 216
 217   // LOOP ON  POSITIVE  TRACKS
 218   for(iTrkP1=0; iTrkP1<nTrksP; iTrkP1++) {
 219     if(fDebug) if(iTrkP1%1==0) printf("  Processing positive track number %d of %d\n",iTrkP1,nTrksP);
 220     // get track from track array
 221     postrack1 = (AliESDtrack*)trksP.UncheckedAt(iTrkP1);
 222
 223     // LOOP ON  NEGATIVE  TRACKS
 224     for(iTrkN1=0; iTrkN1<nTrksN; iTrkN1++) {
 225       if(fDebug) if(iTrkN1%1==0) printf("    Processing negative track number %d of %d\n",iTrkN1,nTrksN);
 226       // get track from tracks array
 227       negtrack1 = (AliESDtrack*)trksN.UncheckedAt(iTrkN1);
 228       // DCA between the two tracks
 229       dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 230       if(dcap1n1>dcaMax) { negtrack1=0; continue; }
 231       // Vertexing
 232       twoTrackArray1->AddAt(postrack1,0);
 233       twoTrackArray1->AddAt(negtrack1,1);
 234       AliESDVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1);
 235       if(vertexp1n1->GetNContributors()!=2) {
 236         if(fDebug) printf("two-track vertexing failed\n");
 237         negtrack1=0; continue;
 238       }
 239       if(fD0toKpi || fJPSItoEle) {
 240         io2Prong = Make2Prong(twoTrackArray1,esd,vertexp1n1,dcap1n1,okD0,okJPSI);
 241         if(okD0)   treeout[itreeD0toKpi].Fill();
 242         if(okJPSI) treeout[itreeJPSItoEle].Fill();
 243         delete io2Prong; io2Prong=NULL;
 244       }
 245
 246       twoTrackArray1->Clear();
 247       if(!f3Prong && !f4Prong)  {
 248         negtrack1=0;
 249         delete vertexp1n1;
 250         continue;
 251       }
 252
 253       // 2nd LOOP  ON  POSITIVE  TRACKS
 254       for(iTrkP2=iTrkP1+1; iTrkP2<nTrksP; iTrkP2++) {
 255         // get track from tracks array
 256         postrack2 = (AliESDtrack*)trksP.UncheckedAt(iTrkP2);
 257         dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 258         if(dcap2n1>dcaMax) { postrack2=0; continue; }
 259         dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
 260         if(dcap1p2>dcaMax) { postrack2=0; continue; }
 261
 262         // Vertexing
 263         twoTrackArray2->AddAt(postrack2,0);
 264         twoTrackArray2->AddAt(negtrack1,1);
 265         AliESDVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2);
 266         if(f3Prong) {
 267           threeTrackArray->AddAt(postrack1,0);
 268           threeTrackArray->AddAt(negtrack1,1);
 269           threeTrackArray->AddAt(postrack2,2);
 270           io3Prong = Make3Prong(threeTrackArray,esd,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
 271           if(ok3Prong) treeout[itree3Prong].Fill();
 272           if(io3Prong) delete io3Prong; io3Prong=NULL;
 273         }
 274         if(f4Prong) {
 275           // 3rd LOOP  ON  NEGATIVE  TRACKS (for 4 prong)
 276           for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) {
 277             // get track from tracks array
 278             negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2);
 279             dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 280             if(dcap1n2>dcaMax) { negtrack2=0; continue; }
 281             // Vertexing
 282             fourTrackArray->AddAt(postrack1,0);
 283             fourTrackArray->AddAt(negtrack1,1);
 284             fourTrackArray->AddAt(postrack2,2);
 285             fourTrackArray->AddAt(negtrack2,3);
 286             io4Prong = Make4Prong(fourTrackArray,esd,vertexp1n1,vertexp2n1,dcap1n1,dcap1n2,dcap2n1,ok4Prong);
 287             if(ok4Prong) treeout[itree4Prong].Fill();
 288             delete io4Prong; io4Prong=NULL;
 289             fourTrackArray->Clear();
 290             negtrack2 = 0;
 291           } // end loop on negative tracks
 292         }
 293         postrack2 = 0;
 294         delete vertexp2n1;
 295       } // end 2nd loop on positive tracks
 296       twoTrackArray2->Clear();
 297
 298       // 2nd LOOP  ON  NEGATIVE  TRACKS
 299       for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) {
 300         // get track from tracks array
 301         negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2);
 302         dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 303         if(dcap1n2>dcaMax) { negtrack2=0; continue; }
 304         dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 305         if(dcan1n2>dcaMax) { negtrack2=0; continue; }
 306
 307         // Vertexing
 308         twoTrackArray2->AddAt(postrack1,0);
 309         twoTrackArray2->AddAt(negtrack2,1);
 310         AliESDVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2);
 311         if(f3Prong) {
 312           threeTrackArray->AddAt(negtrack1,0);
 313           threeTrackArray->AddAt(postrack1,1);
 314           threeTrackArray->AddAt(negtrack2,2);
 315           io3Prong = Make3Prong(threeTrackArray,esd,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
 316           if(ok3Prong) treeout[itree3Prong].Fill();
 317           if(io3Prong) delete io3Prong; io3Prong=NULL;
 318           }
 319         negtrack2 = 0;
 320         delete vertexp1n2;
 321       } // end 2nd loop on negative tracks
 322       twoTrackArray2->Clear();
 323
 324       negtrack1 = 0;
 325       delete vertexp1n1;
 326     } // end 1st loop on negative tracks
 327
 328     postrack1 = 0;
 329   }  // end 1st loop on positive tracks
 330
 331
 332
 333   //printf("delete twoTr 1\n");
 334   twoTrackArray1->Delete(); delete twoTrackArray1;
 335   //printf("delete twoTr 2\n");
 336   twoTrackArray2->Delete(); delete twoTrackArray2;
 337   //printf("delete threeTr 1\n");
 338   threeTrackArray->Clear();
 339   threeTrackArray->Delete(); delete threeTrackArray;
 340   //printf("delete fourTr 1\n");
 341   fourTrackArray->Delete(); delete fourTrackArray;
 342
 343
 344   // create a copy of this class to be written to output file
 345   //AliAnalysisVertexingHF *copy = (AliAnalysisVertexingHF*)this->Clone("AnalysisVertexingHF");
 346
 347   // print statistics
 348   if(fD0toKpi) {
 349     printf(" D0->Kpi: event %d = %d; total = %d;\n",
 350            (Int_t)esd->GetEventNumberInFile(),
 351            (Int_t)treeout[itreeD0toKpi].GetEntries()-initEntriesD0toKpi,
 352            (Int_t)treeout[itreeD0toKpi].GetEntries());
 353   }
 354   if(fJPSItoEle) {
 355     printf(" JPSI->ee: event %d = %d; total = %d;\n",
 356            (Int_t)esd->GetEventNumberInFile(),
 357            (Int_t)treeout[itreeJPSItoEle].GetEntries()-initEntriesJPSItoEle,
 358            (Int_t)treeout[itreeJPSItoEle].GetEntries());
 359   }
 360   if(f3Prong) {
 361     printf(" Charm->3Prong: event %d = %d; total = %d;\n",
 362            (Int_t)esd->GetEventNumberInFile(),
 363            (Int_t)treeout[itree3Prong].GetEntries()-initEntries3Prong,
 364            (Int_t)treeout[itree3Prong].GetEntries());
 365   }
 366   if(f4Prong) {
 367     printf(" Charm->4Prong: event %d = %d; total = %d;\n",
 368            (Int_t)esd->GetEventNumberInFile(),
 369            (Int_t)treeout[itree4Prong].GetEntries()-initEntries4Prong,
 370            (Int_t)treeout[itree4Prong].GetEntries());
 371   }
 372
 373
 374   return;
 375 }
 376 //----------------------------------------------------------------------------
 377 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
 378                                    TObjArray *twoTrackArray1,AliESDEvent *esd,
 379                                    AliESDVertex *secVertexESD,Double_t dca,
 380                                    Bool_t &okD0,Bool_t &okJPSI) const
 381 {
 382   // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
 383   // reconstruction cuts
 384   // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
 385
 386   okD0=kFALSE; okJPSI=kFALSE;
 387
 388   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
 389
 390   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(0);
 391   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(1);
 392
 393   // propagate tracks to secondary vertex, to compute inv. mass
 394   postrack->RelateToVertex(secVertexESD,fBzkG,10.);
 395   negtrack->RelateToVertex(secVertexESD,fBzkG,10.);
 396
 397   Double_t momentum[3];
 398   postrack->GetPxPyPz(momentum);
 399   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 400   negtrack->GetPxPyPz(momentum);
 401   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 402
 403
 404   // invariant mass cut (try to improve coding here..)
 405   Bool_t okMassCut=kFALSE;
 406   if(!okMassCut && fD0toKpi) if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
 407   if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
 408   if(!okMassCut) {
 409     if(fDebug) printf(" candidate didn't pass mass cut\n");
 410     return 0x0;
 411   }
 412
 413
 414   AliESDVertex *primVertex = fV1;
 415   AliESDVertex *ownPrimVertex=0;
 416
 417   // primary vertex from *other* tracks in the event
 418   if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
 419     ownPrimVertex = OwnPrimaryVertex(2,twoTrackArray1,esd);
 420     if(!ownPrimVertex) {
 421       return 0x0;
 422     } else {
 423       if(ownPrimVertex->GetNContributors()<2) {
 424         delete ownPrimVertex;
 425         return 0x0;
 426       } else {
 427         primVertex = ownPrimVertex;
 428       }
 429     }
 430   }
 431
 432   Float_t d0z0[2],covd0z0[3];
 433   postrack->RelateToVertex(primVertex,fBzkG,10.);
 434   postrack->GetImpactParameters(d0z0,covd0z0);
 435   d0[0] = d0z0[0];
 436   d0err[0] = TMath::Sqrt(covd0z0[0]);
 437   negtrack->RelateToVertex(primVertex,fBzkG,10.);
 438   negtrack->GetImpactParameters(d0z0,covd0z0);
 439   d0[1] = d0z0[0];
 440   d0err[1] = TMath::Sqrt(covd0z0[0]);
 441
 442   // create the object AliAODRecoDecayHF2Prong
 443   Double_t pos[3],cov[6];
 444   secVertexESD->GetXYZ(pos); // position
 445   secVertexESD->GetCovMatrix(cov); //covariance matrix
 446   AliAODVertex *secVertexAOD = new AliAODVertex(pos,cov,secVertexESD->GetChi2toNDF());
 447   primVertex->GetXYZ(pos); // position
 448   primVertex->GetCovMatrix(cov); //covariance matrix
 449   AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF());
 450   AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVertexAOD,px,py,pz,d0,d0err,dca);
 451   the2Prong->SetOwnSecondaryVtx(secVertexAOD);//temporary
 452   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
 453
 454   // select D0->Kpi
 455   Int_t checkD0,checkD0bar;
 456   if(fD0toKpi) okD0 = the2Prong->SelectD0(fD0toKpiCuts,checkD0,checkD0bar);
 457   //if(fDebug && fD0toKpi) printf("   %d\n",(Int_t)okD0);
 458   // select J/psi from B
 459   Int_t checkJPSI;
 460   if(fJPSItoEle) okJPSI = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI);
 461   //if(fDebug && fJPSItoEle) printf("   %d\n",(Int_t)okJPSI);
 462
 463
 464   if(okD0 || okJPSI) {
 465     // get PID info from ESD
 466     Double_t esdpid0[5];
 467     postrack->GetESDpid(esdpid0);
 468     Double_t esdpid1[5];
 469     negtrack->GetESDpid(esdpid1);
 470     Double_t esdpid[10];
 471     for(Int_t i=0;i<5;i++) {
 472       esdpid[i]   = esdpid0[i];
 473       esdpid[5+i] = esdpid1[i];
 474     }
 475     the2Prong->SetPID(2,esdpid);
 476   }
 477
 478   if(ownPrimVertex) delete ownPrimVertex;
 479
 480   return the2Prong;
 481 }
 482 //----------------------------------------------------------------------------
 483 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
 484                              TObjArray *threeTrackArray,AliESDEvent *esd,
 485                              AliESDVertex *vertexp1n1,AliESDVertex *vertexp2n1,
 486                              Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
 487                              Bool_t &ok3Prong) const
 488 {
 489   // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
 490   // reconstruction cuts
 491   // E.Bruna, F.Prino
 492
 493   ok3Prong=kFALSE;
 494   Double_t px[3],py[3],pz[3],d0[3],d0err[3];//d0z[3];
 495   Float_t d0z0[2],covd0z0[3];
 496
 497
 498   AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
 499   AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
 500   AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
 501
 502   AliESDVertex *primVertex = fV1;
 503
 504   postrack1->RelateToVertex(primVertex,fBzkG,10.);
 505   negtrack->RelateToVertex(primVertex,fBzkG,10.);
 506   postrack2->RelateToVertex(primVertex,fBzkG,10.);
 507
 508   Double_t momentum[3];
 509   postrack1->GetPxPyPz(momentum);
 510   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 511   negtrack->GetPxPyPz(momentum);
 512   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 513   postrack2->GetPxPyPz(momentum);
 514   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
 515
 516   postrack1->GetImpactParameters(d0z0,covd0z0);
 517   d0[0]=d0z0[0];
 518   d0err[0] = TMath::Sqrt(covd0z0[0]);
 519   negtrack->GetImpactParameters(d0z0,covd0z0);
 520   d0[1]=d0z0[0];
 521   d0err[1] = TMath::Sqrt(covd0z0[0]);
 522   postrack2->GetImpactParameters(d0z0,covd0z0);
 523   d0[2]=d0z0[0];
 524   d0err[2] = TMath::Sqrt(covd0z0[0]);
 525
 526
 527   // invariant mass cut for D+, Ds, Lc
 528   Bool_t okMassCut=kFALSE;
 529   if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
 530   if(!okMassCut) {
 531     if(fDebug) printf(" candidate didn't pass mass cut\n");
 532     return 0x0;
 533   }
 534
 535   //charge
 536   Short_t charge=(Short_t)(postrack1->GetSign()*postrack2->GetSign()*negtrack->GetSign());
 537
 538   AliESDVertex *ownPrimVertex = 0;
 539   // primary vertex from *other* tracks in the event
 540   if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
 541     ownPrimVertex = OwnPrimaryVertex(3,threeTrackArray,esd);
 542     if(!ownPrimVertex) {
 543       return 0x0;
 544     } else {
 545       if(ownPrimVertex->GetNContributors()<2) {
 546         delete ownPrimVertex;
 547         return 0x0;
 548       } else {
 549         primVertex = ownPrimVertex;
 550       }
 551     }
 552   }
 553
 554   // create the object AliAODRecoDecayHF3Prong
 555   AliESDVertex* secVert3Prong = ReconstructSecondaryVertex(threeTrackArray);
 556   Double_t pos[3],cov[6],sigmavert;
 557   secVert3Prong->GetXYZ(pos); // position
 558   secVert3Prong->GetCovMatrix(cov); //covariance matrix
 559   sigmavert=secVert3Prong->GetDispersion();
 560
 561   AliAODVertex *secVert3PrAOD = new AliAODVertex(pos,cov,secVert3Prong->GetChi2toNDF());
 562   primVertex->GetXYZ(pos); // position
 563   primVertex->GetCovMatrix(cov); //covariance matrix
 564   AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF());
 565   Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
 566
 567   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]));
 568   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]));
 569
 570   AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert3PrAOD,px,py,pz,d0,d0err,dca,sigmavert,dist12,dist23,charge);
 571   the3Prong->SetOwnPrimaryVtx(primVertexAOD);
 572   the3Prong->SetOwnSecondaryVtx(secVert3PrAOD);//temporary
 573
 574
 575   // select D+->Kpipi, Ds->KKpi, Lc->pKpi
 576   if(f3Prong) {
 577     ok3Prong = kFALSE;
 578     Int_t ok1,ok2;
 579     if(the3Prong->SelectDplus(fDplusCuts))   ok3Prong = kTRUE;
 580     if(the3Prong->SelectDs(fDsCuts,ok1,ok2)) ok3Prong = kTRUE;
 581     if(the3Prong->SelectLc(fLcCuts,ok1,ok2)) ok3Prong = kTRUE;
 582   }
 583   //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
 584   if(ok3Prong) {
 585     // get PID info from ESD
 586     Double_t esdpid0[5];
 587     postrack1->GetESDpid(esdpid0);
 588     Double_t esdpid1[5];
 589     negtrack->GetESDpid(esdpid1);
 590     Double_t esdpid2[5];
 591     postrack2->GetESDpid(esdpid2);
 592
 593
 594     Double_t esdpid[15];
 595     for(Int_t i=0;i<5;i++) {
 596       esdpid[i]   = esdpid0[i];
 597       esdpid[5+i] = esdpid1[i];
 598       esdpid[10+i] = esdpid2[i];
 599     }
 600     the3Prong->SetPID(3,esdpid);
 601   }
 602
 603   if(ownPrimVertex) delete ownPrimVertex;
 604
 605   return the3Prong;
 606 }
 607 //----------------------------------------------------------------------------
 608 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
 609                              TObjArray *fourTrackArray,AliESDEvent *esd,
 610                              AliESDVertex *vertexp1n1,AliESDVertex *vertexp2n1,
 611                              Double_t dcap1n1,Double_t dcap1n2,Double_t dcap2n1,
 612                              Bool_t &ok4Prong) const
 613 {
 614   // Make 4Prong candidates and check if they pass D0toKpipipi
 615   // reconstruction cuts
 616   // G.E.Bruno, R.Romita
 617
 618   ok4Prong=kFALSE;
 619
 620   Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
 621   //Float_t d0z0[2],covd0z0[3];
 622
 623   px[0]=dcap1n1*dcap1n2*dcap2n1; // TO BE CHANGED (done just to removed compilation warning about dca... not used)
 624
 625   //charge
 626   Short_t charge=0;
 627
 628   AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
 629   AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
 630   AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
 631   AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
 632
 633   AliESDVertex *primVertex = fV1;
 634
 635   postrack1->RelateToVertex(primVertex,fBzkG,10.);
 636   negtrack1->RelateToVertex(primVertex,fBzkG,10.);
 637   postrack2->RelateToVertex(primVertex,fBzkG,10.);
 638   negtrack2->RelateToVertex(primVertex,fBzkG,10.);
 639
 640   Double_t momentum[3];
 641   postrack1->GetPxPyPz(momentum);
 642   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 643   negtrack1->GetPxPyPz(momentum);
 644   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 645   postrack2->GetPxPyPz(momentum);
 646   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
 647   negtrack2->GetPxPyPz(momentum);
 648   px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
 649
 650   // invariant mass cut for rho or D0 (try to improve coding here..)
 651   //Bool_t okMassCut=kFALSE;
 652   //if(!okMassCut) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
 653   //if(!okMassCut) {
 654   //  if(fDebug) printf(" candidate didn't pass mass cut\n");
 655   //  return 0x0;
 656   //}
 657
 658   AliESDVertex *ownPrimVertex = 0;
 659   // primary vertex from *other* tracks in the event
 660   if(fRecoPrimVtxSkippingTrks || fRmTrksFromPrimVtx) {
 661     ownPrimVertex = OwnPrimaryVertex(4,fourTrackArray,esd);
 662     if(!ownPrimVertex) {
 663       return 0x0;
 664     } else {
 665       if(ownPrimVertex->GetNContributors()<2) {
 666         delete ownPrimVertex;
 667         return 0x0;
 668       } else {
 669         primVertex = ownPrimVertex;
 670       }
 671     }
 672   }
 673
 674   // create the object AliAODRecoDecayHF4Prong
 675   AliESDVertex* secVert4Prong = ReconstructSecondaryVertex(fourTrackArray);
 676   Double_t pos[3],cov[6],sigmavert;
 677   secVert4Prong->GetXYZ(pos); // position
 678   secVert4Prong->GetCovMatrix(cov); //covariance matrix
 679   sigmavert=secVert4Prong->GetDispersion();
 680
 681   AliAODVertex *secVert4PrAOD = new AliAODVertex(pos,cov,secVert4Prong->GetChi2toNDF());
 682   primVertex->GetXYZ(pos); // position
 683   primVertex->GetCovMatrix(cov); //covariance matrix
 684   AliAODVertex *primVertexAOD = new AliAODVertex(pos,cov,primVertex->GetChi2toNDF());
 685   //Double_t dca[6]={dcap1n1,dcap2n1,dcap1p2,0.,0.,0.}; //
 686   Double_t dca[6]={0.,0.,0.,0.,0.,0.}; //  modify it
 687
 688   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]));
 689   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]));
 690   Double_t dist14=0.; // to be implemented
 691   Double_t dist34=0.; // to be implemented
 692
 693   //AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert4PrAOD,px,py,pz,d0,d0err,dca,sigmavert,dist12,dist23,charge);
 694   AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert4PrAOD,px,py,pz,d0,d0err,dca,dist12,dist23,dist14,dist34,charge);
 695   the4Prong->SetOwnPrimaryVtx(primVertexAOD);
 696   the4Prong->SetOwnSecondaryVtx(secVert4PrAOD);//temporary
 697
 698
 699   // use the following two lines once AliAODRecoDecayHF4Prong::SelectD0 is available
 700   // select D0->Kpipipi
 701   //Int_t checkD0,checkD0bar;
 702   // ok4Prong=the4Prong->SelectD0(fD04pCuts,checkD0,checkD0bar);
 703   ok4Prong=kFALSE;  //for the time being ...
 704
 705
 706   // get PID info from ESD
 707   Double_t esdpid0[5];
 708   postrack1->GetESDpid(esdpid0);
 709   Double_t esdpid1[5];
 710   negtrack1->GetESDpid(esdpid1);
 711   Double_t esdpid2[5];
 712   postrack2->GetESDpid(esdpid2);
 713   Double_t esdpid3[5];
 714   negtrack2->GetESDpid(esdpid3);
 715
 716   Double_t esdpid[20];
 717   for(Int_t i=0;i<5;i++) {
 718     esdpid[i]   = esdpid0[i];
 719     esdpid[5+i] = esdpid1[i];
 720     esdpid[10+i] = esdpid2[i];
 721     esdpid[15+i] = esdpid3[i];
 722   }
 723   the4Prong->SetPID(4,esdpid);
 724
 725   if(ownPrimVertex) delete ownPrimVertex;
 726
 727   return the4Prong;
 728 }
 729 //-----------------------------------------------------------------------------
 730 AliESDVertex* AliAnalysisVertexingHF::OwnPrimaryVertex(Int_t ntrks,
 731                                                        TObjArray *trkArray,
 732                                                        AliESDEvent *esd) const
 733 {
 734   // Returns primary vertex specific to this candidate
 735
 736   AliVertexerTracks *vertexer1 = new AliVertexerTracks(esd->GetMagneticField());
 737   AliESDVertex *ownPrimVertex = 0;
 738
 739   // recalculating the vertex
 740   if(fRecoPrimVtxSkippingTrks) {
 741     if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
 742       Float_t diamondcovxy[3];
 743       esd->GetDiamondCovXY(diamondcovxy);
 744       Double_t pos[3]={esd->GetDiamondX(),esd->GetDiamondY(),0.};
 745       Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.};
 746       AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
 747       vertexer1->SetVtxStart(diamond);
 748       delete diamond; diamond=NULL;
 749       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
 750         vertexer1->SetOnlyFitter();
 751     }
 752     Int_t skipped[10];
 753     AliESDtrack *t = 0;
 754     for(Int_t i=0; i<ntrks; i++) {
 755       t = (AliESDtrack*)trkArray->UncheckedAt(i);
 756       skipped[i] = (Int_t)t->GetID();
 757     }
 758     vertexer1->SetSkipTracks(ntrks,skipped);
 759     ownPrimVertex = (AliESDVertex*)vertexer1->FindPrimaryVertex(esd);
 760   }
 761
 762   // removing the prongs tracks
 763   if(fRmTrksFromPrimVtx) {
 764     TTree *rmTree = new TTree("trksTree","tree with tracks to be removed");
 765     AliESDtrack *esdTrack = 0;
 766     rmTree->Branch("tracks","AliESDtrack",&esdTrack);
 767     AliESDtrack *t = 0;
 768     for(Int_t i=0; i<ntrks; i++) {
 769       t = (AliESDtrack*)trkArray->UncheckedAt(i);
 770       esdTrack = new AliESDtrack(*t);
 771       rmTree->Fill();
 772       delete esdTrack;
 773     }
 774     Float_t diamondxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
 775     ownPrimVertex = vertexer1->RemoveTracksFromVertex(fV1,rmTree,diamondxy);
 776     delete rmTree; rmTree=NULL;
 777   }
 778
 779   delete vertexer1; vertexer1=NULL;
 780
 781   return ownPrimVertex;
 782 }
 783 //-----------------------------------------------------------------------------
 784 void AliAnalysisVertexingHF::PrintStatus() const {
 785   // Print parameters being used
 786
 787   printf("Preselections:\n");
 788   printf("    fITSrefit   = %d\n",(Int_t)fITSrefit);
 789   printf("    fBothSPD   = %d\n",(Int_t)fBothSPD);
 790   printf("    fMinITSCls   = %d\n",fMinITSCls);
 791   printf("    fMinPtCut   = %f GeV/c\n",fMinPtCut);
 792   printf("    fMind0rphiCut   = %f cm\n",fMind0rphiCut);
 793   if(fSecVtxWithKF) {
 794     printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
 795   } else {
 796     printf("Secondary vertex with AliVertexerTracks\n");
 797   }
 798   if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
 799   if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
 800   if(fD0toKpi) {
 801     printf("Reconstruct D0->Kpi candidates with cuts:\n");
 802     printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
 803     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
 804     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
 805     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
 806     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
 807     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
 808     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
 809     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
 810     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
 811   }
 812   if(fJPSItoEle) {
 813     printf("Reconstruct J/psi from B candidates with cuts:\n");
 814     printf("    |M-MJPSI| [GeV]    < %f\n",fBtoJPSICuts[0]);
 815     printf("    dca    [cm]  < %f\n",fBtoJPSICuts[1]);
 816     printf("    cosThetaStar     < %f\n",fBtoJPSICuts[2]);
 817     printf("    pTP     [GeV/c]    > %f\n",fBtoJPSICuts[3]);
 818     printf("    pTN    [GeV/c]    > %f\n",fBtoJPSICuts[4]);
 819     printf("    |d0P|  [cm]  < %f\n",fBtoJPSICuts[5]);
 820     printf("    |d0N| [cm]  < %f\n",fBtoJPSICuts[6]);
 821     printf("    d0d0  [cm^2] < %f\n",fBtoJPSICuts[7]);
 822     printf("    cosThetaPoint    > %f\n",fBtoJPSICuts[8]);
 823   }
 824   if(f3Prong) {
 825     printf("Reconstruct 3 prong candidates.\n");
 826     printf("  D+->Kpipi cuts:\n");
 827     printf("    |M-MD+| [GeV]    < %f\n",fDplusCuts[0]);
 828     printf("    pTK     [GeV/c]    > %f\n",fDplusCuts[1]);
 829     printf("    pTPi    [GeV/c]    > %f\n",fDplusCuts[2]);
 830     printf("    |d0K|  [cm]  > %f\n",fDplusCuts[3]);
 831     printf("    |d0Pi| [cm]  > %f\n",fDplusCuts[4]);
 832     printf("    dist12    [cm]  < %f\n",fDplusCuts[5]);
 833     printf("    sigmavert [cm]   < %f\n",fDplusCuts[6]);
 834     printf("    dist prim-sec [cm] > %f\n",fDplusCuts[7]);
 835     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDplusCuts[8]);
 836     printf("    cosThetaPoint    > %f\n",fDplusCuts[9]);
 837     printf("    Sum d0^2 [cm^2]  > %f\n",fDplusCuts[10]);
 838     printf("    dca cut [cm]  < %f\n",fDplusCuts[11]);
 839     printf("  Ds->KKpi cuts:\n");
 840     printf("    |M-MDs| [GeV]    < %f\n",fDsCuts[0]);
 841     printf("  Lc->pKpi cuts:\n");
 842     printf("    |M-MD+| [GeV]    < %f\n",fLcCuts[0]);
 843   }
 844
 845   return;
 846 }
 847 //-----------------------------------------------------------------------------
 848 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
 849                                              Int_t nprongs,
 850                                              Double_t *px,
 851                                              Double_t *py,
 852                                              Double_t *pz) const {
 853   // Check invariant mass cut
 854
 855   Short_t dummycharge=0;
 856   Double_t *dummyd0 = new Double_t[nprongs];
 857   AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
 858   delete [] dummyd0;
 859
 860   UInt_t pdg2[2],pdg3[3];
 861   Double_t mPDG,minv;
 862
 863   Bool_t retval=kFALSE;
 864   switch (decay)
 865     {
 866     case 0:                  // D0->Kpi
 867       pdg2[0]=211; pdg2[1]=321;
 868       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 869       minv = rd->InvMass(nprongs,pdg2);
 870       if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
 871       pdg2[0]=321; pdg2[1]=211;
 872       minv = rd->InvMass(nprongs,pdg2);
 873       if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
 874       break;
 875     case 1:                  // JPSI->ee
 876       pdg2[0]=11; pdg2[1]=11;
 877       mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
 878       minv = rd->InvMass(nprongs,pdg2);
 879       if(TMath::Abs(minv-mPDG)<fBtoJPSICuts[0]) retval=kTRUE;
 880       break;
 881     case 2:                  // D+->Kpipi
 882       pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
 883       mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
 884       minv = rd->InvMass(nprongs,pdg3);
 885       if(TMath::Abs(minv-mPDG)<fDplusCuts[0]) retval=kTRUE;
 886                             // Ds+->KKpi
 887       pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
 888       mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
 889       minv = rd->InvMass(nprongs,pdg3);
 890       if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
 891       pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
 892       minv = rd->InvMass(nprongs,pdg3);
 893       if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
 894                             // Lc->pKpi
 895       pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
 896       mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
 897       minv = rd->InvMass(nprongs,pdg3);
 898       if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
 899       pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
 900       minv = rd->InvMass(nprongs,pdg3);
 901       if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
 902       break;
 903     default:
 904       break;
 905     }
 906
 907   delete rd;
 908
 909   return retval;
 910 }
 911 //-----------------------------------------------------------------------------
 912 void AliAnalysisVertexingHF::SelectTracks(AliESDEvent *esd,
 913                                           TObjArray &trksP,Int_t &nTrksP,
 914                                           TObjArray &trksN,Int_t &nTrksN) const
 915 {
 916   // Fill two TObjArrays with positive and negative tracks and
 917   // apply single-track preselection
 918
 919   nTrksP=0,nTrksN=0;
 920
 921   Int_t entries = (Int_t)esd->GetNumberOfTracks();
 922
 923   // transfer ITS tracks from ESD to arrays and to a tree
 924   for(Int_t i=0; i<entries; i++) {
 925
 926     AliESDtrack *esdtrack = esd->GetTrack(i);
 927     UInt_t status = esdtrack->GetStatus();
 928
 929     // require refit in ITS
 930     if(fITSrefit && !(status&AliESDtrack::kITSrefit)) {
 931       if(fDebug) printf("track %d is not kITSrefit\n",i);
 932       continue;
 933     }
 934
 935     // require minimum # of ITS points
 936     if(esdtrack->GetNcls(0)<fMinITSCls)  {
 937       if(fDebug) printf("track %d has %d ITS cls\n",i,esdtrack->GetNcls(0));
 938       continue;
 939     }
 940     // require points on the 2 pixel layers
 941     if(fBothSPD)
 942       if(!TESTBIT(esdtrack->GetITSClusterMap(),0) ||
 943          !TESTBIT(esdtrack->GetITSClusterMap(),1)) continue;
 944
 945     // single track selection
 946     if(!SingleTrkCuts(*esdtrack)) continue;
 947
 948     if(esdtrack->GetSign()<0) { // negative track
 949       trksN.AddLast(esdtrack);
 950       nTrksN++;
 951     } else {                 // positive track
 952       trksP.AddLast(esdtrack);
 953       nTrksP++;
 954     }
 955
 956   } // loop on ESD tracks
 957
 958   return;
 959 }
 960 //-----------------------------------------------------------------------------
 961 void AliAnalysisVertexingHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
 962                                    Double_t cut2,Double_t cut3,Double_t cut4,
 963                                    Double_t cut5,Double_t cut6,
 964                                    Double_t cut7,Double_t cut8)
 965 {
 966   // Set the cuts for D0 selection
 967   fD0toKpiCuts[0] = cut0;
 968   fD0toKpiCuts[1] = cut1;
 969   fD0toKpiCuts[2] = cut2;
 970   fD0toKpiCuts[3] = cut3;
 971   fD0toKpiCuts[4] = cut4;
 972   fD0toKpiCuts[5] = cut5;
 973   fD0toKpiCuts[6] = cut6;
 974   fD0toKpiCuts[7] = cut7;
 975   fD0toKpiCuts[8] = cut8;
 976
 977   return;
 978 }
 979 //-----------------------------------------------------------------------------
 980 void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9])
 981 {
 982   // Set the cuts for D0 selection
 983
 984   for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
 985
 986   return;
 987 }
 988 //-----------------------------------------------------------------------------
 989 void AliAnalysisVertexingHF::SetBtoJPSICuts(Double_t cut0,Double_t cut1,
 990                                    Double_t cut2,Double_t cut3,Double_t cut4,
 991                                    Double_t cut5,Double_t cut6,
 992                                    Double_t cut7,Double_t cut8)
 993 {
 994   // Set the cuts for J/psi from B selection
 995   fBtoJPSICuts[0] = cut0;
 996   fBtoJPSICuts[1] = cut1;
 997   fBtoJPSICuts[2] = cut2;
 998   fBtoJPSICuts[3] = cut3;
 999   fBtoJPSICuts[4] = cut4;
1000   fBtoJPSICuts[5] = cut5;
1001   fBtoJPSICuts[6] = cut6;
1002   fBtoJPSICuts[7] = cut7;
1003   fBtoJPSICuts[8] = cut8;
1004
1005   return;
1006 }
1007 //-----------------------------------------------------------------------------
1008 void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9])
1009 {
1010   // Set the cuts for J/psi from B selection
1011
1012   for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i];
1013
1014   return;
1015 }
1016 //-----------------------------------------------------------------------------
1017 void AliAnalysisVertexingHF::SetDplusCuts(Double_t cut0,Double_t cut1,
1018                                    Double_t cut2,Double_t cut3,Double_t cut4,
1019                                    Double_t cut5,Double_t cut6,
1020                                    Double_t cut7,Double_t cut8,
1021                                    Double_t cut9,Double_t cut10,Double_t cut11)
1022 {
1023   // Set the cuts for Dplus->Kpipi selection
1024   fDplusCuts[0] = cut0;
1025   fDplusCuts[1] = cut1;
1026   fDplusCuts[2] = cut2;
1027   fDplusCuts[3] = cut3;
1028   fDplusCuts[4] = cut4;
1029   fDplusCuts[5] = cut5;
1030   fDplusCuts[6] = cut6;
1031   fDplusCuts[7] = cut7;
1032   fDplusCuts[8] = cut8;
1033   fDplusCuts[9] = cut9;
1034   fDplusCuts[10] = cut10;
1035   fDplusCuts[11] = cut11;
1036
1037   return;
1038 }
1039 //-----------------------------------------------------------------------------
1040 void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12])
1041 {
1042   // Set the cuts for Dplus->Kpipi selection
1043
1044   for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i];
1045
1046   return;
1047 }
1048 //-----------------------------------------------------------------------------
1049 void AliAnalysisVertexingHF::SetDsCuts(Double_t cut0)
1050 {
1051   // Set the cuts for Ds->KKpi selection
1052   fDsCuts[0] = cut0;
1053
1054   return;
1055 }
1056 //-----------------------------------------------------------------------------
1057 void AliAnalysisVertexingHF::SetDsCuts(const Double_t cuts[1])
1058 {
1059   // Set the cuts for Ds->KKpi selection
1060
1061   for(Int_t i=0; i<1; i++) fDsCuts[i] = cuts[i];
1062
1063   return;
1064 }
1065 //-----------------------------------------------------------------------------
1066 void AliAnalysisVertexingHF::SetLcCuts(Double_t cut0)
1067 {
1068   // Set the cuts for Lc->pKpi selection
1069   fLcCuts[0] = cut0;
1070
1071   return;
1072 }
1073 //-----------------------------------------------------------------------------
1074 void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[1])
1075 {
1076   // Set the cuts for Lc->pKpi selection
1077
1078   for(Int_t i=0; i<1; i++) fLcCuts[i] = cuts[i];
1079
1080   return;
1081 }
1082 //-----------------------------------------------------------------------------
1083 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack& trk) const
1084 {
1085   // Check if track passes some kinematical cuts
1086
1087   if(TMath::Abs(1./trk.GetParameter()[4]) < fMinPtCut) {
1088     printf("pt %f\n",1./trk.GetParameter()[4]);
1089     return kFALSE;
1090   }
1091   trk.RelateToVertex(fV1,fBzkG,10.);
1092   Float_t d0z0[2],covd0z0[3];
1093   trk.GetImpactParameters(d0z0,covd0z0);
1094   if(TMath::Abs(d0z0[0]) < fMind0rphiCut) {
1095     printf("d0rphi %f\n",TMath::Abs(d0z0[0]));
1096     return kFALSE;
1097   }
1098
1099   return kTRUE;
1100 }
1101 //-----------------------------------------------------------------------------
1102 AliESDVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray) const
1103 {
1104   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1105
1106   AliESDVertex *vertex = 0;
1107
1108   if(!fSecVtxWithKF) { // AliVertexerTracks
1109
1110     AliVertexerTracks *vertexer2 = new AliVertexerTracks(fBzkG);
1111     vertexer2->SetDebug(0);
1112     vertexer2->SetVtxStart(fV1);
1113     vertex = (AliESDVertex*)vertexer2->VertexForSelectedTracks(trkArray);
1114     delete vertexer2;
1115
1116   } else { // Kalman Filter vertexer (AliKFParticle)
1117
1118     AliKFParticle::SetField(fBzkG);
1119
1120     AliKFParticle vertexKF;
1121
1122     Int_t nTrks = trkArray->GetEntriesFast();
1123     for(Int_t i=0; i<nTrks; i++) {
1124       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1125       AliKFParticle daughterKF(*esdTrack,211);
1126       vertexKF.AddDaughter(daughterKF);
1127     }
1128     vertex = new AliESDVertex();
1129     vertexKF.CopyToESDVertex(*vertex);
1130
1131   }
1132
1133   return vertex;
1134 }
1135 //-----------------------------------------------------------------------------
1136
1137
1138