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