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