PWG3/vertexingHF/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 stored in the AOD as objects deriving from AliAODRecoDecay.
  19 // To be used as a task of AliAnalysisManager by means of the interface
  20 // class AliAnalysisTaskSEVertexingHF.
  21 // An example of usage in the macro AliAnalysisTaskSEVertexingHFTest.C.
  22 //
  23 //  Origin: E.Bruna, G.E.Bruno, A.Dainese, F.Prino, R.Romita
  24 //  Contact: andrea.dainese@lnl.infn.it
  25 //----------------------------------------------------------------------------
  26 #include <TFile.h>
  27 #include <TDatabasePDG.h>
  28 #include <TString.h>
  29 #include "AliLog.h"
  30 #include "AliVEvent.h"
  31 #include "AliVVertex.h"
  32 #include "AliVTrack.h"
  33 #include "AliVertexerTracks.h"
  34 #include "AliKFParticle.h"
  35 #include "AliESDEvent.h"
  36 #include "AliESDVertex.h"
  37 #include "AliESDtrack.h"
  38 #include "AliAODEvent.h"
  39 #include "AliAODRecoDecay.h"
  40 #include "AliAODRecoDecayHF.h"
  41 #include "AliAODRecoDecayHF2Prong.h"
  42 #include "AliAODRecoDecayHF3Prong.h"
  43 #include "AliAODRecoDecayHF4Prong.h"
  44 #include "AliAnalysisVertexingHF.h"
  45
  46 ClassImp(AliAnalysisVertexingHF)
  47
  48 //----------------------------------------------------------------------------
  49 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
  50 fInputAOD(kFALSE),
  51 fAODMap(0),
  52 fBzkG(0.),
  53 fSecVtxWithKF(kFALSE),
  54 fRecoPrimVtxSkippingTrks(kFALSE),
  55 fRmTrksFromPrimVtx(kFALSE),
  56 fV1(0x0),
  57 fD0toKpi(kTRUE),
  58 fJPSItoEle(kTRUE),
  59 f3Prong(kTRUE),
  60 f4Prong(kTRUE),
  61 fITSrefit(kFALSE),
  62 fBothSPD(kTRUE),
  63 fMinITSCls(5),
  64 fMinPtCut(0.),
  65 fMind0rphiCut(0.)
  66 {
  67   // Default constructor
  68
  69   SetD0toKpiCuts();
  70   SetBtoJPSICuts();
  71   SetDplusCuts();
  72   SetDsCuts();
  73   SetLcCuts();
  74 }
  75 //--------------------------------------------------------------------------
  76 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
  77 TNamed(source),
  78 fInputAOD(source.fInputAOD),
  79 fAODMap(source.fAODMap),
  80 fBzkG(source.fBzkG),
  81 fSecVtxWithKF(source.fSecVtxWithKF),
  82 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
  83 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
  84 fV1(source.fV1),
  85 fD0toKpi(source.fD0toKpi),
  86 fJPSItoEle(source.fJPSItoEle),
  87 f3Prong(source.f3Prong),
  88 f4Prong(source.f4Prong),
  89 fITSrefit(source.fITSrefit),
  90 fBothSPD(source.fBothSPD),
  91 fMinITSCls(source.fMinITSCls),
  92 fMinPtCut(source.fMinPtCut),
  93 fMind0rphiCut(source.fMind0rphiCut)
  94 {
  95   //
  96   // Copy constructor
  97   //
  98   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
  99   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
 100   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
 101   for(Int_t i=0; i<13; i++)  fDsCuts[i]=source.fDsCuts[i];
 102   for(Int_t i=0; i<12; i++)  fLcCuts[i]=source.fLcCuts[i];
 103 }
 104 //--------------------------------------------------------------------------
 105 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
 106 {
 107   //
 108   // assignment operator
 109   //
 110   if(&source == this) return *this;
 111   fInputAOD = source.fInputAOD;
 112   fBzkG = source.fBzkG;
 113   fSecVtxWithKF = source.fSecVtxWithKF;
 114   fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
 115   fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
 116   fV1 = source.fV1;
 117   fD0toKpi = source.fD0toKpi;
 118   fJPSItoEle = source.fJPSItoEle;
 119   f3Prong = source.f3Prong;
 120   f4Prong = source.f4Prong;
 121   fITSrefit = source.fITSrefit;
 122   fBothSPD = source.fBothSPD;
 123   fMinITSCls = source.fMinITSCls;
 124   fMinPtCut = source.fMinPtCut;
 125   fMind0rphiCut = source.fMind0rphiCut;
 126
 127   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
 128   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
 129   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
 130   for(Int_t i=0; i<13; i++)  fDsCuts[i]=source.fDsCuts[i];
 131   for(Int_t i=0; i<12; i++)  fLcCuts[i]=source.fLcCuts[i];
 132
 133   return *this;
 134 }
 135 //----------------------------------------------------------------------------
 136 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
 137   // Destructor
 138   if(fV1) { delete fV1; fV1=0; }
 139 }
 140 //----------------------------------------------------------------------------
 141 void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
 142                                             TClonesArray *aodVerticesHFTClArr,
 143                                             TClonesArray *aodD0toKpiTClArr,
 144                                             TClonesArray *aodJPSItoEleTClArr,
 145                                             TClonesArray *aodCharm3ProngTClArr,
 146                                             TClonesArray *aodCharm4ProngTClArr)
 147 {
 148   // Find heavy-flavour vertex candidates
 149   // Input:  ESD or AOD
 150   // Output: AOD (additional branches added)
 151
 152
 153   TString evtype = event->IsA()->GetName();
 154   fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
 155
 156   if(!aodVerticesHFTClArr) {
 157     printf("ERROR: no aodVerticesHFTClArr");
 158     return;
 159   }
 160   if(fD0toKpi && !aodD0toKpiTClArr) {
 161     printf("ERROR: no aodD0toKpiTClArr");
 162     return;
 163   }
 164   if(fJPSItoEle && !aodJPSItoEleTClArr) {
 165     printf("ERROR: no aodJPSItoEleTClArr");
 166     return;
 167   }
 168   if(f3Prong && !aodCharm3ProngTClArr) {
 169     printf("ERROR: no aodCharm3ProngTClArr");
 170     return;
 171   }
 172   if(f4Prong && !aodCharm4ProngTClArr) {
 173     printf("ERROR: no aodCharm4ProngTClArr");
 174     return;
 175   }
 176
 177   // delete candidates from previous event and create references
 178   Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0;
 179   aodVerticesHFTClArr->Delete();
 180   iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
 181   TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
 182   if(fD0toKpi)   {
 183     aodD0toKpiTClArr->Delete();
 184     iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
 185   }
 186   if(fJPSItoEle) {
 187     aodJPSItoEleTClArr->Delete();
 188     iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
 189   }
 190   if(f3Prong) {
 191     aodCharm3ProngTClArr->Delete();
 192     i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
 193   }
 194   if(f4Prong) {
 195     aodCharm4ProngTClArr->Delete();
 196     i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
 197   }
 198   TClonesArray &aodD0toKpiRef     = *aodD0toKpiTClArr;
 199   TClonesArray &aodJPSItoEleRef   = *aodJPSItoEleTClArr;
 200   TClonesArray &aodCharm3ProngRef = *aodCharm3ProngTClArr;
 201   TClonesArray &aodCharm4ProngRef = *aodCharm4ProngTClArr;
 202
 203
 204   AliAODRecoDecayHF2Prong *io2Prong = 0;
 205   AliAODRecoDecayHF3Prong *io3Prong = 0;
 206   AliAODRecoDecayHF4Prong *io4Prong = 0;
 207
 208   Int_t    iTrkP1,iTrkP2,iTrkN1,iTrkN2,trkEntries;
 209   Int_t    nTrksP=0,nTrksN=0;
 210   Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2;
 211   Bool_t   okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
 212   AliESDtrack *postrack1 = 0;
 213   AliESDtrack *postrack2 = 0;
 214   AliESDtrack *negtrack1 = 0;
 215   AliESDtrack *negtrack2 = 0;
 216   Double_t dcaMax = fD0toKpiCuts[1];
 217   if(dcaMax < fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1];
 218   if(dcaMax < fDplusCuts[11])  dcaMax=fDplusCuts[11];
 219   AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
 220
 221
 222   // get Bz
 223   fBzkG = (Double_t)event->GetMagneticField();
 224
 225   trkEntries = (Int_t)event->GetNumberOfTracks();
 226   AliDebug(1,Form(" Number of tracks: %d",trkEntries));
 227
 228   if(trkEntries<2 || !event->GetPrimaryVertex()) {
 229     return;
 230   }
 231
 232   // call function that applies sigle-track selection,
 233   // separates positives and negatives,
 234   // and retrieves primary vertex
 235   TObjArray trksP(trkEntries/2);
 236   TObjArray trksN(trkEntries/2);
 237   SelectTracksAndCopyVertex(event,trksP,nTrksP,trksN,nTrksN);
 238
 239   AliDebug(1,Form(" Pos. tracks: %d    Neg. tracks: %d",nTrksP,nTrksN));
 240
 241   TObjArray *twoTrackArray1  = new TObjArray(2);
 242   TObjArray *twoTrackArray2  = new TObjArray(2);
 243   TObjArray *threeTrackArray = new TObjArray(3);
 244   TObjArray *fourTrackArray  = new TObjArray(4);
 245
 246   Double_t dispersion;
 247
 248   // LOOP ON  POSITIVE  TRACKS
 249   for(iTrkP1=0; iTrkP1<nTrksP; iTrkP1++) {
 250     if(iTrkP1%1==0) AliDebug(1,Form("  Processing positive track number %d of %d",iTrkP1,nTrksP));
 251     // get track from tracks array
 252     postrack1 = (AliESDtrack*)trksP.UncheckedAt(iTrkP1);
 253
 254     // LOOP ON  NEGATIVE  TRACKS
 255     for(iTrkN1=0; iTrkN1<nTrksN; iTrkN1++) {
 256       if(iTrkN1%1==0) AliDebug(1,Form("    Processing negative track number %d of %d",iTrkN1,nTrksN));
 257       // get track from tracks array
 258       negtrack1 = (AliESDtrack*)trksN.UncheckedAt(iTrkN1);
 259       // back to primary vertex
 260       postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 261       negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 262       // DCA between the two tracks
 263       dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 264       if(dcap1n1>dcaMax) { negtrack1=0; continue; }
 265       // Vertexing
 266       twoTrackArray1->AddAt(postrack1,0);
 267       twoTrackArray1->AddAt(negtrack1,1);
 268       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
 269       if(!vertexp1n1) {
 270         twoTrackArray1->Clear();
 271         negtrack1=0;
 272         continue;
 273       }
 274       if(fD0toKpi || fJPSItoEle) {
 275         io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI);
 276         if(okD0 || okJPSI) {
 277           AliAODVertex *v = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
 278           if(fInputAOD) AddDaughterRefs(v,event,twoTrackArray1);
 279           if(okD0) {
 280             AliAODRecoDecayHF2Prong *rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
 281             rd->SetSecondaryVtx(v);
 282             v->SetParent(rd);
 283           }
 284           if(okJPSI) {
 285             AliAODRecoDecayHF2Prong *rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
 286             rd->SetSecondaryVtx(v);
 287             if(!okD0) v->SetParent(rd); // do something better here...
 288           }
 289         }
 290         io2Prong=NULL;
 291       }
 292
 293       twoTrackArray1->Clear();
 294       if(!f3Prong && !f4Prong)  {
 295         negtrack1=0;
 296         delete vertexp1n1;
 297         continue;
 298       }
 299
 300
 301       // 2nd LOOP  ON  POSITIVE  TRACKS
 302       for(iTrkP2=iTrkP1+1; iTrkP2<nTrksP; iTrkP2++) {
 303         // get track from tracks array
 304         postrack2 = (AliESDtrack*)trksP.UncheckedAt(iTrkP2);
 305         // back to primary vertex
 306         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 307         postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 308         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 309         dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 310         if(dcap2n1>dcaMax) { postrack2=0; continue; }
 311         dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
 312         if(dcap1p2>dcaMax) { postrack2=0; continue; }
 313
 314         // Vertexing
 315         twoTrackArray2->AddAt(postrack2,0);
 316         twoTrackArray2->AddAt(negtrack1,1);
 317         AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
 318         if(!vertexp2n1) {
 319           twoTrackArray2->Clear();
 320           postrack2=0;
 321           continue;
 322         }
 323         if(f3Prong) {
 324           threeTrackArray->AddAt(postrack1,0);
 325           threeTrackArray->AddAt(negtrack1,1);
 326           threeTrackArray->AddAt(postrack2,2);
 327           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
 328           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
 329           if(ok3Prong) {
 330             AliAODVertex *v = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
 331             if(fInputAOD) AddDaughterRefs(v,event,threeTrackArray);
 332             AliAODRecoDecayHF3Prong *rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
 333             rd->SetSecondaryVtx(v);
 334             v->SetParent(rd);
 335           }
 336           if(io3Prong) io3Prong=NULL;
 337         }
 338         if(f4Prong) {
 339           // 3rd LOOP  ON  NEGATIVE  TRACKS (for 4 prong)
 340           for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) {
 341             // get track from tracks array
 342             negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2);
 343             // back to primary vertex
 344             postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 345             postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 346             negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 347             negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 348             dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 349             if(dcap1n2>dcaMax) { negtrack2=0; continue; }
 350             // Vertexing
 351             fourTrackArray->AddAt(postrack1,0);
 352             fourTrackArray->AddAt(negtrack1,1);
 353             fourTrackArray->AddAt(postrack2,2);
 354             fourTrackArray->AddAt(negtrack2,3);
 355             AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
 356             io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp2n1,dcap1n1,dcap1n2,dcap2n1,ok4Prong);
 357             if(ok4Prong) {
 358               AliAODVertex *v = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
 359               if(fInputAOD) AddDaughterRefs(v,event,fourTrackArray);
 360               AliAODRecoDecayHF4Prong *rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
 361               rd->SetSecondaryVtx(v);
 362               v->SetParent(rd);
 363             }
 364             if(io4Prong) io4Prong=NULL;
 365             fourTrackArray->Clear();
 366             negtrack2 = 0;
 367           } // end loop on negative tracks
 368         }
 369         postrack2 = 0;
 370         delete vertexp2n1;
 371       } // end 2nd loop on positive tracks
 372       twoTrackArray2->Clear();
 373
 374       // 2nd LOOP  ON  NEGATIVE  TRACKS
 375       for(iTrkN2=iTrkN1+1; iTrkN2<nTrksN; iTrkN2++) {
 376         // get track from tracks array
 377         negtrack2 = (AliESDtrack*)trksN.UncheckedAt(iTrkN2);
 378         // back to primary vertex
 379         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 380         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 381         negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 382         dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 383         if(dcap1n2>dcaMax) { negtrack2=0; continue; }
 384         dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 385         if(dcan1n2>dcaMax) { negtrack2=0; continue; }
 386
 387         // Vertexing
 388         twoTrackArray2->AddAt(postrack1,0);
 389         twoTrackArray2->AddAt(negtrack2,1);
 390
 391         AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
 392         if(!vertexp1n2) {
 393           twoTrackArray2->Clear();
 394           negtrack2=0;
 395           continue;
 396         }
 397         if(f3Prong) {
 398           threeTrackArray->AddAt(negtrack1,0);
 399           threeTrackArray->AddAt(postrack1,1);
 400           threeTrackArray->AddAt(negtrack2,2);
 401           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
 402           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
 403           if(ok3Prong) {
 404             AliAODVertex *v = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
 405             if(fInputAOD) AddDaughterRefs(v,event,threeTrackArray);
 406             AliAODRecoDecayHF3Prong *rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
 407             rd->SetSecondaryVtx(v);
 408             v->SetParent(rd);
 409           }
 410           if(io3Prong) io3Prong=NULL;
 411         }
 412         negtrack2 = 0;
 413         delete vertexp1n2;
 414       } // end 2nd loop on negative tracks
 415       twoTrackArray2->Clear();
 416
 417       negtrack1 = 0;
 418       delete vertexp1n1;
 419     } // end 1st loop on negative tracks
 420
 421     postrack1 = 0;
 422   }  // end 1st loop on positive tracks
 423
 424
 425   if(fD0toKpi) {
 426     AliDebug(1,Form(" D0->Kpi in event = %d;",
 427                     (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
 428   }
 429   if(fJPSItoEle) {
 430     AliDebug(1,Form(" JPSI->ee in event = %d;",
 431                     (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
 432   }
 433   if(f3Prong) {
 434     AliDebug(1,Form(" Charm->3Prong in event = %d;",
 435                     (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
 436   }
 437   if(f4Prong) {
 438     AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
 439                     (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
 440   }
 441
 442
 443   twoTrackArray1->Delete();  delete twoTrackArray1;
 444   twoTrackArray2->Delete();  delete twoTrackArray2;
 445   threeTrackArray->Clear();
 446   threeTrackArray->Delete(); delete threeTrackArray;
 447   fourTrackArray->Delete();  delete fourTrackArray;
 448
 449   if(fInputAOD) {
 450     trksP.Delete();
 451     trksN.Delete();
 452   }
 453
 454   return;
 455 }
 456 //----------------------------------------------------------------------------
 457 void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,AliVEvent *event,
 458                                              TObjArray *trkArray) const
 459 {
 460   // Add the AOD tracks as daughters of the vertex (TRef)
 461
 462   Int_t nTrks = trkArray->GetEntriesFast();
 463
 464   AliESDtrack *esdTrack = 0;
 465   AliAODTrack *aodTrack = 0;
 466   Int_t id;
 467
 468   for(Int_t i=0; i<nTrks; i++) {
 469     esdTrack = (AliESDtrack*)trkArray->UncheckedAt(i);
 470     id = (Int_t)esdTrack->GetID();
 471     aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
 472     v->AddDaughter(aodTrack);
 473   }
 474
 475   return;
 476 }
 477 //----------------------------------------------------------------------------
 478 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
 479                                    TObjArray *twoTrackArray1,AliVEvent *event,
 480                                    AliAODVertex *secVert,Double_t dca,
 481                                    Bool_t &okD0,Bool_t &okJPSI) const
 482 {
 483   // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
 484   // reconstruction cuts
 485   // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
 486
 487   okD0=kFALSE; okJPSI=kFALSE;
 488
 489   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
 490
 491   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(0);
 492   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray1->UncheckedAt(1);
 493
 494   // propagate tracks to secondary vertex, to compute inv. mass
 495   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
 496   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
 497
 498   Double_t momentum[3];
 499   postrack->GetPxPyPz(momentum);
 500   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 501   negtrack->GetPxPyPz(momentum);
 502   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 503
 504
 505   // invariant mass cut (try to improve coding here..)
 506   Bool_t okMassCut=kFALSE;
 507   if(!okMassCut && fD0toKpi)   if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
 508   if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
 509   if(!okMassCut) {
 510     AliDebug(2," candidate didn't pass mass cut");
 511     return 0x0;
 512   }
 513
 514   // primary vertex to be used by this candidate
 515   AliAODVertex *primVertexAOD  = OwnPrimaryVertex(twoTrackArray1,event);
 516   if(!primVertexAOD) return 0x0;
 517
 518   Double_t d0z0[2],covd0z0[3];
 519   postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 520   d0[0] = d0z0[0];
 521   d0err[0] = TMath::Sqrt(covd0z0[0]);
 522   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 523   d0[1] = d0z0[0];
 524   d0err[1] = TMath::Sqrt(covd0z0[0]);
 525
 526   // create the object AliAODRecoDecayHF2Prong
 527   AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
 528   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
 529   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
 530   the2Prong->SetProngIDs(2,id);
 531
 532
 533   // select D0->Kpi
 534   Int_t checkD0,checkD0bar;
 535   if(fD0toKpi) okD0 = the2Prong->SelectD0(fD0toKpiCuts,checkD0,checkD0bar);
 536   //if(fDebug && fD0toKpi) printf("   %d\n",(Int_t)okD0);
 537   // select J/psi from B
 538   Int_t checkJPSI;
 539   if(fJPSItoEle) okJPSI = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI);
 540   //if(fDebug && fJPSItoEle) printf("   %d\n",(Int_t)okJPSI);
 541
 542   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) the2Prong->UnsetOwnPrimaryVtx();
 543
 544   if(okD0 || okJPSI) {
 545     // get PID info from ESD
 546     Double_t esdpid0[5];
 547     postrack->GetESDpid(esdpid0);
 548     Double_t esdpid1[5];
 549     negtrack->GetESDpid(esdpid1);
 550     Double_t esdpid[10];
 551     for(Int_t i=0;i<5;i++) {
 552       esdpid[i]   = esdpid0[i];
 553       esdpid[5+i] = esdpid1[i];
 554     }
 555     the2Prong->SetPID(2,esdpid);
 556   }
 557
 558   return the2Prong;
 559 }
 560 //----------------------------------------------------------------------------
 561 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
 562                              TObjArray *threeTrackArray,AliVEvent *event,
 563                              AliAODVertex *secVert,Double_t dispersion,
 564                              AliAODVertex *vertexp1n1,AliAODVertex *vertexp2n1,
 565                              Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
 566                              Bool_t &ok3Prong) const
 567 {
 568   // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
 569   // reconstruction cuts
 570   // E.Bruna, F.Prino
 571
 572
 573   ok3Prong=kFALSE;
 574   if(!secVert) return 0x0;
 575
 576   Double_t px[3],py[3],pz[3],d0[3],d0err[3];
 577   Double_t momentum[3];
 578
 579
 580   AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
 581   AliESDtrack *negtrack  = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
 582   AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
 583
 584   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
 585   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
 586   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
 587   postrack1->GetPxPyPz(momentum);
 588   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 589   negtrack->GetPxPyPz(momentum);
 590   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 591   postrack2->GetPxPyPz(momentum);
 592   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
 593
 594   // invariant mass cut for D+, Ds, Lc
 595   Bool_t okMassCut=kFALSE;
 596   if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
 597   if(!okMassCut) {
 598     AliDebug(2," candidate didn't pass mass cut");
 599     return 0x0;
 600   }
 601
 602   // primary vertex to be used by this candidate
 603   AliAODVertex *primVertexAOD  = OwnPrimaryVertex(threeTrackArray,event);
 604   if(!primVertexAOD) return 0x0;
 605
 606   Double_t d0z0[2],covd0z0[3];
 607   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 608   d0[0]=d0z0[0];
 609   d0err[0] = TMath::Sqrt(covd0z0[0]);
 610   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 611   d0[1]=d0z0[0];
 612   d0err[1] = TMath::Sqrt(covd0z0[0]);
 613   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 614   d0[2]=d0z0[0];
 615   d0err[2] = TMath::Sqrt(covd0z0[0]);
 616
 617
 618   // create the object AliAODRecoDecayHF3Prong
 619   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
 620   Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
 621   Double_t dist12=TMath::Sqrt((vertexp1n1->GetX()-pos[0])*(vertexp1n1->GetX()-pos[0])+(vertexp1n1->GetY()-pos[1])*(vertexp1n1->GetY()-pos[1])+(vertexp1n1->GetZ()-pos[2])*(vertexp1n1->GetZ()-pos[2]));
 622   Double_t dist23=TMath::Sqrt((vertexp2n1->GetX()-pos[0])*(vertexp2n1->GetX()-pos[0])+(vertexp2n1->GetY()-pos[1])*(vertexp2n1->GetY()-pos[1])+(vertexp2n1->GetZ()-pos[2])*(vertexp2n1->GetZ()-pos[2]));
 623   Short_t charge=(Short_t)(postrack1->Charge()*postrack2->Charge()*negtrack->Charge());
 624   AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
 625   the3Prong->SetOwnPrimaryVtx(primVertexAOD);
 626   UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
 627   the3Prong->SetProngIDs(3,id);
 628
 629
 630   // select D+->Kpipi, Ds->KKpi, Lc->pKpi
 631   if(f3Prong) {
 632     ok3Prong = kFALSE;
 633     Int_t ok1,ok2;
 634     if(the3Prong->SelectDplus(fDplusCuts))   ok3Prong = kTRUE;
 635     if(the3Prong->SelectDs(fDsCuts,ok1,ok2)) ok3Prong = kTRUE;
 636     if(the3Prong->SelectLc(fLcCuts,ok1,ok2)) ok3Prong = kTRUE;
 637   }
 638   //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
 639
 640   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) the3Prong->UnsetOwnPrimaryVtx();
 641
 642   if(ok3Prong) {
 643     // get PID info from ESD
 644     Double_t esdpid0[5];
 645     postrack1->GetESDpid(esdpid0);
 646     Double_t esdpid1[5];
 647     negtrack->GetESDpid(esdpid1);
 648     Double_t esdpid2[5];
 649     postrack2->GetESDpid(esdpid2);
 650
 651     Double_t esdpid[15];
 652     for(Int_t i=0;i<5;i++) {
 653       esdpid[i]   = esdpid0[i];
 654       esdpid[5+i] = esdpid1[i];
 655       esdpid[10+i] = esdpid2[i];
 656     }
 657     the3Prong->SetPID(3,esdpid);
 658   }
 659
 660   return the3Prong;
 661 }
 662 //----------------------------------------------------------------------------
 663 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
 664                              TObjArray *fourTrackArray,AliVEvent *event,
 665                              AliAODVertex *secVert,
 666                              AliAODVertex *vertexp1n1,AliAODVertex *vertexp2n1,
 667                              Double_t dcap1n1,Double_t dcap1n2,Double_t dcap2n1,
 668                              Bool_t &ok4Prong) const
 669 {
 670   // Make 4Prong candidates and check if they pass D0toKpipipi
 671   // reconstruction cuts
 672   // G.E.Bruno, R.Romita
 673
 674   ok4Prong=kFALSE;
 675   if(!secVert) return 0x0;
 676
 677   Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
 678
 679   px[0]=dcap1n1*dcap1n2*dcap2n1; // TO BE CHANGED (done just to removed compilation warning about dca... not used)
 680
 681   AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
 682   AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
 683   AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
 684   AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
 685
 686   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
 687   negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
 688   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
 689   negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
 690
 691   Double_t momentum[3];
 692   postrack1->GetPxPyPz(momentum);
 693   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
 694   negtrack1->GetPxPyPz(momentum);
 695   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
 696   postrack2->GetPxPyPz(momentum);
 697   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
 698   negtrack2->GetPxPyPz(momentum);
 699   px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
 700
 701   // invariant mass cut for rho or D0 (try to improve coding here..)
 702   //Bool_t okMassCut=kFALSE;
 703   //if(!okMassCut) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
 704   //if(!okMassCut) {
 705   //  if(fDebug) printf(" candidate didn't pass mass cut\n");
 706   //  return 0x0;
 707   //}
 708
 709   // primary vertex to be used by this candidate
 710   AliAODVertex *primVertexAOD  = OwnPrimaryVertex(fourTrackArray,event);
 711   if(!primVertexAOD) return 0x0;
 712
 713   /*
 714     Double_t d0z0[2],covd0z0[3];
 715     postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 716     negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 717     postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 718     negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
 719   */
 720
 721   // create the object AliAODRecoDecayHF4Prong
 722   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
 723   Double_t dca[6]={0.,0.,0.,0.,0.,0.}; //  modify it
 724   Double_t dist12=TMath::Sqrt((vertexp1n1->GetX()-pos[0])*(vertexp1n1->GetX()-pos[0])+(vertexp1n1->GetY()-pos[1])*(vertexp1n1->GetY()-pos[1])+(vertexp1n1->GetZ()-pos[2])*(vertexp1n1->GetZ()-pos[2]));
 725   Double_t dist23=TMath::Sqrt((vertexp2n1->GetX()-pos[0])*(vertexp2n1->GetX()-pos[0])+(vertexp2n1->GetY()-pos[1])*(vertexp2n1->GetY()-pos[1])+(vertexp2n1->GetZ()-pos[2])*(vertexp2n1->GetZ()-pos[2]));
 726   Double_t dist14=0.; // to be implemented
 727   Double_t dist34=0.; // to be implemented
 728   Short_t charge=0;
 729   AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist23,dist14,dist34,charge);
 730   the4Prong->SetOwnPrimaryVtx(primVertexAOD);
 731   UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
 732   the4Prong->SetProngIDs(4,id);
 733
 734
 735   // use the following two lines once AliAODRecoDecayHF4Prong::SelectD0 is available
 736   // select D0->Kpipipi
 737   //Int_t checkD0,checkD0bar;
 738   // ok4Prong=the4Prong->SelectD0(fD04pCuts,checkD0,checkD0bar);
 739   ok4Prong=kFALSE;  //for the time being ...
 740
 741   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) the4Prong->UnsetOwnPrimaryVtx();
 742
 743   // get PID info from ESD
 744   Double_t esdpid0[5];
 745   postrack1->GetESDpid(esdpid0);
 746   Double_t esdpid1[5];
 747   negtrack1->GetESDpid(esdpid1);
 748   Double_t esdpid2[5];
 749   postrack2->GetESDpid(esdpid2);
 750   Double_t esdpid3[5];
 751   negtrack2->GetESDpid(esdpid3);
 752
 753   Double_t esdpid[20];
 754   for(Int_t i=0;i<5;i++) {
 755     esdpid[i]   = esdpid0[i];
 756     esdpid[5+i] = esdpid1[i];
 757     esdpid[10+i] = esdpid2[i];
 758     esdpid[15+i] = esdpid3[i];
 759   }
 760   the4Prong->SetPID(4,esdpid);
 761
 762   return the4Prong;
 763 }
 764 //-----------------------------------------------------------------------------
 765 AliAODVertex* AliAnalysisVertexingHF::OwnPrimaryVertex(TObjArray *trkArray,
 766                                                        AliVEvent *event) const
 767 {
 768   // Returns primary vertex specific to this candidate
 769
 770   AliESDVertex *vertexESD = 0;
 771   AliAODVertex *vertexAOD = 0;
 772
 773   Int_t nTrks = trkArray->GetEntriesFast();
 774
 775   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
 776     // primary vertex from the input event
 777
 778     vertexESD = new AliESDVertex(*fV1);
 779
 780   } else {
 781     // primary vertex specific to this candidate
 782
 783     AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
 784
 785     if(fRecoPrimVtxSkippingTrks) {
 786       // recalculating the vertex
 787
 788       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
 789         Float_t diamondcovxy[3];
 790         event->GetDiamondCovXY(diamondcovxy);
 791         Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
 792         Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.};
 793         AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
 794         vertexer->SetVtxStart(diamond);
 795         delete diamond; diamond=NULL;
 796         if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
 797           vertexer->SetOnlyFitter();
 798       }
 799       Int_t skipped[10];
 800       AliESDtrack *t = 0;
 801       for(Int_t i=0; i<nTrks; i++) {
 802         t = (AliESDtrack*)trkArray->UncheckedAt(i);
 803         skipped[i] = (Int_t)t->GetID();
 804       }
 805       vertexer->SetSkipTracks(nTrks,skipped);
 806       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
 807
 808     } else if(fRmTrksFromPrimVtx) {
 809       // removing the prongs tracks
 810
 811       TObjArray rmArray(nTrks);
 812       UShort_t *rmId = new UShort_t[nTrks];
 813       AliESDtrack *esdTrack = 0;
 814       AliESDtrack *t = 0;
 815       for(Int_t i=0; i<nTrks; i++) {
 816         t = (AliESDtrack*)trkArray->UncheckedAt(i);
 817         esdTrack = new AliESDtrack(*t);
 818         rmArray.AddLast(esdTrack);
 819         rmId[i]=(UShort_t)esdTrack->GetID();
 820       }
 821       Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
 822       vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
 823       delete [] rmId; rmId=NULL;
 824       rmArray.Delete();
 825
 826     }
 827
 828     if(!vertexESD) return vertexAOD;
 829     if(vertexESD->GetNContributors()<=0) {
 830       AliDebug(2,"vertexing failed");
 831       delete vertexESD; vertexESD=NULL;
 832       return vertexAOD;
 833     }
 834
 835     delete vertexer; vertexer=NULL;
 836
 837   }
 838
 839   // convert to AliAODVertex
 840   Double_t pos[3],cov[6],chi2perNDF;
 841   vertexESD->GetXYZ(pos); // position
 842   vertexESD->GetCovMatrix(cov); //covariance matrix
 843   chi2perNDF = vertexESD->GetChi2toNDF();
 844   delete vertexESD; vertexESD=NULL;
 845
 846   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
 847
 848   return vertexAOD;
 849 }
 850 //-----------------------------------------------------------------------------
 851 void AliAnalysisVertexingHF::PrintStatus() const {
 852   // Print parameters being used
 853
 854   printf("Preselections:\n");
 855   printf("    fITSrefit   = %d\n",(Int_t)fITSrefit);
 856   printf("    fBothSPD   = %d\n",(Int_t)fBothSPD);
 857   printf("    fMinITSCls   = %d\n",fMinITSCls);
 858   printf("    fMinPtCut   = %f GeV/c\n",fMinPtCut);
 859   printf("    fMind0rphiCut   = %f cm\n",fMind0rphiCut);
 860   if(fSecVtxWithKF) {
 861     printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
 862   } else {
 863     printf("Secondary vertex with AliVertexerTracks\n");
 864   }
 865   if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
 866   if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
 867   if(fD0toKpi) {
 868     printf("Reconstruct D0->Kpi candidates with cuts:\n");
 869     printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
 870     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
 871     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
 872     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
 873     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
 874     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
 875     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
 876     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
 877     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
 878   }
 879   if(fJPSItoEle) {
 880     printf("Reconstruct J/psi from B candidates with cuts:\n");
 881     printf("    |M-MJPSI| [GeV]    < %f\n",fBtoJPSICuts[0]);
 882     printf("    dca    [cm]  < %f\n",fBtoJPSICuts[1]);
 883     printf("    cosThetaStar     < %f\n",fBtoJPSICuts[2]);
 884     printf("    pTP     [GeV/c]    > %f\n",fBtoJPSICuts[3]);
 885     printf("    pTN    [GeV/c]    > %f\n",fBtoJPSICuts[4]);
 886     printf("    |d0P|  [cm]  < %f\n",fBtoJPSICuts[5]);
 887     printf("    |d0N| [cm]  < %f\n",fBtoJPSICuts[6]);
 888     printf("    d0d0  [cm^2] < %f\n",fBtoJPSICuts[7]);
 889     printf("    cosThetaPoint    > %f\n",fBtoJPSICuts[8]);
 890   }
 891   if(f3Prong) {
 892     printf("Reconstruct 3 prong candidates.\n");
 893     printf("  D+->Kpipi cuts:\n");
 894     printf("    |M-MD+| [GeV]    < %f\n",fDplusCuts[0]);
 895     printf("    pTK     [GeV/c]    > %f\n",fDplusCuts[1]);
 896     printf("    pTPi    [GeV/c]    > %f\n",fDplusCuts[2]);
 897     printf("    |d0K|  [cm]  > %f\n",fDplusCuts[3]);
 898     printf("    |d0Pi| [cm]  > %f\n",fDplusCuts[4]);
 899     printf("    dist12    [cm]  < %f\n",fDplusCuts[5]);
 900     printf("    sigmavert [cm]   < %f\n",fDplusCuts[6]);
 901     printf("    dist prim-sec [cm] > %f\n",fDplusCuts[7]);
 902     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDplusCuts[8]);
 903     printf("    cosThetaPoint    > %f\n",fDplusCuts[9]);
 904     printf("    Sum d0^2 [cm^2]  > %f\n",fDplusCuts[10]);
 905     printf("    dca cut [cm]  < %f\n",fDplusCuts[11]);
 906     printf("  Ds->KKpi cuts:\n");
 907     printf("    |M-MDs| [GeV]    < %f\n",fDsCuts[0]);
 908     printf("    pTK     [GeV/c]    > %f\n",fDsCuts[1]);
 909     printf("    pTPi    [GeV/c]    > %f\n",fDsCuts[2]);
 910     printf("    |d0K|  [cm]  > %f\n",fDsCuts[3]);
 911     printf("    |d0Pi| [cm]  > %f\n",fDsCuts[4]);
 912     printf("    dist12    [cm]  < %f\n",fDsCuts[5]);
 913     printf("    sigmavert [cm]   < %f\n",fDsCuts[6]);
 914     printf("    dist prim-sec [cm] > %f\n",fDsCuts[7]);
 915     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDsCuts[8]);
 916     printf("    cosThetaPoint    > %f\n",fDsCuts[9]);
 917     printf("    Sum d0^2 [cm^2]  > %f\n",fDsCuts[10]);
 918     printf("    dca cut [cm]  < %f\n",fDsCuts[11]);
 919     printf("    Inv. Mass  phi/K0* [GeV]  < %f\n",fDsCuts[12]);
 920     printf("  Lc->pKpi cuts:\n");
 921     printf("    |M-MLc| [GeV]    < %f\n",fLcCuts[0]);
 922     printf("    pTP     [GeV/c]    > %f\n",fLcCuts[1]);
 923     printf("    pTPi and pTK [GeV/c]    > %f\n",fLcCuts[2]);
 924     printf("    |d0P|  [cm]  > %f\n",fLcCuts[3]);
 925     printf("    |d0Pi| and |d0K| [cm]  > %f\n",fLcCuts[4]);
 926     printf("    dist12    [cm]  < %f\n",fLcCuts[5]);
 927     printf("    sigmavert [cm]   < %f\n",fLcCuts[6]);
 928     printf("    dist prim-sec [cm] > %f\n",fLcCuts[7]);
 929     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fLcCuts[8]);
 930     printf("    cosThetaPoint    > %f\n",fLcCuts[9]);
 931     printf("    Sum d0^2 [cm^2]  > %f\n",fLcCuts[10]);
 932     printf("    dca cut [cm]  < %f\n",fLcCuts[11]);
 933     printf("  Ds->KKpi cuts:\n");
 934   }
 935
 936   return;
 937 }
 938 //-----------------------------------------------------------------------------
 939 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
 940                                                                  Double_t &dispersion) const
 941 {
 942   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
 943
 944   AliESDVertex *vertexESD = 0;
 945   AliAODVertex *vertexAOD = 0;
 946
 947   if(!fSecVtxWithKF) { // AliVertexerTracks
 948
 949     AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
 950     vertexer->SetVtxStart(fV1);
 951     vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
 952     delete vertexer; vertexer=NULL;
 953
 954     if(!vertexESD) return vertexAOD;
 955
 956     if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
 957       AliDebug(2,"vertexing failed");
 958       delete vertexESD; vertexESD=NULL;
 959       return vertexAOD;
 960     }
 961
 962   } else { // Kalman Filter vertexer (AliKFParticle)
 963
 964     AliKFParticle::SetField(fBzkG);
 965
 966     AliKFParticle vertexKF;
 967
 968     Int_t nTrks = trkArray->GetEntriesFast();
 969     for(Int_t i=0; i<nTrks; i++) {
 970       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
 971       AliKFParticle daughterKF(*esdTrack,211);
 972       vertexKF.AddDaughter(daughterKF);
 973     }
 974     vertexESD = new AliESDVertex();
 975     vertexKF.CopyToESDVertex(*vertexESD);
 976
 977   }
 978
 979   // convert to AliAODVertex
 980   Double_t pos[3],cov[6],chi2perNDF;
 981   vertexESD->GetXYZ(pos); // position
 982   vertexESD->GetCovMatrix(cov); //covariance matrix
 983   chi2perNDF = vertexESD->GetChi2toNDF();
 984   dispersion = vertexESD->GetDispersion();
 985   delete vertexESD; vertexESD=NULL;
 986
 987   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
 988   vertexAOD->SetType(AliAODVertex::kUndef); // to be changed
 989
 990   return vertexAOD;
 991 }
 992 //-----------------------------------------------------------------------------
 993 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
 994                                              Int_t nprongs,
 995                                              Double_t *px,
 996                                              Double_t *py,
 997                                              Double_t *pz) const {
 998   // Check invariant mass cut
 999
1000   Short_t dummycharge=0;
1001   Double_t *dummyd0 = new Double_t[nprongs];
1002   for(Int_t ip=0;ip<nprongs;ip++) dummyd0[ip]=0.;
1003   AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
1004   delete [] dummyd0;
1005
1006   UInt_t pdg2[2],pdg3[3];
1007   Double_t mPDG,minv;
1008
1009   Bool_t retval=kFALSE;
1010   switch (decay)
1011     {
1012     case 0:                  // D0->Kpi
1013       pdg2[0]=211; pdg2[1]=321;
1014       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1015       minv = rd->InvMass(nprongs,pdg2);
1016       if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1017       pdg2[0]=321; pdg2[1]=211;
1018       minv = rd->InvMass(nprongs,pdg2);
1019       if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1020       break;
1021     case 1:                  // JPSI->ee
1022       pdg2[0]=11; pdg2[1]=11;
1023       mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
1024       minv = rd->InvMass(nprongs,pdg2);
1025       if(TMath::Abs(minv-mPDG)<fBtoJPSICuts[0]) retval=kTRUE;
1026       break;
1027     case 2:                  // D+->Kpipi
1028       pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
1029       mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
1030       minv = rd->InvMass(nprongs,pdg3);
1031       if(TMath::Abs(minv-mPDG)<fDplusCuts[0]) retval=kTRUE;
1032                             // Ds+->KKpi
1033       pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
1034       mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1035       minv = rd->InvMass(nprongs,pdg3);
1036       if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1037       pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
1038       minv = rd->InvMass(nprongs,pdg3);
1039       if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1040                             // Lc->pKpi
1041       pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
1042       mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1043       minv = rd->InvMass(nprongs,pdg3);
1044       if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1045       pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
1046       minv = rd->InvMass(nprongs,pdg3);
1047       if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1048       break;
1049     default:
1050       break;
1051     }
1052
1053   delete rd;
1054
1055   return retval;
1056 }
1057 //-----------------------------------------------------------------------------
1058 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(AliVEvent *event,
1059                                           TObjArray &trksP,Int_t &nTrksP,
1060                                           TObjArray &trksN,Int_t &nTrksN)
1061 {
1062   // Fill two TObjArrays with positive and negative tracks and
1063   // apply single-track preselection.
1064   // Create the AliESDVertex object (convert from AliAODVertex if necessary)
1065   // In case of AOD input, also fill fAODMap for track index<->ID
1066
1067   const AliVVertex *vprimary = event->GetPrimaryVertex();
1068
1069   if(fV1) { delete fV1; fV1=NULL; }
1070   if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
1071
1072   Int_t nindices=0;
1073   UShort_t *indices = 0;
1074   Double_t pos[3],cov[6];
1075
1076   if(!fInputAOD) { // ESD
1077     fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
1078   } else {         // AOD
1079     vprimary->GetXYZ(pos);
1080     vprimary->GetCovarianceMatrix(cov);
1081     fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1082     indices = new UShort_t[event->GetNumberOfTracks()];
1083     fAODMap = new Int_t[100000];
1084   }
1085
1086
1087     for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++) {
1088
1089     }
1090
1091
1092   nTrksP=0,nTrksN=0;
1093
1094   Int_t entries = (Int_t)event->GetNumberOfTracks();
1095
1096   // transfer ITS tracks from event to arrays
1097   for(Int_t i=0; i<entries; i++) {
1098
1099     AliVTrack *track = (AliVTrack*)event->GetTrack(i);
1100
1101     if(fInputAOD) {
1102       AliAODTrack *aodt = (AliAODTrack*)track;
1103       if(aodt->GetUsedForPrimVtxFit()) {
1104         indices[nindices]=aodt->GetID(); nindices++;
1105       }
1106       fAODMap[(Int_t)aodt->GetID()] = i;
1107     }
1108
1109     // require refit in ITS
1110     UInt_t status = track->GetStatus();
1111     if(fITSrefit && !(status&AliESDtrack::kITSrefit)) {
1112       AliDebug(2,Form("track %d is not kITSrefit",i));
1113       continue;
1114     }
1115
1116     Int_t ncls0=0;
1117     for(Int_t l=0;l<6;l++) if(TESTBIT(track->GetITSClusterMap(),l)) ncls0++;
1118
1119     // require minimum # of ITS points
1120     if(ncls0<fMinITSCls)  {
1121       AliDebug(2,Form("track %d has %d ITS cls",i,ncls0));
1122       continue;
1123     }
1124     // require points on the 2 pixel layers
1125     if(fBothSPD)
1126       if(!TESTBIT(track->GetITSClusterMap(),0) ||
1127          !TESTBIT(track->GetITSClusterMap(),1)) continue;
1128
1129     AliESDtrack *esdt = 0;
1130     if(!fInputAOD) {
1131       esdt = (AliESDtrack*)track;
1132     } else {
1133       esdt = new AliESDtrack(track);
1134     }
1135
1136     // single track selection
1137     if(!SingleTrkCuts(*esdt)) {
1138       if(fInputAOD) delete esdt;
1139       esdt = NULL;
1140       continue;
1141     }
1142
1143     if(track->Charge()<0) { // negative track
1144       trksN.AddLast(esdt);
1145       nTrksN++;
1146     } else {                 // positive track
1147       trksP.AddLast(esdt);
1148       nTrksP++;
1149     }
1150
1151   } // end loop on tracks
1152
1153   // primary vertex from AOD
1154   if(fInputAOD) {
1155     delete fV1; fV1=NULL;
1156     vprimary->GetXYZ(pos);
1157     vprimary->GetCovarianceMatrix(cov);
1158     Double_t chi2toNDF = vprimary->GetChi2perNDF();
1159     Int_t ncontr=nindices;
1160     if(vprimary->GetTitle()=="VertexerTracksWithContraint") ncontr += 1;
1161     Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
1162     fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
1163     fV1->SetTitle(vprimary->GetTitle());
1164     fV1->SetIndices(nindices,indices);
1165   }
1166   if(indices) { delete [] indices; indices=NULL; }
1167
1168
1169   return;
1170 }
1171 //-----------------------------------------------------------------------------
1172 void AliAnalysisVertexingHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
1173                                    Double_t cut2,Double_t cut3,Double_t cut4,
1174                                    Double_t cut5,Double_t cut6,
1175                                    Double_t cut7,Double_t cut8)
1176 {
1177   // Set the cuts for D0 selection
1178   fD0toKpiCuts[0] = cut0;
1179   fD0toKpiCuts[1] = cut1;
1180   fD0toKpiCuts[2] = cut2;
1181   fD0toKpiCuts[3] = cut3;
1182   fD0toKpiCuts[4] = cut4;
1183   fD0toKpiCuts[5] = cut5;
1184   fD0toKpiCuts[6] = cut6;
1185   fD0toKpiCuts[7] = cut7;
1186   fD0toKpiCuts[8] = cut8;
1187
1188   return;
1189 }
1190 //-----------------------------------------------------------------------------
1191 void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9])
1192 {
1193   // Set the cuts for D0 selection
1194
1195   for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
1196
1197   return;
1198 }
1199 //-----------------------------------------------------------------------------
1200 void AliAnalysisVertexingHF::SetBtoJPSICuts(Double_t cut0,Double_t cut1,
1201                                    Double_t cut2,Double_t cut3,Double_t cut4,
1202                                    Double_t cut5,Double_t cut6,
1203                                    Double_t cut7,Double_t cut8)
1204 {
1205   // Set the cuts for J/psi from B selection
1206   fBtoJPSICuts[0] = cut0;
1207   fBtoJPSICuts[1] = cut1;
1208   fBtoJPSICuts[2] = cut2;
1209   fBtoJPSICuts[3] = cut3;
1210   fBtoJPSICuts[4] = cut4;
1211   fBtoJPSICuts[5] = cut5;
1212   fBtoJPSICuts[6] = cut6;
1213   fBtoJPSICuts[7] = cut7;
1214   fBtoJPSICuts[8] = cut8;
1215
1216   return;
1217 }
1218 //-----------------------------------------------------------------------------
1219 void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9])
1220 {
1221   // Set the cuts for J/psi from B selection
1222
1223   for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i];
1224
1225   return;
1226 }
1227 //-----------------------------------------------------------------------------
1228 void AliAnalysisVertexingHF::SetDplusCuts(Double_t cut0,Double_t cut1,
1229                                    Double_t cut2,Double_t cut3,Double_t cut4,
1230                                    Double_t cut5,Double_t cut6,
1231                                    Double_t cut7,Double_t cut8,
1232                                    Double_t cut9,Double_t cut10,Double_t cut11)
1233 {
1234   // Set the cuts for Dplus->Kpipi selection
1235   fDplusCuts[0] = cut0;
1236   fDplusCuts[1] = cut1;
1237   fDplusCuts[2] = cut2;
1238   fDplusCuts[3] = cut3;
1239   fDplusCuts[4] = cut4;
1240   fDplusCuts[5] = cut5;
1241   fDplusCuts[6] = cut6;
1242   fDplusCuts[7] = cut7;
1243   fDplusCuts[8] = cut8;
1244   fDplusCuts[9] = cut9;
1245   fDplusCuts[10] = cut10;
1246   fDplusCuts[11] = cut11;
1247
1248   return;
1249 }
1250 //-----------------------------------------------------------------------------
1251 void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12])
1252 {
1253   // Set the cuts for Dplus->Kpipi selection
1254
1255   for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i];
1256
1257   return;
1258 }
1259 //-----------------------------------------------------------------------------
1260 void AliAnalysisVertexingHF::SetDsCuts(Double_t cut0,Double_t cut1,
1261                                    Double_t cut2,Double_t cut3,Double_t cut4,
1262                                    Double_t cut5,Double_t cut6,
1263                                    Double_t cut7,Double_t cut8,
1264                                    Double_t cut9,Double_t cut10,
1265                                    Double_t cut11,Double_t cut12)
1266 {
1267   // Set the cuts for Ds->KKpi selection
1268   fDsCuts[0] = cut0;
1269   fDsCuts[1] = cut1;
1270   fDsCuts[2] = cut2;
1271   fDsCuts[3] = cut3;
1272   fDsCuts[4] = cut4;
1273   fDsCuts[5] = cut5;
1274   fDsCuts[6] = cut6;
1275   fDsCuts[7] = cut7;
1276   fDsCuts[8] = cut8;
1277   fDsCuts[9] = cut9;
1278   fDsCuts[10] = cut10;
1279   fDsCuts[11] = cut11;
1280   fDsCuts[12] = cut12;
1281
1282   return;
1283 }
1284 //-----------------------------------------------------------------------------
1285 void AliAnalysisVertexingHF::SetDsCuts(const Double_t cuts[13])
1286 {
1287   // Set the cuts for Ds->KKpi selection
1288
1289   for(Int_t i=0; i<13; i++) fDsCuts[i] = cuts[i];
1290
1291   return;
1292 }
1293 //-----------------------------------------------------------------------------
1294 void AliAnalysisVertexingHF::SetLcCuts(Double_t cut0,Double_t cut1,
1295                                    Double_t cut2,Double_t cut3,Double_t cut4,
1296                                    Double_t cut5,Double_t cut6,
1297                                    Double_t cut7,Double_t cut8,
1298                                    Double_t cut9,Double_t cut10,Double_t cut11)
1299 {
1300   // Set the cuts for Lc->pKpi selection
1301   fLcCuts[0] = cut0;
1302   fLcCuts[1] = cut1;
1303   fLcCuts[2] = cut2;
1304   fLcCuts[3] = cut3;
1305   fLcCuts[4] = cut4;
1306   fLcCuts[5] = cut5;
1307   fLcCuts[6] = cut6;
1308   fLcCuts[7] = cut7;
1309   fLcCuts[8] = cut8;
1310   fLcCuts[9] = cut9;
1311   fLcCuts[10] = cut10;
1312   fLcCuts[11] = cut11;
1313
1314   return;
1315 }
1316 //-----------------------------------------------------------------------------
1317 void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[12])
1318 {
1319   // Set the cuts for Lc->pKpi selection
1320
1321   for(Int_t i=0; i<12; i++) fLcCuts[i] = cuts[i];
1322
1323   return;
1324 }
1325 //-----------------------------------------------------------------------------
1326 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack& trk) const
1327 {
1328   // Check if track passes some kinematical cuts
1329
1330   if(trk.Pt() < fMinPtCut) {
1331     //printf("pt %f\n",1./trk.GetParameter()[4]);
1332     return kFALSE;
1333   }
1334   Double_t d0z0[2],covd0z0[3];
1335   trk.PropagateToDCA(fV1,fBzkG,kVeryBig,d0z0,covd0z0);
1336   if(TMath::Abs(d0z0[0]) < fMind0rphiCut) {
1337     //printf("d0rphi %f\n",TMath::Abs(d0z0[0]));
1338     return kFALSE;
1339   }
1340
1341   return kTRUE;
1342 }
1343 //-----------------------------------------------------------------------------