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