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 //  Contact: andrea.dainese@pd.infn.it
  24 //  Contributors: E.Bruna, G.E.Bruno, A.Dainese, C.Di Gliglio,
  25 //                F.Prino, R.Romita, X.M.Zhang
  26 //----------------------------------------------------------------------------
  27 #include <Riostream.h>
  28 #include <TFile.h>
  29 #include <TDatabasePDG.h>
  30 #include <TString.h>
  31 #include <TList.h>
  32 #include "AliLog.h"
  33 #include "AliVEvent.h"
  34 #include "AliVVertex.h"
  35 #include "AliVTrack.h"
  36 #include "AliVertexerTracks.h"
  37 #include "AliKFVertex.h"
  38 #include "AliESDEvent.h"
  39 #include "AliESDVertex.h"
  40 #include "AliExternalTrackParam.h"
  41 #include "AliNeutralTrackParam.h"
  42 #include "AliESDtrack.h"
  43 #include "AliESDtrackCuts.h"
  44 #include "AliAODEvent.h"
  45 #include "AliAODRecoDecay.h"
  46 #include "AliAODRecoDecayHF.h"
  47 #include "AliAODRecoDecayHF2Prong.h"
  48 #include "AliAODRecoDecayHF3Prong.h"
  49 #include "AliAODRecoDecayHF4Prong.h"
  50 #include "AliAODRecoCascadeHF.h"
  51 #include "AliRDHFCutsD0toKpi.h"
  52 #include "AliRDHFCutsJpsitoee.h"
  53 #include "AliRDHFCutsDplustoKpipi.h"
  54 #include "AliRDHFCutsDstoKKpi.h"
  55 #include "AliRDHFCutsLctopKpi.h"
  56 #include "AliRDHFCutsLctoV0.h"
  57 #include "AliRDHFCutsD0toKpipipi.h"
  58 #include "AliAnalysisFilter.h"
  59 #include "AliAnalysisVertexingHF.h"
  60 #include "AliMixedEvent.h"
  61 #include "AliESDv0.h"
  62 #include "AliAODv0.h"
  63
  64 ClassImp(AliAnalysisVertexingHF)
  65
  66 //----------------------------------------------------------------------------
  67 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
  68 fInputAOD(kFALSE),
  69 fAODMapSize(0),
  70 fAODMap(0),
  71 fBzkG(0.),
  72 fSecVtxWithKF(kFALSE),
  73 fRecoPrimVtxSkippingTrks(kFALSE),
  74 fRmTrksFromPrimVtx(kFALSE),
  75 fV1(0x0),
  76 fD0toKpi(kTRUE),
  77 fJPSItoEle(kTRUE),
  78 f3Prong(kTRUE),
  79 f4Prong(kTRUE),
  80 fDstar(kTRUE),
  81 fCascades(kTRUE),
  82 fLikeSign(kFALSE),
  83 fMixEvent(kFALSE),
  84 fTrackFilter(0x0),
  85 fTrackFilterSoftPi(0x0),
  86 fCutsD0toKpi(0x0),
  87 fCutsJpsitoee(0x0),
  88 fCutsDplustoKpipi(0x0),
  89 fCutsDstoKKpi(0x0),
  90 fCutsLctopKpi(0x0),
  91 fCutsLctoV0(0x0),
  92 fCutsD0toKpipipi(0x0),
  93 fCutsD0fromDstar(0x0),
  94 fListOfCuts(0x0),
  95 fFindVertexForDstar(kTRUE),
  96 fFindVertexForCascades(kTRUE)
  97 {
  98   // Default constructor
  99
 100   SetD0toKpiCuts();
 101   SetBtoJPSICuts();
 102   SetDplusCuts();
 103   SetDsCuts();
 104   SetLcCuts();
 105   SetDstarCuts();
 106   SetD0to4ProngsCuts();
 107   SetLctoV0Cuts();
 108 }
 109 //--------------------------------------------------------------------------
 110 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
 111 TNamed(source),
 112 fInputAOD(source.fInputAOD),
 113 fAODMapSize(source.fAODMapSize),
 114 fAODMap(source.fAODMap),
 115 fBzkG(source.fBzkG),
 116 fSecVtxWithKF(source.fSecVtxWithKF),
 117 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
 118 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
 119 fV1(source.fV1),
 120 fD0toKpi(source.fD0toKpi),
 121 fJPSItoEle(source.fJPSItoEle),
 122 f3Prong(source.f3Prong),
 123 f4Prong(source.f4Prong),
 124 fDstar(source.fDstar),
 125 fCascades(source.fCascades),
 126 fLikeSign(source.fLikeSign),
 127 fMixEvent(source.fMixEvent),
 128 fTrackFilter(source.fTrackFilter),
 129 fTrackFilterSoftPi(source.fTrackFilterSoftPi),
 130 fCutsD0toKpi(source.fCutsD0toKpi),
 131 fCutsJpsitoee(source.fCutsJpsitoee),
 132 fCutsDplustoKpipi(source.fCutsDplustoKpipi),
 133 fCutsDstoKKpi(source.fCutsDstoKKpi),
 134 fCutsLctopKpi(source.fCutsLctopKpi),
 135 fCutsLctoV0(source.fCutsLctoV0),
 136 fCutsD0toKpipipi(source.fCutsD0toKpipipi),
 137 fCutsD0fromDstar(source.fCutsD0fromDstar),
 138 fListOfCuts(source.fListOfCuts),
 139 fFindVertexForDstar(source.fFindVertexForDstar),
 140 fFindVertexForCascades(source.fFindVertexForCascades)
 141 {
 142   //
 143   // Copy constructor
 144   //
 145   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
 146   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
 147   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
 148   for(Int_t i=0; i<14; i++) fDsCuts[i]=source.fDsCuts[i];
 149   for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
 150   for(Int_t i=0; i<8; i++)  fLctoV0Cuts[i]=source.fLctoV0Cuts[i];
 151   for(Int_t i=0; i<5; i++)  fDstarCuts[i]=source.fDstarCuts[i];
 152   for(Int_t i=0; i<9; i++)  fD0to4ProngsCuts[i]=source.fD0to4ProngsCuts[i];
 153 }
 154 //--------------------------------------------------------------------------
 155 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
 156 {
 157   //
 158   // assignment operator
 159   //
 160   if(&source == this) return *this;
 161   fInputAOD = source.fInputAOD;
 162   fAODMapSize = source.fAODMapSize;
 163   fBzkG = source.fBzkG;
 164   fSecVtxWithKF = source.fSecVtxWithKF;
 165   fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
 166   fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
 167   fV1 = source.fV1;
 168   fD0toKpi = source.fD0toKpi;
 169   fJPSItoEle = source.fJPSItoEle;
 170   f3Prong = source.f3Prong;
 171   f4Prong = source.f4Prong;
 172   fDstar = source.fDstar;
 173   fCascades = source.fCascades;
 174   fLikeSign = source.fLikeSign;
 175   fMixEvent = source.fMixEvent;
 176   fTrackFilter = source.fTrackFilter;
 177   fTrackFilterSoftPi = source.fTrackFilterSoftPi;
 178   fCutsD0toKpi = source.fCutsD0toKpi;
 179   fCutsJpsitoee = source.fCutsJpsitoee;
 180   fCutsDplustoKpipi = source.fCutsDplustoKpipi;
 181   fCutsDstoKKpi = source.fCutsDstoKKpi;
 182   fCutsLctopKpi = source.fCutsLctopKpi;
 183   fCutsLctoV0 = source.fCutsLctoV0;
 184   fCutsD0toKpipipi = source.fCutsD0toKpipipi;
 185   fCutsD0fromDstar = source.fCutsD0fromDstar;
 186   fListOfCuts = source.fListOfCuts;
 187   fFindVertexForDstar = source.fFindVertexForDstar;
 188   fFindVertexForCascades = source.fFindVertexForCascades;
 189
 190   for(Int_t i=0; i<9; i++)  fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
 191   for(Int_t i=0; i<9; i++)  fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
 192   for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
 193   for(Int_t i=0; i<14; i++) fDsCuts[i]=source.fDsCuts[i];
 194   for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
 195   for(Int_t i=0; i<8; i++)  fLctoV0Cuts[i]=source.fLctoV0Cuts[i];
 196   for(Int_t i=0; i<5; i++)  fDstarCuts[i]=source.fDstarCuts[i];
 197   for(Int_t i=0; i<9; i++)  fD0to4ProngsCuts[i]=source.fD0to4ProngsCuts[i];
 198
 199   return *this;
 200 }
 201 //----------------------------------------------------------------------------
 202 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
 203   // Destructor
 204   if(fV1) { delete fV1; fV1=0; }
 205   if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
 206   if(fTrackFilterSoftPi) { delete fTrackFilterSoftPi; fTrackFilterSoftPi=0; }
 207   if(fCutsD0toKpi) { delete fCutsD0toKpi; fCutsD0toKpi=0; }
 208   if(fCutsJpsitoee) { delete fCutsJpsitoee; fCutsJpsitoee=0; }
 209   if(fCutsDplustoKpipi) { delete fCutsDplustoKpipi; fCutsDplustoKpipi=0; }
 210   if(fCutsDstoKKpi) { delete fCutsDstoKKpi; fCutsDstoKKpi=0; }
 211   if(fCutsLctopKpi) { delete fCutsLctopKpi; fCutsLctopKpi=0; }
 212   if(fCutsLctoV0) { delete fCutsLctoV0; fCutsLctoV0=0; }
 213   if(fCutsD0toKpipipi) { delete fCutsD0toKpipipi; fCutsD0toKpipipi=0; }
 214   if(fCutsD0fromDstar) { delete fCutsD0fromDstar; fCutsD0fromDstar=0; }
 215   if(fAODMap) { delete fAODMap; fAODMap=0; }
 216 }
 217 //----------------------------------------------------------------------------
 218 TList *AliAnalysisVertexingHF::FillListOfCuts() {
 219   // Fill list of analysis cuts
 220
 221   TList *list = new TList();
 222   list->SetOwner();
 223   list->SetName("ListOfCuts");
 224
 225   if(fCutsD0toKpi) {
 226     AliRDHFCutsD0toKpi *cutsD0toKpi = new AliRDHFCutsD0toKpi(*fCutsD0toKpi);
 227     list->Add(cutsD0toKpi);
 228   }
 229   if(fCutsJpsitoee) {
 230     AliRDHFCutsJpsitoee *cutsJpsitoee = new AliRDHFCutsJpsitoee(*fCutsJpsitoee);
 231     list->Add(cutsJpsitoee);
 232   }
 233   if(fCutsDplustoKpipi) {
 234     AliRDHFCutsDplustoKpipi *cutsDplustoKpipi = new AliRDHFCutsDplustoKpipi(*fCutsDplustoKpipi);
 235     list->Add(cutsDplustoKpipi);
 236   }
 237   if(fCutsDstoKKpi) {
 238     AliRDHFCutsDstoKKpi *cutsDstoKKpi = new AliRDHFCutsDstoKKpi(*fCutsDstoKKpi);
 239     list->Add(cutsDstoKKpi);
 240   }
 241   if(fCutsLctopKpi) {
 242     AliRDHFCutsLctopKpi *cutsLctopKpi = new AliRDHFCutsLctopKpi(*fCutsLctopKpi);
 243     list->Add(cutsLctopKpi);
 244   }
 245   if(fCutsLctoV0){
 246     AliRDHFCutsLctoV0 *cutsLctoV0 = new AliRDHFCutsLctoV0(*fCutsLctoV0);
 247     list->Add(cutsLctoV0);
 248   }
 249   if(fCutsD0toKpipipi) {
 250     AliRDHFCutsD0toKpipipi *cutsD0toKpipipi = new AliRDHFCutsD0toKpipipi(*fCutsD0toKpipipi);
 251     list->Add(cutsD0toKpipipi);
 252   }
 253   if(fCutsD0fromDstar) {
 254     AliRDHFCutsD0toKpi *cutsD0fromDstar = new AliRDHFCutsD0toKpi(*fCutsD0fromDstar);
 255     list->Add(cutsD0fromDstar);
 256   }
 257
 258   // keep a pointer to the list
 259   fListOfCuts = list;
 260
 261   return list;
 262 }
 263 //----------------------------------------------------------------------------
 264 void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
 265                                             TClonesArray *aodVerticesHFTClArr,
 266                                             TClonesArray *aodD0toKpiTClArr,
 267                                             TClonesArray *aodJPSItoEleTClArr,
 268                                             TClonesArray *aodCharm3ProngTClArr,
 269                                             TClonesArray *aodCharm4ProngTClArr,
 270                                             TClonesArray *aodDstarTClArr,
 271                                             TClonesArray *aodCascadesTClArr,
 272                                             TClonesArray *aodLikeSign2ProngTClArr,
 273                                             TClonesArray *aodLikeSign3ProngTClArr)
 274 {
 275   // Find heavy-flavour vertex candidates
 276   // Input:  ESD or AOD
 277   // Output: AOD (additional branches added)
 278
 279   if(!fMixEvent){
 280     TString evtype = event->IsA()->GetName();
 281     fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
 282   } // if we do mixing AliVEvent is a AliMixedEvent
 283
 284   if(fInputAOD) {
 285     AliDebug(2,"Creating HF candidates from AOD");
 286   } else {
 287     AliDebug(2,"Creating HF candidates from ESD");
 288   }
 289
 290   if(!aodVerticesHFTClArr) {
 291     printf("ERROR: no aodVerticesHFTClArr");
 292     return;
 293   }
 294   if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
 295     printf("ERROR: no aodD0toKpiTClArr");
 296     return;
 297   }
 298   if(fJPSItoEle && !aodJPSItoEleTClArr) {
 299     printf("ERROR: no aodJPSItoEleTClArr");
 300     return;
 301   }
 302   if(f3Prong && !aodCharm3ProngTClArr) {
 303     printf("ERROR: no aodCharm3ProngTClArr");
 304     return;
 305   }
 306   if(f4Prong && !aodCharm4ProngTClArr) {
 307     printf("ERROR: no aodCharm4ProngTClArr");
 308     return;
 309   }
 310   if(fDstar && !aodDstarTClArr) {
 311     printf("ERROR: no aodDstarTClArr");
 312     return;
 313   }
 314   if(fCascades && !aodCascadesTClArr){
 315     printf("ERROR: no aodCascadesTClArr ");
 316     return;
 317   }
 318   if(fLikeSign && !aodLikeSign2ProngTClArr) {
 319     printf("ERROR: no aodLikeSign2ProngTClArr");
 320     return;
 321   }
 322   if(fLikeSign && f3Prong && !aodLikeSign3ProngTClArr) {
 323     printf("ERROR: no aodLikeSign2ProngTClArr");
 324     return;
 325   }
 326
 327   // delete candidates from previous event and create references
 328   Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0,iCascades=0,iLikeSign2Prong=0,iLikeSign3Prong=0;
 329   aodVerticesHFTClArr->Delete();
 330   iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
 331   TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
 332   if(fD0toKpi || fDstar)   {
 333     aodD0toKpiTClArr->Delete();
 334     iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
 335   }
 336   if(fJPSItoEle) {
 337     aodJPSItoEleTClArr->Delete();
 338     iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
 339   }
 340   if(f3Prong) {
 341     aodCharm3ProngTClArr->Delete();
 342     i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
 343   }
 344   if(f4Prong) {
 345     aodCharm4ProngTClArr->Delete();
 346     i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
 347   }
 348   if(fDstar) {
 349     aodDstarTClArr->Delete();
 350     iDstar = aodDstarTClArr->GetEntriesFast();
 351   }
 352   if(fCascades) {
 353     aodCascadesTClArr->Delete();
 354     iCascades = aodCascadesTClArr->GetEntriesFast();
 355   }
 356   if(fLikeSign) {
 357     aodLikeSign2ProngTClArr->Delete();
 358     iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast();
 359   }
 360   if(fLikeSign && f3Prong) {
 361     aodLikeSign3ProngTClArr->Delete();
 362     iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast();
 363   }
 364
 365   TClonesArray &aodD0toKpiRef        = *aodD0toKpiTClArr;
 366   TClonesArray &aodJPSItoEleRef      = *aodJPSItoEleTClArr;
 367   TClonesArray &aodCharm3ProngRef    = *aodCharm3ProngTClArr;
 368   TClonesArray &aodCharm4ProngRef    = *aodCharm4ProngTClArr;
 369   TClonesArray &aodDstarRef          = *aodDstarTClArr;
 370   TClonesArray &aodCascadesRef       = *aodCascadesTClArr;
 371   TClonesArray &aodLikeSign2ProngRef = *aodLikeSign2ProngTClArr;
 372   TClonesArray &aodLikeSign3ProngRef = *aodLikeSign3ProngTClArr;
 373
 374
 375   AliAODRecoDecayHF2Prong *io2Prong  = 0;
 376   AliAODRecoDecayHF3Prong *io3Prong  = 0;
 377   AliAODRecoDecayHF4Prong *io4Prong  = 0;
 378   AliAODRecoCascadeHF     *ioCascade = 0;
 379
 380   Int_t    iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries,iv0,nv0;
 381   Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaV0,dcaCasc;
 382   Bool_t   okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
 383   Bool_t   okDstar=kFALSE,okD0fromDstar=kFALSE;
 384   Bool_t   okCascades=kFALSE;
 385   AliESDtrack *postrack1 = 0;
 386   AliESDtrack *postrack2 = 0;
 387   AliESDtrack *negtrack1 = 0;
 388   AliESDtrack *negtrack2 = 0;
 389   AliESDtrack *trackPi   = 0;
 390 //   AliESDtrack *posV0track = 0;
 391 //   AliESDtrack *negV0track = 0;
 392   /*
 393   Double_t dcaMax = fD0toKpiCuts[1];
 394   if(dcaMax < fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1];
 395   if(dcaMax < fDplusCuts[11])  dcaMax=fDplusCuts[11];
 396   if(dcaMax < fD0to4ProngsCuts[1])  dcaMax=fD0to4ProngsCuts[1];
 397   */
 398   Float_t dcaMax = fCutsD0toKpi->GetDCACut();
 399   if(fCutsJpsitoee) dcaMax=TMath::Max(dcaMax,fCutsJpsitoee->GetDCACut());
 400   if(fCutsDplustoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDplustoKpipi->GetDCACut());
 401   if(fCutsDstoKKpi) dcaMax=TMath::Max(dcaMax,fCutsDstoKKpi->GetDCACut());
 402   if(fCutsLctopKpi) dcaMax=TMath::Max(dcaMax,fCutsLctopKpi->GetDCACut());
 403   if(fCutsD0toKpipipi) dcaMax=TMath::Max(dcaMax,fCutsD0toKpipipi->GetDCACut());
 404
 405   AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
 406
 407
 408   // get Bz
 409   fBzkG = (Double_t)event->GetMagneticField();
 410
 411   trkEntries = (Int_t)event->GetNumberOfTracks();
 412   AliDebug(1,Form(" Number of tracks: %d",trkEntries));
 413
 414   nv0 = (Int_t)event->GetNumberOfV0s();
 415   AliDebug(1,Form(" Number of V0s: %d",nv0));
 416
 417   if( trkEntries<2 && (trkEntries<1 || nv0<1) ) {
 418     AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
 419     return;
 420   }
 421
 422   // event selection
 423   if(!fCutsD0toKpi->IsEventSelected(event)) return;
 424
 425   // call function that applies sigle-track selection,
 426   // for displaced tracks and soft pions (both charges) for D*,
 427   // and retrieves primary vertex
 428   TObjArray seleTrksArray(trkEntries);
 429   UChar_t  *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi
 430   Int_t     nSeleTrks=0;
 431   Int_t *evtNumber    = new Int_t[trkEntries];
 432   SelectTracksAndCopyVertex(event,seleTrksArray,nSeleTrks,seleFlags,evtNumber);
 433
 434   AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
 435
 436   TObjArray *twoTrackArray1    = new TObjArray(2);
 437   TObjArray *twoTrackArray2    = new TObjArray(2);
 438   TObjArray *twoTrackArrayV0   = new TObjArray(2);
 439   TObjArray *twoTrackArrayCasc = new TObjArray(2);
 440   TObjArray *threeTrackArray   = new TObjArray(3);
 441   TObjArray *fourTrackArray    = new TObjArray(4);
 442
 443   Double_t dispersion;
 444   Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
 445
 446   AliAODRecoDecayHF   *rd = 0;
 447   AliAODRecoCascadeHF *rc = 0;
 448   AliAODv0            *V0 = 0;
 449   AliESDv0         *esdV0 = 0;
 450
 451   // LOOP ON  POSITIVE  TRACKS
 452   for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
 453
 454     if(iTrkP1%1==0) AliDebug(1,Form("  1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));
 455
 456     // get track from tracks array
 457     postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
 458
 459     if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
 460
 461     // LOOP ON v0s here
 462     //
 463     if(fCascades)
 464       for(iv0=0; iv0<nv0; iv0++){
 465
 466         AliDebug(1,Form("   loop on v0s for track number %d and v0 number %d",iTrkP1,iv0));
 467
 468         // Get the V0
 469         if(fInputAOD) V0 = ((AliAODEvent*)event)->GetV0(iv0);
 470         else {
 471           esdV0 = ((AliESDEvent*)event)->GetV0(iv0);
 472         }
 473         if ( (!V0 || !V0->IsA()->InheritsFrom("AliAODv0") ) &&
 474              (!esdV0 || !esdV0->IsA()->InheritsFrom("AliESDv0") ) )
 475           continue;
 476
 477
 478         // Get the tracks that form the V0
 479         //  ( parameters at primary vertex )
 480         //   and define an AliExternalTrackParam out of them
 481         AliExternalTrackParam * posV0track;
 482         AliExternalTrackParam * negV0track;
 483
 484         if(fInputAOD){
 485           AliAODTrack *posVV0track = (AliAODTrack*)(V0->GetDaughter(0));
 486           AliAODTrack *negVV0track = (AliAODTrack*)(V0->GetDaughter(1));
 487           if( !posVV0track || !negVV0track ) continue;
 488           //
 489           // Apply some basic V0 daughter criteria
 490           //
 491           // bachelor must not be a v0-track
 492           if (posVV0track->GetID() == postrack1->GetID() ||
 493               negVV0track->GetID() == postrack1->GetID()) continue;
 494           // reject like-sign v0
 495           if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
 496           // avoid ghost TPC tracks
 497           if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
 498              !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
 499           // Get AliExternalTrackParam out of the AliAODTracks
 500           Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign;
 501           posVV0track->PxPyPz(pxpypz);                    posVV0track->XvYvZv(xyz);
 502           posVV0track->GetCovarianceXYZPxPyPz(cv);        sign=posVV0track->Charge();
 503           posV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
 504           negVV0track->PxPyPz(pxpypz);                    negVV0track->XvYvZv(xyz);
 505           negVV0track->GetCovarianceXYZPxPyPz(cv);        sign=negVV0track->Charge();
 506           negV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
 507         }
 508         else {
 509           AliESDtrack *posVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetPindex() ));
 510           AliESDtrack *negVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetNindex() ));
 511           if( !posVV0track || !negVV0track ) continue;
 512           //
 513           // Apply some basic V0 daughter criteria
 514           //
 515           // bachelor must not be a v0-track
 516           if (posVV0track->GetID() == postrack1->GetID() ||
 517               negVV0track->GetID() == postrack1->GetID()) continue;
 518           // reject like-sign v0
 519           if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
 520           // avoid ghost TPC tracks
 521           if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
 522              !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
 523           //  reject kinks (only necessary on AliESDtracks)
 524           if (posVV0track->GetKinkIndex(0)>0  || negVV0track->GetKinkIndex(0)>0) continue;
 525           // Get AliExternalTrackParam out of the AliESDtracks
 526           posV0track = dynamic_cast<AliExternalTrackParam*>(posVV0track);
 527           negV0track = dynamic_cast<AliExternalTrackParam*>(negVV0track);
 528         }
 529         if( !posV0track || !negV0track ){
 530           AliDebug(1,Form(" Couldn't get the V0 daughters"));
 531           continue;
 532         }
 533
 534         // fill in the v0 two-external-track-param array
 535         twoTrackArrayV0->AddAt(posV0track,0);
 536         twoTrackArrayV0->AddAt(negV0track,1);
 537
 538         // Define the AODv0 from ESDv0 if reading ESDs
 539         if(!fInputAOD) V0 = TransformESDv0toAODv0(esdV0,twoTrackArrayV0);
 540
 541         // Get the V0 dca
 542         dcaV0 = V0->DcaV0Daughters();
 543
 544         // Define the V0 (neutral) track
 545         AliNeutralTrackParam *trackV0;
 546         if(fInputAOD) {
 547           const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(V0);
 548           if(!trackVV0) continue;
 549           trackV0 = new AliNeutralTrackParam(trackVV0);
 550         }
 551         else{
 552           Double_t xyz[3], pxpypz[3];
 553           esdV0->XvYvZv(xyz);
 554           esdV0->PxPyPz(pxpypz);
 555           Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
 556           trackV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
 557         }
 558         if(!trackV0){
 559           AliDebug(1, Form("Couldn't define the V0 as a neutral track !! \n"));
 560           continue;
 561         }
 562         // Fill in the object array to create the cascade
 563         twoTrackArrayCasc->AddAt(postrack1,0);
 564         twoTrackArrayCasc->AddAt(trackV0,1);
 565
 566         // Compute the cascade vertex
 567         AliAODVertex *vertexCasc = 0;
 568         if(fFindVertexForCascades) {
 569           // DCA between the two tracks
 570           dcaCasc = postrack1->GetDCA(trackV0,fBzkG,xdummy,ydummy);
 571           // Vertexing+
 572           vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
 573         } else {
 574           // assume Cascade decays at the primary vertex
 575           Double_t pos[3],cov[6],chi2perNDF;
 576           fV1->GetXYZ(pos);
 577           fV1->GetCovMatrix(cov);
 578           chi2perNDF = fV1->GetChi2toNDF();
 579           vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
 580           dcaCasc = 0.;
 581         }
 582         if(!vertexCasc) {
 583           twoTrackArrayCasc->Clear();
 584           continue;
 585         }
 586
 587         // Create and store the Cascade if passed the cuts
 588         ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,V0,dcaCasc,okCascades);
 589         if(okCascades && ioCascade) {
 590           AliDebug(1,Form("Storing a cascade object... "));
 591           // add the vertex and the cascade to the AOD
 592           AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
 593           rc = new(aodCascadesRef[iCascades++])AliAODRecoCascadeHF(*ioCascade);
 594           rc->SetSecondaryVtx(vCasc);
 595           vCasc->SetParent(rc);
 596           rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
 597           if(!fInputAOD) vCasc->AddDaughter(V0); // just to fill ref #0 ??
 598           AddRefs(vCasc,rc,event,twoTrackArrayCasc); // add the track (proton)
 599           vCasc->AddDaughter(V0); // fill the 2prong V0
 600         }
 601
 602         // Clean up
 603         twoTrackArrayV0->Clear();
 604         twoTrackArrayCasc->Clear();
 605         if(ioCascade) { delete ioCascade; ioCascade=NULL; }
 606         if(vertexCasc) { delete vertexCasc; vertexCasc=NULL; }
 607       }
 608
 609
 610     // If there is less than 2 particles exit
 611     if(trkEntries<2) {
 612       AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
 613       return;
 614     }
 615
 616     if(postrack1->Charge()<0 && !fLikeSign) continue;
 617
 618     // LOOP ON  NEGATIVE  TRACKS
 619     for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
 620
 621       if(iTrkN1%1==0) AliDebug(1,Form("    1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));
 622
 623       if(iTrkN1==iTrkP1) continue;
 624
 625       // get track from tracks array
 626       negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
 627
 628       if(negtrack1->Charge()>0 && !fLikeSign) continue;
 629
 630       if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
 631
 632       if(fMixEvent) {
 633         if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
 634       }
 635
 636       if(postrack1->Charge()==negtrack1->Charge()) { // like-sign
 637         isLikeSign2Prong=kTRUE;
 638         if(!fLikeSign)    continue;
 639         if(iTrkN1<iTrkP1) continue; // this is needed to avoid double-counting of like-sign
 640       } else { // unlike-sign
 641         isLikeSign2Prong=kFALSE;
 642         if(postrack1->Charge()<0 || negtrack1->Charge()>0) continue;  // this is needed to avoid double-counting of unlike-sign
 643         if(fMixEvent) {
 644           if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
 645         }
 646
 647       }
 648
 649       // back to primary vertex
 650       postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 651       negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 652
 653       // DCA between the two tracks
 654       dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 655       if(dcap1n1>dcaMax) { negtrack1=0; continue; }
 656
 657       // Vertexing
 658       twoTrackArray1->AddAt(postrack1,0);
 659       twoTrackArray1->AddAt(negtrack1,1);
 660       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
 661       if(!vertexp1n1) {
 662         twoTrackArray1->Clear();
 663         negtrack1=0;
 664         continue;
 665       }
 666
 667       // 2 prong candidate
 668       if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) {
 669
 670         io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
 671
 672         if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI) || (isLikeSign2Prong && (okD0 || okJPSI))) {
 673           // add the vertex and the decay to the AOD
 674           AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
 675           if(!isLikeSign2Prong) {
 676             if(okD0) {
 677               rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
 678               rd->SetSecondaryVtx(v2Prong);
 679               v2Prong->SetParent(rd);
 680               AddRefs(v2Prong,rd,event,twoTrackArray1);
 681             }
 682             if(okJPSI) {
 683               rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
 684               rd->SetSecondaryVtx(v2Prong);
 685               if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
 686               AddRefs(v2Prong,rd,event,twoTrackArray1);
 687             }
 688           } else { // isLikeSign2Prong
 689             rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
 690             rd->SetSecondaryVtx(v2Prong);
 691             v2Prong->SetParent(rd);
 692             AddRefs(v2Prong,rd,event,twoTrackArray1);
 693           }
 694         }
 695         // D* candidates
 696         if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
 697           // write references in io2Prong
 698           if(fInputAOD) {
 699             AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
 700           } else {
 701             vertexp1n1->AddDaughter(postrack1);
 702             vertexp1n1->AddDaughter(negtrack1);
 703           }
 704           io2Prong->SetSecondaryVtx(vertexp1n1);
 705           //printf("--->  %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
 706           // create a track from the D0
 707           AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
 708
 709           // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
 710           for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
 711
 712             if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
 713
 714             if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
 715
 716             if(fMixEvent) {
 717               if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] ||
 718                  evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] ||
 719                  evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
 720             }
 721
 722             if(iTrkSoftPi%1==0) AliDebug(1,Form("    1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));
 723
 724             trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
 725             if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
 726
 727             // get track from tracks array
 728             trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
 729             trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
 730
 731             twoTrackArrayCasc->AddAt(trackPi,0);
 732             twoTrackArrayCasc->AddAt(trackD0,1);
 733
 734             AliAODVertex *vertexCasc = 0;
 735
 736             if(fFindVertexForDstar) {
 737               // DCA between the two tracks
 738               dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
 739               // Vertexing
 740               vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
 741             } else {
 742               // assume Dstar decays at the primary vertex
 743               Double_t pos[3],cov[6],chi2perNDF;
 744               fV1->GetXYZ(pos);
 745               fV1->GetCovMatrix(cov);
 746               chi2perNDF = fV1->GetChi2toNDF();
 747               vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
 748               dcaCasc = 0.;
 749             }
 750             if(!vertexCasc) {
 751               twoTrackArrayCasc->Clear();
 752               trackPi=0;
 753               continue;
 754             }
 755
 756             ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
 757             if(okDstar) {
 758               // add the D0 to the AOD (if not already done)
 759               if(!okD0) {
 760                 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
 761                 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
 762                 rd->SetSecondaryVtx(v2Prong);
 763                 v2Prong->SetParent(rd);
 764                 AddRefs(v2Prong,rd,event,twoTrackArray1);
 765                 okD0=kTRUE; // this is done to add it only once
 766               }
 767               // add the vertex and the cascade to the AOD
 768               AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
 769               rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
 770               rc->SetSecondaryVtx(vCasc);
 771               vCasc->SetParent(rc);
 772               if(!fInputAOD) vCasc->AddDaughter(rd); // just to fill ref #0
 773               AddRefs(vCasc,rc,event,twoTrackArrayCasc);
 774               vCasc->AddDaughter(rd); // add the D0 (in ref #1)
 775             }
 776             twoTrackArrayCasc->Clear();
 777             trackPi=0;
 778             if(ioCascade) {delete ioCascade; ioCascade=NULL;}
 779             delete vertexCasc; vertexCasc=NULL;
 780           } // end loop on soft pi tracks
 781
 782           if(trackD0) {delete trackD0; trackD0=NULL;}
 783
 784         }
 785         if(io2Prong) {delete io2Prong; io2Prong=NULL;}
 786       }
 787
 788       twoTrackArray1->Clear();
 789       if( (!f3Prong && !f4Prong) ||
 790           (isLikeSign2Prong && !f3Prong) ) {
 791         negtrack1=0;
 792         delete vertexp1n1;
 793         continue;
 794       }
 795
 796
 797       // 2nd LOOP  ON  POSITIVE  TRACKS
 798       for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
 799
 800         if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
 801
 802         if(iTrkP2%1==0) AliDebug(1,Form("    2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));
 803
 804         // get track from tracks array
 805         postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
 806
 807         if(postrack2->Charge()<0) continue;
 808
 809         if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
 810
 811         if(fMixEvent) {
 812           if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
 813              evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
 814              evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
 815         }
 816
 817         if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
 818           if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
 819             isLikeSign3Prong=kTRUE;
 820           } else { // not ok
 821             continue;
 822           }
 823         } else { // normal triplet (+-+)
 824           isLikeSign3Prong=kFALSE;
 825           if(fMixEvent) {
 826             if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
 827                evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
 828                evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
 829           }
 830         }
 831
 832         // back to primary vertex
 833         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 834         postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 835         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 836         //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
 837
 838         dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 839         if(dcap2n1>dcaMax) { postrack2=0; continue; }
 840         dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
 841         if(dcap1p2>dcaMax) { postrack2=0; continue; }
 842
 843         // Vertexing
 844         twoTrackArray2->AddAt(postrack2,0);
 845         twoTrackArray2->AddAt(negtrack1,1);
 846         AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
 847         if(!vertexp2n1) {
 848           twoTrackArray2->Clear();
 849           postrack2=0;
 850           continue;
 851         }
 852
 853         // 3 prong candidates
 854         if(f3Prong) {
 855           if(postrack2->Charge()>0) {
 856             threeTrackArray->AddAt(postrack1,0);
 857             threeTrackArray->AddAt(negtrack1,1);
 858             threeTrackArray->AddAt(postrack2,2);
 859           } else {
 860             threeTrackArray->AddAt(negtrack1,0);
 861             threeTrackArray->AddAt(postrack1,1);
 862             threeTrackArray->AddAt(postrack2,2);
 863           }
 864
 865           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
 866           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
 867           if(ok3Prong) {
 868             AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
 869             if(!isLikeSign3Prong) {
 870               rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
 871               rd->SetSecondaryVtx(v3Prong);
 872               v3Prong->SetParent(rd);
 873               AddRefs(v3Prong,rd,event,threeTrackArray);
 874             } else { // isLikeSign3Prong
 875               rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
 876               rd->SetSecondaryVtx(v3Prong);
 877               v3Prong->SetParent(rd);
 878               AddRefs(v3Prong,rd,event,threeTrackArray);
 879             }
 880           }
 881           if(io3Prong) {delete io3Prong; io3Prong=NULL;}
 882           if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
 883         }
 884
 885         // 4 prong candidates
 886         if(f4Prong
 887            // don't make 4 prong with like-sign pairs and triplets
 888            && !isLikeSign2Prong && !isLikeSign3Prong
 889            // track-to-track dca cuts already now
 890            && dcap1n1 < fD0to4ProngsCuts[1]
 891            && dcap2n1 < fD0to4ProngsCuts[1]) {
 892
 893           // back to primary vertex
 894           postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 895           postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 896           negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 897           // Vertexing for these 3 (can be taken from above?)
 898           threeTrackArray->AddAt(postrack1,0);
 899           threeTrackArray->AddAt(negtrack1,1);
 900           threeTrackArray->AddAt(postrack2,2);
 901           AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
 902
 903           // 3rd LOOP  ON  NEGATIVE  TRACKS (for 4 prong)
 904           for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
 905
 906             if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
 907
 908             if(iTrkN2%1==0) AliDebug(1,Form("    3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
 909
 910             // get track from tracks array
 911             negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
 912
 913             if(negtrack2->Charge()>0) continue;
 914
 915             if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
 916             if(fMixEvent){
 917               if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
 918                  evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
 919                  evtNumber[iTrkP2]==evtNumber[iTrkN2] ||
 920                  evtNumber[iTrkP1]==evtNumber[iTrkN1] ||
 921                  evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
 922                  evtNumber[iTrkN1]==evtNumber[iTrkP2]) continue;
 923             }
 924
 925             // back to primary vertex
 926             postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 927             postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 928             negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 929             negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 930             dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 931             if(dcap1n2 > fD0to4ProngsCuts[1]) { negtrack2=0; continue; }
 932             dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 933             if(dcap2n2 > fD0to4ProngsCuts[1]) { negtrack2=0; continue; }
 934
 935             // Vertexing
 936             fourTrackArray->AddAt(postrack1,0);
 937             fourTrackArray->AddAt(negtrack1,1);
 938             fourTrackArray->AddAt(postrack2,2);
 939             fourTrackArray->AddAt(negtrack2,3);
 940
 941             AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
 942             io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
 943             if(ok4Prong) {
 944               AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
 945               rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
 946               rd->SetSecondaryVtx(v4Prong);
 947               v4Prong->SetParent(rd);
 948               AddRefs(v4Prong,rd,event,fourTrackArray);
 949             }
 950
 951             if(io4Prong) {delete io4Prong; io4Prong=NULL;}
 952             if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;}
 953             fourTrackArray->Clear();
 954             negtrack2 = 0;
 955
 956           } // end loop on negative tracks
 957
 958           threeTrackArray->Clear();
 959           delete vertexp1n1p2;
 960
 961         }
 962
 963         postrack2 = 0;
 964         delete vertexp2n1;
 965
 966       } // end 2nd loop on positive tracks
 967
 968       twoTrackArray2->Clear();
 969
 970       // 2nd LOOP  ON  NEGATIVE  TRACKS (for 3 prong -+-)
 971       for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
 972
 973         if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
 974
 975         if(iTrkN2%1==0) AliDebug(1,Form("    2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
 976
 977         // get track from tracks array
 978         negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
 979
 980         if(negtrack2->Charge()>0) continue;
 981
 982         if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
 983
 984         if(fMixEvent) {
 985           if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
 986              evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
 987              evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
 988         }
 989
 990         if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
 991           if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
 992             isLikeSign3Prong=kTRUE;
 993           } else { // not ok
 994             continue;
 995           }
 996         } else { // normal triplet (-+-)
 997           isLikeSign3Prong=kFALSE;
 998         }
 999
1000         // back to primary vertex
1001         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1002         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1003         negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1004         //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
1005
1006         dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1007         if(dcap1n2>dcaMax) { negtrack2=0; continue; }
1008         dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1009         if(dcan1n2>dcaMax) { negtrack2=0; continue; }
1010
1011         // Vertexing
1012         twoTrackArray2->AddAt(postrack1,0);
1013         twoTrackArray2->AddAt(negtrack2,1);
1014
1015         AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1016         if(!vertexp1n2) {
1017           twoTrackArray2->Clear();
1018           negtrack2=0;
1019           continue;
1020         }
1021
1022         if(f3Prong) {
1023           threeTrackArray->AddAt(negtrack1,0);
1024           threeTrackArray->AddAt(postrack1,1);
1025           threeTrackArray->AddAt(negtrack2,2);
1026           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1027           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
1028           if(ok3Prong) {
1029             AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1030             if(!isLikeSign3Prong) {
1031               rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1032               rd->SetSecondaryVtx(v3Prong);
1033               v3Prong->SetParent(rd);
1034               AddRefs(v3Prong,rd,event,threeTrackArray);
1035             } else { // isLikeSign3Prong
1036               rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1037               rd->SetSecondaryVtx(v3Prong);
1038               v3Prong->SetParent(rd);
1039               AddRefs(v3Prong,rd,event,threeTrackArray);
1040             }
1041           }
1042           if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1043           if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1044         }
1045
1046         negtrack2 = 0;
1047         delete vertexp1n2;
1048
1049       } // end 2nd loop on negative tracks
1050
1051       twoTrackArray2->Clear();
1052
1053       negtrack1 = 0;
1054       delete vertexp1n1;
1055     } // end 1st loop on negative tracks
1056
1057     postrack1 = 0;
1058   }  // end 1st loop on positive tracks
1059
1060
1061   AliDebug(1,Form(" Total HF vertices in event = %d;",
1062                   (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
1063   if(fD0toKpi) {
1064     AliDebug(1,Form(" D0->Kpi in event = %d;",
1065                     (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
1066   }
1067   if(fJPSItoEle) {
1068     AliDebug(1,Form(" JPSI->ee in event = %d;",
1069                     (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
1070   }
1071   if(f3Prong) {
1072     AliDebug(1,Form(" Charm->3Prong in event = %d;",
1073                     (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
1074   }
1075   if(f4Prong) {
1076     AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
1077                     (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
1078   }
1079   if(fDstar) {
1080     AliDebug(1,Form(" D*->D0pi in event = %d;\n",
1081                     (Int_t)aodDstarTClArr->GetEntriesFast()));
1082   }
1083   if(fCascades){
1084     AliDebug(1,Form(" cascades -> v0 + track in event = %d;\n",
1085                     (Int_t)aodCascadesTClArr->GetEntriesFast()));
1086   }
1087   if(fLikeSign) {
1088     AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
1089                     (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
1090   }
1091   if(fLikeSign && f3Prong) {
1092     AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
1093                     (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
1094   }
1095
1096
1097   twoTrackArray1->Delete();  delete twoTrackArray1;
1098   twoTrackArray2->Delete();  delete twoTrackArray2;
1099   twoTrackArrayCasc->Delete();  delete twoTrackArrayCasc;
1100   twoTrackArrayV0->Delete();  delete twoTrackArrayV0;
1101   threeTrackArray->Clear();
1102   threeTrackArray->Delete(); delete threeTrackArray;
1103   fourTrackArray->Delete();  delete fourTrackArray;
1104   delete [] seleFlags; seleFlags=NULL;
1105   if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
1106
1107   if(fInputAOD) {
1108     seleTrksArray.Delete();
1109     if(fAODMap) { delete fAODMap; fAODMap=NULL; }
1110   }
1111
1112   return;
1113 }
1114 //----------------------------------------------------------------------------
1115 void AliAnalysisVertexingHF::AddRefs(AliAODVertex *v,AliAODRecoDecayHF *rd,
1116                                      const AliVEvent *event,
1117                                      const TObjArray *trkArray) const
1118 {
1119   // Add the AOD tracks as daughters of the vertex (TRef)
1120   // Also add the references to the primary vertex and to the cuts
1121
1122   if(fInputAOD) {
1123     AddDaughterRefs(v,event,trkArray);
1124     rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1125   }
1126
1127   /*
1128   rd->SetListOfCutsRef((TList*)fListOfCuts);
1129   //fListOfCuts->Print();
1130   cout<<fListOfCuts<<endl;
1131   TList *l=(TList*)rd->GetListOfCuts();
1132   cout<<l<<endl;
1133   if(l) {l->Print(); }else{printf("error\n");}
1134   */
1135
1136   return;
1137 }
1138 //----------------------------------------------------------------------------
1139 void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,
1140                                              const AliVEvent *event,
1141                                              const TObjArray *trkArray) const
1142 {
1143   // Add the AOD tracks as daughters of the vertex (TRef)
1144
1145   Int_t nDg = v->GetNDaughters();
1146   TObject *dg = 0;
1147   if(nDg) dg = v->GetDaughter(0);
1148
1149   if(dg) return; // daughters already added
1150
1151   Int_t nTrks = trkArray->GetEntriesFast();
1152
1153   AliExternalTrackParam *track = 0;
1154   AliAODTrack *aodTrack = 0;
1155   Int_t id;
1156
1157   for(Int_t i=0; i<nTrks; i++) {
1158     track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1159     id = (Int_t)track->GetID();
1160     //printf("---> %d\n",id);
1161     if(id<0) continue; // this track is a AliAODRecoDecay
1162     aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
1163     v->AddDaughter(aodTrack);
1164   }
1165
1166   return;
1167 }
1168 //----------------------------------------------------------------------------
1169 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1170                                    TObjArray *twoTrackArray,AliVEvent *event,
1171                                    AliAODVertex *secVert,
1172                                    AliAODRecoDecayHF2Prong *rd2Prong,
1173                                    Double_t dca,
1174                                    Bool_t &okDstar) const
1175 {
1176   // Make the cascade as a 2Prong decay and check if it passes Dstar
1177   // reconstruction cuts
1178
1179   okDstar = kFALSE;
1180
1181   Bool_t dummy1,dummy2,dummy3;
1182
1183   // We use Make2Prong to construct the AliAODRecoCascadeHF
1184   // (which inherits from AliAODRecoDecayHF2Prong)
1185   AliAODRecoCascadeHF *theCascade =
1186     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1187                                      dummy1,dummy2,dummy3);
1188   if(!theCascade) return 0x0;
1189
1190   // charge
1191   AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1192   theCascade->SetCharge(trackPi->Charge());
1193
1194   //--- selection cuts
1195   //
1196   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1197   tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
1198   tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
1199   AliAODVertex *primVertexAOD=0;
1200   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1201     // take event primary vertex
1202     primVertexAOD = PrimaryVertex();
1203     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1204     rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
1205   }
1206   // select D*->D0pi
1207   if(fDstar) {
1208     Bool_t testD0=kTRUE;
1209     okDstar = tmpCascade->SelectDstar(fDstarCuts,fD0fromDstarCuts,testD0);
1210   }
1211   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1212   tmpCascade->UnsetOwnPrimaryVtx();
1213   delete tmpCascade; tmpCascade=NULL;
1214   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1215     rd2Prong->UnsetOwnPrimaryVtx();
1216   }
1217   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1218   //---
1219
1220   return theCascade;
1221 }
1222
1223
1224 //----------------------------------------------------------------------------
1225 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1226                                    TObjArray *twoTrackArray,AliVEvent *event,
1227                                    AliAODVertex *secVert,
1228                                    AliAODv0 *v0,
1229                                    Double_t dca,
1230                                    Bool_t &okCascades) const
1231 {
1232   //
1233   // Make the cascade as a 2Prong decay and check if it passes
1234   // cascades reconstruction cuts
1235
1236   //  AliDebug(2,Form("         building the cascade"));
1237   okCascades= kFALSE;
1238   Bool_t dummy1,dummy2,dummy3;
1239
1240   // We use Make2Prong to construct the AliAODRecoCascadeHF
1241   // (which inherits from AliAODRecoDecayHF2Prong)
1242   AliAODRecoCascadeHF *theCascade =
1243     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1244                                      dummy1,dummy2,dummy3);
1245   if(!theCascade) return 0x0;
1246
1247   // bachelor track and charge
1248   AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1249   theCascade->SetCharge(trackBachelor->Charge());
1250
1251   //--- selection cuts
1252   //
1253   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1254   tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
1255   tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
1256   AliAODVertex *primVertexAOD=0;
1257   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1258     // take event primary vertex
1259     primVertexAOD = PrimaryVertex();
1260     if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
1261     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1262   }
1263
1264   // select Cascades
1265   bool okLcksp=0, okLcLpi=0;
1266   if(fCascades && fInputAOD){
1267     if(fCutsLctoV0) {
1268       okCascades = (bool)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1269       if(okCascades==1) okLcksp=1;
1270       if(okCascades==2) okLcLpi=1;
1271       if(okCascades==3) { okLcksp=1; okLcLpi=1;}
1272     }
1273     else okCascades = tmpCascade->SelectLctoV0(fLctoV0Cuts,okLcksp,okLcLpi);
1274   }
1275   else { AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied")); okCascades=true; }// no cuts implemented from ESDs
1276   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1277   tmpCascade->UnsetOwnPrimaryVtx();
1278   delete tmpCascade; tmpCascade=NULL;
1279   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1280   //---
1281
1282   return theCascade;
1283 }
1284
1285 //-----------------------------------------------------------------------------
1286 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
1287                                    TObjArray *twoTrackArray,AliVEvent *event,
1288                                    AliAODVertex *secVert,Double_t dca,
1289                                    Bool_t &okD0,Bool_t &okJPSI,
1290                                    Bool_t &okD0fromDstar) const
1291 {
1292   // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
1293   // reconstruction cuts
1294   // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
1295
1296   okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
1297
1298   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1299   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1300   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
1301
1302   // propagate tracks to secondary vertex, to compute inv. mass
1303   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1304   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1305
1306   Double_t momentum[3];
1307   postrack->GetPxPyPz(momentum);
1308   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1309   negtrack->GetPxPyPz(momentum);
1310   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1311
1312
1313   // invariant mass cut (try to improve coding here..)
1314   Bool_t okMassCut=kFALSE;
1315   if(!okMassCut && fD0toKpi)   if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
1316   if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
1317   if(!okMassCut && fDstar)     if(SelectInvMass(3,2,px,py,pz)) okMassCut=kTRUE;
1318   if(!okMassCut) {
1319     AliDebug(2," candidate didn't pass mass cut");
1320     return 0x0;
1321   }
1322   // primary vertex to be used by this candidate
1323   AliAODVertex *primVertexAOD  = PrimaryVertex(twoTrackArray,event);
1324   if(!primVertexAOD) return 0x0;
1325
1326
1327   Double_t d0z0[2],covd0z0[3];
1328   postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1329   d0[0] = d0z0[0];
1330   d0err[0] = TMath::Sqrt(covd0z0[0]);
1331   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1332   d0[1] = d0z0[0];
1333   d0err[1] = TMath::Sqrt(covd0z0[0]);
1334
1335   // create the object AliAODRecoDecayHF2Prong
1336   AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
1337   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
1338   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1339   the2Prong->SetProngIDs(2,id);
1340   delete primVertexAOD; primVertexAOD=NULL;
1341
1342
1343   if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
1344     // select D0->Kpi
1345     //Int_t checkD0,checkD0bar;
1346     //if(fD0toKpi)   okD0 = the2Prong->SelectD0(fD0toKpiCuts,checkD0,checkD0bar);
1347     if(fD0toKpi)   okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1348     //if(fDebug && fD0toKpi) printf("   %d\n",(Int_t)okD0);
1349     // select J/psi from B
1350     //Int_t checkJPSI;
1351     //if(fJPSItoEle) okJPSI        = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI);
1352     if(fJPSItoEle)   okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1353     //if(fDebug && fJPSItoEle) printf("   %d\n",(Int_t)okJPSI);
1354     // select D0->Kpi from Dstar
1355     //if(fDstar)     okD0fromDstar = the2Prong->SelectD0(fD0fromDstarCuts,checkD0,checkD0bar);
1356     if(fDstar)   okD0fromDstar = (Bool_t)fCutsD0fromDstar->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1357     //if(fDebug && fDstar) printf("   %d\n",(Int_t)okD0fromDstar);
1358   }
1359
1360
1361   // remove the primary vertex (was used only for selection)
1362   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1363     the2Prong->UnsetOwnPrimaryVtx();
1364   }
1365
1366   // get PID info from ESD
1367   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1368   if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1369   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1370   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1371   Double_t esdpid[10];
1372   for(Int_t i=0;i<5;i++) {
1373     esdpid[i]   = esdpid0[i];
1374     esdpid[5+i] = esdpid1[i];
1375   }
1376   the2Prong->SetPID(2,esdpid);
1377
1378   return the2Prong;
1379 }
1380 //----------------------------------------------------------------------------
1381 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
1382                              TObjArray *threeTrackArray,AliVEvent *event,
1383                              AliAODVertex *secVert,Double_t dispersion,
1384                              const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
1385                              Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1386                              Bool_t &ok3Prong) const
1387 {
1388   // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1389   // reconstruction cuts
1390   // E.Bruna, F.Prino
1391
1392
1393   ok3Prong=kFALSE;
1394   if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
1395
1396   Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1397   Double_t momentum[3];
1398
1399
1400   AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1401   AliESDtrack *negtrack  = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1402   AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1403
1404   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1405   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1406   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1407   postrack1->GetPxPyPz(momentum);
1408   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1409   negtrack->GetPxPyPz(momentum);
1410   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1411   postrack2->GetPxPyPz(momentum);
1412   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1413
1414   // invariant mass cut for D+, Ds, Lc
1415   Bool_t okMassCut=kFALSE;
1416   if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
1417   if(!okMassCut) {
1418     AliDebug(2," candidate didn't pass mass cut");
1419     return 0x0;
1420   }
1421
1422   // primary vertex to be used by this candidate
1423   AliAODVertex *primVertexAOD  = PrimaryVertex(threeTrackArray,event);
1424   if(!primVertexAOD) return 0x0;
1425
1426   Double_t d0z0[2],covd0z0[3];
1427   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1428   d0[0]=d0z0[0];
1429   d0err[0] = TMath::Sqrt(covd0z0[0]);
1430   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1431   d0[1]=d0z0[0];
1432   d0err[1] = TMath::Sqrt(covd0z0[0]);
1433   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1434   d0[2]=d0z0[0];
1435   d0err[2] = TMath::Sqrt(covd0z0[0]);
1436
1437
1438   // create the object AliAODRecoDecayHF3Prong
1439   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1440   Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1441   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]));
1442   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]));
1443   Short_t charge=(Short_t)(postrack1->Charge()*postrack2->Charge()*negtrack->Charge());
1444   AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1445   the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1446   UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1447   the3Prong->SetProngIDs(3,id);
1448
1449   delete primVertexAOD; primVertexAOD=NULL;
1450
1451   // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1452   if(f3Prong) {
1453     ok3Prong = kFALSE;
1454     //Int_t ok1,ok2;
1455     //Int_t dum1,dum2;
1456     //if(the3Prong->SelectDplus(fDplusCuts))   ok3Prong = kTRUE;
1457     //if(the3Prong->SelectDs(fDsCuts,ok1,ok2,dum1,dum2)) ok3Prong = kTRUE;
1458     //if(the3Prong->SelectLc(fLcCuts,ok1,ok2)) ok3Prong = kTRUE;
1459
1460     if(fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1461     if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1462     if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1463
1464   }
1465   //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1466
1467   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1468     the3Prong->UnsetOwnPrimaryVtx();
1469   }
1470
1471   // get PID info from ESD
1472   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1473   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1474   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1475   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1476   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1477   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1478
1479   Double_t esdpid[15];
1480   for(Int_t i=0;i<5;i++) {
1481     esdpid[i]    = esdpid0[i];
1482     esdpid[5+i]  = esdpid1[i];
1483     esdpid[10+i] = esdpid2[i];
1484   }
1485   the3Prong->SetPID(3,esdpid);
1486
1487   return the3Prong;
1488 }
1489 //----------------------------------------------------------------------------
1490 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1491                              TObjArray *fourTrackArray,AliVEvent *event,
1492                              AliAODVertex *secVert,
1493                              const AliAODVertex *vertexp1n1,
1494                              const AliAODVertex *vertexp1n1p2,
1495                              Double_t dcap1n1,Double_t dcap1n2,
1496                              Double_t dcap2n1,Double_t dcap2n2,
1497                              Bool_t &ok4Prong) const
1498 {
1499   // Make 4Prong candidates and check if they pass D0toKpipipi
1500   // reconstruction cuts
1501   // G.E.Bruno, R.Romita
1502
1503   ok4Prong=kFALSE;
1504   if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
1505
1506   Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1507
1508   AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1509   AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1510   AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1511   AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1512
1513   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1514   negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1515   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1516   negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1517
1518   Double_t momentum[3];
1519   postrack1->GetPxPyPz(momentum);
1520   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1521   negtrack1->GetPxPyPz(momentum);
1522   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1523   postrack2->GetPxPyPz(momentum);
1524   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1525   negtrack2->GetPxPyPz(momentum);
1526   px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1527
1528   // invariant mass cut for rho or D0 (try to improve coding here..)
1529   Bool_t okMassCut=kFALSE;
1530   if(!okMassCut && TMath::Abs(fD0to4ProngsCuts[8])<1.e-13){      //no PID, to be implemented with PID
1531     if(SelectInvMass(4,4,px,py,pz)) okMassCut=kTRUE;
1532   }
1533   if(!okMassCut) {
1534     //if(fDebug) printf(" candidate didn't pass mass cut\n");
1535     //printf(" candidate didn't pass mass cut\n");
1536     return 0x0;
1537   }
1538
1539   // primary vertex to be used by this candidate
1540   AliAODVertex *primVertexAOD  = PrimaryVertex(fourTrackArray,event);
1541   if(!primVertexAOD) return 0x0;
1542
1543   Double_t d0z0[2],covd0z0[3];
1544   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1545   d0[0]=d0z0[0];
1546   d0err[0] = TMath::Sqrt(covd0z0[0]);
1547   negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1548   d0[1]=d0z0[0];
1549   d0err[1] = TMath::Sqrt(covd0z0[0]);
1550   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1551   d0[2]=d0z0[0];
1552   d0err[2] = TMath::Sqrt(covd0z0[0]);
1553   negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1554   d0[3]=d0z0[0];
1555   d0err[3] = TMath::Sqrt(covd0z0[0]);
1556
1557
1558   // create the object AliAODRecoDecayHF4Prong
1559   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1560   Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
1561   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]));
1562   Double_t dist3=TMath::Sqrt((vertexp1n1p2->GetX()-pos[0])*(vertexp1n1p2->GetX()-pos[0])+(vertexp1n1p2->GetY()-pos[1])*(vertexp1n1p2->GetY()-pos[1])+(vertexp1n1p2->GetZ()-pos[2])*(vertexp1n1p2->GetZ()-pos[2]));
1563   Double_t dist4=TMath::Sqrt((secVert->GetX()-pos[0])*(secVert->GetX()-pos[0])+(secVert->GetY()-pos[1])*(secVert->GetY()-pos[1])+(secVert->GetZ()-pos[2])*(secVert->GetZ()-pos[2]));
1564   Short_t charge=0;
1565   AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
1566   the4Prong->SetOwnPrimaryVtx(primVertexAOD);
1567   UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
1568   the4Prong->SetProngIDs(4,id);
1569
1570   delete primVertexAOD; primVertexAOD=NULL;
1571
1572   //Int_t checkD0,checkD0bar;
1573   //ok4Prong=the4Prong->SelectD0(fD0to4ProngsCuts,checkD0,checkD0bar);
1574   ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
1575
1576
1577   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1578     the4Prong->UnsetOwnPrimaryVtx();
1579   }
1580
1581
1582   // get PID info from ESD
1583   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1584   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1585   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1586   if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
1587   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1588   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1589   Double_t esdpid3[5]={0.,0.,0.,0.,0.};
1590   if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
1591
1592   Double_t esdpid[20];
1593   for(Int_t i=0;i<5;i++) {
1594     esdpid[i]    = esdpid0[i];
1595     esdpid[5+i]  = esdpid1[i];
1596     esdpid[10+i] = esdpid2[i];
1597     esdpid[15+i] = esdpid3[i];
1598   }
1599   the4Prong->SetPID(4,esdpid);
1600
1601   return the4Prong;
1602 }
1603 //-----------------------------------------------------------------------------
1604 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
1605                                                     AliVEvent *event) const
1606 {
1607   // Returns primary vertex to be used for this candidate
1608
1609   AliESDVertex *vertexESD = 0;
1610   AliAODVertex *vertexAOD = 0;
1611
1612
1613   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1614     // primary vertex from the input event
1615
1616     vertexESD = new AliESDVertex(*fV1);
1617
1618   } else {
1619     // primary vertex specific to this candidate
1620
1621     Int_t nTrks = trkArray->GetEntriesFast();
1622     AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1623
1624     if(fRecoPrimVtxSkippingTrks) {
1625       // recalculating the vertex
1626
1627       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1628         Float_t diamondcovxy[3];
1629         event->GetDiamondCovXY(diamondcovxy);
1630         Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1631         Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1632         AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1633         vertexer->SetVtxStart(diamond);
1634         delete diamond; diamond=NULL;
1635         if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1636           vertexer->SetOnlyFitter();
1637       }
1638       Int_t skipped[1000];
1639       Int_t nTrksToSkip=0,id;
1640       AliExternalTrackParam *t = 0;
1641       for(Int_t i=0; i<nTrks; i++) {
1642         t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1643         id = (Int_t)t->GetID();
1644         if(id<0) continue;
1645         skipped[nTrksToSkip++] = id;
1646       }
1647       // TEMPORARY FIX
1648       // For AOD, skip also tracks without covariance matrix
1649       if(fInputAOD) {
1650         Double_t covtest[21];
1651         for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1652           AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1653           if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1654             id = (Int_t)vtrack->GetID();
1655             if(id<0) continue;
1656             skipped[nTrksToSkip++] = id;
1657           }
1658         }
1659       }
1660       //
1661       vertexer->SetSkipTracks(nTrksToSkip,skipped);
1662       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1663
1664     } else if(fRmTrksFromPrimVtx) {
1665       // removing the prongs tracks
1666
1667       TObjArray rmArray(nTrks);
1668       UShort_t *rmId = new UShort_t[nTrks];
1669       AliESDtrack *esdTrack = 0;
1670       AliESDtrack *t = 0;
1671       for(Int_t i=0; i<nTrks; i++) {
1672         t = (AliESDtrack*)trkArray->UncheckedAt(i);
1673         esdTrack = new AliESDtrack(*t);
1674         rmArray.AddLast(esdTrack);
1675         if(esdTrack->GetID()>=0) {
1676           rmId[i]=(UShort_t)esdTrack->GetID();
1677         } else {
1678           rmId[i]=9999;
1679         }
1680       }
1681       Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
1682       vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1683       delete [] rmId; rmId=NULL;
1684       rmArray.Delete();
1685
1686     }
1687
1688     if(!vertexESD) return vertexAOD;
1689     if(vertexESD->GetNContributors()<=0) {
1690       AliDebug(2,"vertexing failed");
1691       delete vertexESD; vertexESD=NULL;
1692       return vertexAOD;
1693     }
1694
1695     delete vertexer; vertexer=NULL;
1696
1697   }
1698
1699   // convert to AliAODVertex
1700   Double_t pos[3],cov[6],chi2perNDF;
1701   vertexESD->GetXYZ(pos); // position
1702   vertexESD->GetCovMatrix(cov); //covariance matrix
1703   chi2perNDF = vertexESD->GetChi2toNDF();
1704   delete vertexESD; vertexESD=NULL;
1705
1706   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1707
1708   return vertexAOD;
1709 }
1710 //-----------------------------------------------------------------------------
1711 void AliAnalysisVertexingHF::PrintStatus() const {
1712   // Print parameters being used
1713
1714   //printf("Preselections:\n");
1715   //   fTrackFilter->Dump();
1716   if(fSecVtxWithKF) {
1717     printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1718   } else {
1719     printf("Secondary vertex with AliVertexerTracks\n");
1720   }
1721   if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1722   if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1723   if(fD0toKpi) {
1724     printf("Reconstruct D0->Kpi candidates with cuts:\n");
1725     if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
1726     /*
1727     printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
1728     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
1729     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
1730     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
1731     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
1732     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
1733     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
1734     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
1735     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
1736     */
1737   }
1738   if(fDstar) {
1739     if(fFindVertexForDstar) {
1740       printf("    Reconstruct a secondary vertex for the D*\n");
1741     } else {
1742       printf("    Assume the D* comes from the primary vertex\n");
1743     }
1744     printf("    |M-MD*| [GeV]    < %f\n",fDstarCuts[0]);
1745     printf("    |M_Kpipi-M_Kpi-(MD*-MD0)| [GeV]  < %f\n",fDstarCuts[1]);
1746     printf("    pTpisoft [GeV/c]    > %f\n",fDstarCuts[2]);
1747     printf("    pTpisoft [GeV/c]    < %f\n",fDstarCuts[3]);
1748     printf("    Theta(pisoft,D0plane) < %f\n",fDstarCuts[4]);
1749     printf("Reconstruct D*->D0pi candidates with cuts:\n");
1750     printf("   D0 from D* cuts:\n");
1751     if(fCutsD0fromDstar) fCutsD0fromDstar->PrintAll();
1752     /*
1753     printf("    |M-MD0| [GeV]    < %f\n",fD0fromDstarCuts[0]);
1754     printf("    dca    [cm]  < %f\n",fD0fromDstarCuts[1]);
1755     printf("    cosThetaStar     < %f\n",fD0fromDstarCuts[2]);
1756     printf("    pTK     [GeV/c]    > %f\n",fD0fromDstarCuts[3]);
1757     printf("    pTpi    [GeV/c]    > %f\n",fD0fromDstarCuts[4]);
1758     printf("    |d0K|  [cm]  < %f\n",fD0fromDstarCuts[5]);
1759     printf("    |d0pi| [cm]  < %f\n",fD0fromDstarCuts[6]);
1760     printf("    d0d0  [cm^2] < %f\n",fD0fromDstarCuts[7]);
1761     printf("    cosThetaPoint    > %f\n",fD0fromDstarCuts[8]);
1762     */
1763   }
1764   if(fJPSItoEle) {
1765     printf("Reconstruct J/psi from B candidates with cuts:\n");
1766     if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
1767     /*
1768     printf("    |M-MJPSI| [GeV]    < %f\n",fBtoJPSICuts[0]);
1769     printf("    dca    [cm]  < %f\n",fBtoJPSICuts[1]);
1770     printf("    cosThetaStar     < %f\n",fBtoJPSICuts[2]);
1771     printf("    pTP     [GeV/c]    > %f\n",fBtoJPSICuts[3]);
1772     printf("    pTN    [GeV/c]    > %f\n",fBtoJPSICuts[4]);
1773     printf("    |d0P|  [cm]  < %f\n",fBtoJPSICuts[5]);
1774     printf("    |d0N| [cm]  < %f\n",fBtoJPSICuts[6]);
1775     printf("    d0d0  [cm^2] < %f\n",fBtoJPSICuts[7]);
1776     printf("    cosThetaPoint    > %f\n",fBtoJPSICuts[8]);
1777     */
1778   }
1779   if(f3Prong) {
1780     printf("Reconstruct 3 prong candidates.\n");
1781     printf("  D+->Kpipi cuts:\n");
1782     if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
1783     /*
1784     printf("    |M-MD+| [GeV]    < %f\n",fDplusCuts[0]);
1785     printf("    pTK     [GeV/c]    > %f\n",fDplusCuts[1]);
1786     printf("    pTPi    [GeV/c]    > %f\n",fDplusCuts[2]);
1787     printf("    |d0K|  [cm]  > %f\n",fDplusCuts[3]);
1788     printf("    |d0Pi| [cm]  > %f\n",fDplusCuts[4]);
1789     printf("    dist12    [cm]  < %f\n",fDplusCuts[5]);
1790     printf("    sigmavert [cm]   < %f\n",fDplusCuts[6]);
1791     printf("    dist prim-sec [cm] > %f\n",fDplusCuts[7]);
1792     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDplusCuts[8]);
1793     printf("    cosThetaPoint    > %f\n",fDplusCuts[9]);
1794     printf("    Sum d0^2 [cm^2]  > %f\n",fDplusCuts[10]);
1795     printf("    dca cut [cm]  < %f\n",fDplusCuts[11]);
1796     */
1797     printf("  Ds->KKpi cuts:\n");
1798     if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
1799     /*
1800     printf("    |M-MDs| [GeV]    < %f\n",fDsCuts[0]);
1801     printf("    pTK     [GeV/c]    > %f\n",fDsCuts[1]);
1802     printf("    pTPi    [GeV/c]    > %f\n",fDsCuts[2]);
1803     printf("    |d0K|  [cm]  > %f\n",fDsCuts[3]);
1804     printf("    |d0Pi| [cm]  > %f\n",fDsCuts[4]);
1805     printf("    dist12    [cm]  < %f\n",fDsCuts[5]);
1806     printf("    sigmavert [cm]   < %f\n",fDsCuts[6]);
1807     printf("    dist prim-sec [cm] > %f\n",fDsCuts[7]);
1808     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fDsCuts[8]);
1809     printf("    cosThetaPoint    > %f\n",fDsCuts[9]);
1810     printf("    Sum d0^2 [cm^2]  > %f\n",fDsCuts[10]);
1811     printf("    dca cut [cm]  < %f\n",fDsCuts[11]);
1812     printf("    Inv. Mass  phi [GeV]  < %f\n",fDsCuts[12]);
1813     printf("    Inv. Mass  K0* [GeV]  < %f\n",fDsCuts[13]);
1814     */
1815     printf("  Lc->pKpi cuts:\n");
1816     if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
1817     /*
1818     printf("    |M-MLc| [GeV]    < %f\n",fLcCuts[0]);
1819     printf("    pTP     [GeV/c]    > %f\n",fLcCuts[1]);
1820     printf("    pTPi and pTK [GeV/c]    > %f\n",fLcCuts[2]);
1821     printf("    |d0P|  [cm]  > %f\n",fLcCuts[3]);
1822     printf("    |d0Pi| and |d0K| [cm]  > %f\n",fLcCuts[4]);
1823     printf("    dist12    [cm]  < %f\n",fLcCuts[5]);
1824     printf("    sigmavert [cm]   < %f\n",fLcCuts[6]);
1825     printf("    dist prim-sec [cm] > %f\n",fLcCuts[7]);
1826     printf("    pM=Max{pT1,pT2,pT3} [GeV/c] > %f\n",fLcCuts[8]);
1827     printf("    cosThetaPoint    > %f\n",fLcCuts[9]);
1828     printf("    Sum d0^2 [cm^2]  > %f\n",fLcCuts[10]);
1829     printf("    dca cut [cm]  < %f\n",fLcCuts[11]);
1830     */
1831   }
1832   if(f4Prong) {
1833     printf("Reconstruct 4 prong candidates.\n");
1834     printf("  D0->Kpipipi cuts:\n");
1835     if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
1836   }
1837   if(fCascades) {
1838     printf("Reconstruct cascades candidates formed with v0s.\n");
1839     printf("  Lc -> k0s P & Lc -> L Pi cuts:\n");
1840     if(fCutsLctoV0) fCutsLctoV0->PrintAll();
1841   }
1842
1843   return;
1844 }
1845 //-----------------------------------------------------------------------------
1846 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
1847                                                                  Double_t &dispersion,Bool_t useTRefArray) const
1848 {
1849   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1850
1851   AliESDVertex *vertexESD = 0;
1852   AliAODVertex *vertexAOD = 0;
1853
1854   if(!fSecVtxWithKF) { // AliVertexerTracks
1855
1856     AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1857     vertexer->SetVtxStart(fV1);
1858     vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1859     delete vertexer; vertexer=NULL;
1860
1861     if(!vertexESD) return vertexAOD;
1862
1863     if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
1864       AliDebug(2,"vertexing failed");
1865       delete vertexESD; vertexESD=NULL;
1866       return vertexAOD;
1867     }
1868
1869   } else { // Kalman Filter vertexer (AliKFParticle)
1870
1871     AliKFParticle::SetField(fBzkG);
1872
1873     AliKFVertex vertexKF;
1874
1875     Int_t nTrks = trkArray->GetEntriesFast();
1876     for(Int_t i=0; i<nTrks; i++) {
1877       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1878       AliKFParticle daughterKF(*esdTrack,211);
1879       vertexKF.AddDaughter(daughterKF);
1880     }
1881     vertexESD = new AliESDVertex(vertexKF.Parameters(),
1882                                  vertexKF.CovarianceMatrix(),
1883                                  vertexKF.GetChi2(),
1884                                  vertexKF.GetNContributors());
1885
1886   }
1887
1888   // convert to AliAODVertex
1889   Double_t pos[3],cov[6],chi2perNDF;
1890   vertexESD->GetXYZ(pos); // position
1891   vertexESD->GetCovMatrix(cov); //covariance matrix
1892   chi2perNDF = vertexESD->GetChi2toNDF();
1893   dispersion = vertexESD->GetDispersion();
1894   delete vertexESD; vertexESD=NULL;
1895
1896   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
1897   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
1898
1899   return vertexAOD;
1900 }
1901 //-----------------------------------------------------------------------------
1902 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1903                                              Int_t nprongs,
1904                                              Double_t *px,
1905                                              Double_t *py,
1906                                              Double_t *pz) const {
1907   // Check invariant mass cut
1908
1909   Short_t dummycharge=0;
1910   Double_t *dummyd0 = new Double_t[nprongs];
1911   for(Int_t ip=0;ip<nprongs;ip++) dummyd0[ip]=0.;
1912   AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
1913   delete [] dummyd0; dummyd0=NULL;
1914
1915   UInt_t pdg2[2],pdg3[3],pdg4[4];
1916   Double_t mPDG,minv;
1917
1918   Bool_t retval=kFALSE;
1919   switch (decay)
1920     {
1921     case 0:                  // D0->Kpi
1922       pdg2[0]=211; pdg2[1]=321;
1923       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1924       minv = rd->InvMass(nprongs,pdg2);
1925       //if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1926       if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
1927       pdg2[0]=321; pdg2[1]=211;
1928       minv = rd->InvMass(nprongs,pdg2);
1929       //if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1930       if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
1931       break;
1932     case 1:                  // JPSI->ee
1933       pdg2[0]=11; pdg2[1]=11;
1934       mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
1935       minv = rd->InvMass(nprongs,pdg2);
1936       //if(TMath::Abs(minv-mPDG)<fBtoJPSICuts[0]) retval=kTRUE;
1937       if(TMath::Abs(minv-mPDG)<fCutsJpsitoee->GetMassCut()) retval=kTRUE;
1938       break;
1939     case 2:                  // D+->Kpipi
1940       pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
1941       mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
1942       minv = rd->InvMass(nprongs,pdg3);
1943       //if(TMath::Abs(minv-mPDG)<fDplusCuts[0]) retval=kTRUE;
1944       if(TMath::Abs(minv-mPDG)<fCutsDplustoKpipi->GetMassCut()) retval=kTRUE;
1945                             // Ds+->KKpi
1946       pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
1947       mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1948       minv = rd->InvMass(nprongs,pdg3);
1949       //if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1950       if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
1951       pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
1952       minv = rd->InvMass(nprongs,pdg3);
1953       //if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1954       if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
1955                             // Lc->pKpi
1956       pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
1957       mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1958       minv = rd->InvMass(nprongs,pdg3);
1959       //if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1960       if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
1961       pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
1962       minv = rd->InvMass(nprongs,pdg3);
1963       //if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1964       if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
1965       break;
1966     case 3:                  // D*->D0pi
1967       pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
1968       mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
1969       minv = rd->InvMass(nprongs,pdg2);
1970       if(TMath::Abs(minv-mPDG)<fDstarCuts[0]) retval=kTRUE;
1971       break;
1972     case 4:                 // D0->Kpipipi without PID
1973       pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
1974       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1975       minv = rd->InvMass(nprongs,pdg4);
1976       //if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1977       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1978       pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
1979       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1980       minv = rd->InvMass(nprongs,pdg4);
1981       //if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1982       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1983       pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
1984       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1985       minv = rd->InvMass(nprongs,pdg4);
1986       //if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1987       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1988       pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
1989       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1990       minv = rd->InvMass(nprongs,pdg4);
1991       //if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1992       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1993       break;
1994     default:
1995       printf("SelectInvMass(): wrong decay selection\n");
1996       break;
1997     }
1998
1999   delete rd; rd=NULL;
2000
2001   return retval;
2002 }
2003 //-----------------------------------------------------------------------------
2004 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
2005                                    TObjArray &seleTrksArray,Int_t &nSeleTrks,
2006                                    UChar_t *seleFlags,Int_t *evtNumber)
2007 {
2008   // Apply single-track preselection.
2009   // Fill a TObjArray with selected tracks (for displaced vertices or
2010   // soft pion from D*). Selection flag stored in seleFlags.
2011   // Create the AliESDVertex object (convert from AliAODVertex if necessary)
2012   // In case of AOD input, also fill fAODMap for track index<->ID
2013
2014   const AliVVertex *vprimary = event->GetPrimaryVertex();
2015
2016   if(fV1) { delete fV1; fV1=NULL; }
2017   if(fAODMap) { delete fAODMap; fAODMap=NULL; }
2018
2019   Int_t nindices=0;
2020   UShort_t *indices = 0;
2021   Double_t pos[3],cov[6];
2022
2023   if(!fInputAOD) { // ESD
2024     fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
2025   } else {         // AOD
2026     vprimary->GetXYZ(pos);
2027     vprimary->GetCovarianceMatrix(cov);
2028     fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2029     indices = new UShort_t[event->GetNumberOfTracks()];
2030     fAODMapSize = 100000;
2031     fAODMap = new Int_t[fAODMapSize];
2032   }
2033
2034
2035   Int_t entries = (Int_t)event->GetNumberOfTracks();
2036   Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
2037   nSeleTrks=0;
2038
2039   // transfer ITS tracks from event to arrays
2040   for(Int_t i=0; i<entries; i++) {
2041     AliVTrack *track;
2042     track = (AliVTrack*)event->GetTrack(i);
2043
2044     // skip pure ITS SA tracks
2045     if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2046
2047     // TEMPORARY: check that the cov matrix is there
2048     Double_t covtest[21];
2049     if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2050     //
2051
2052     if(fInputAOD) {
2053       AliAODTrack *aodt = (AliAODTrack*)track;
2054       if(aodt->GetUsedForPrimVtxFit()) {
2055         indices[nindices]=aodt->GetID(); nindices++;
2056       }
2057       fAODMap[(Int_t)aodt->GetID()] = i;
2058     }
2059
2060     AliESDtrack *esdt = 0;
2061
2062     if(!fInputAOD) {
2063       esdt = (AliESDtrack*)track;
2064     } else {
2065       esdt = new AliESDtrack(track);
2066     }
2067
2068     // single track selection
2069     okDisplaced=kFALSE; okSoftPi=kFALSE;
2070     if(fMixEvent){
2071       evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
2072       const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
2073       Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
2074       eventVtx->GetXYZ(vtxPos);
2075       vprimary->GetXYZ(primPos);
2076       eventVtx->GetCovarianceMatrix(primCov);
2077       for(Int_t ind=0;ind<3;ind++){
2078         trasl[ind]=vtxPos[ind]-primPos[ind];
2079       }
2080
2081       Bool_t isTransl=esdt->Translate(trasl,primCov);
2082       if(!isTransl) {
2083         delete esdt;
2084         esdt = NULL;
2085         continue;
2086       }
2087     }
2088
2089     if(SingleTrkCuts(esdt,okDisplaced,okSoftPi)) {
2090       seleTrksArray.AddLast(esdt);
2091       seleFlags[nSeleTrks]=0;
2092       if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2093       if(okSoftPi)    SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2094       nSeleTrks++;
2095     } else {
2096       if(fInputAOD) delete esdt;
2097       esdt = NULL;
2098       continue;
2099     }
2100
2101   } // end loop on tracks
2102
2103   // primary vertex from AOD
2104   if(fInputAOD) {
2105     delete fV1; fV1=NULL;
2106     vprimary->GetXYZ(pos);
2107     vprimary->GetCovarianceMatrix(cov);
2108     Double_t chi2toNDF = vprimary->GetChi2perNDF();
2109     Int_t ncontr=nindices;
2110     if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
2111     Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
2112     fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2113     fV1->SetTitle(vprimary->GetTitle());
2114     fV1->SetIndices(nindices,indices);
2115   }
2116   if(indices) { delete [] indices; indices=NULL; }
2117
2118
2119   return;
2120 }
2121 //-----------------------------------------------------------------------------
2122 void AliAnalysisVertexingHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
2123                                    Double_t cut2,Double_t cut3,Double_t cut4,
2124                                    Double_t cut5,Double_t cut6,
2125                                    Double_t cut7,Double_t cut8)
2126 {
2127   // Set the cuts for D0 selection
2128   fD0toKpiCuts[0] = cut0;
2129   fD0toKpiCuts[1] = cut1;
2130   fD0toKpiCuts[2] = cut2;
2131   fD0toKpiCuts[3] = cut3;
2132   fD0toKpiCuts[4] = cut4;
2133   fD0toKpiCuts[5] = cut5;
2134   fD0toKpiCuts[6] = cut6;
2135   fD0toKpiCuts[7] = cut7;
2136   fD0toKpiCuts[8] = cut8;
2137
2138   return;
2139 }
2140 //-----------------------------------------------------------------------------
2141 void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9])
2142 {
2143   // Set the cuts for D0 selection
2144
2145   for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
2146
2147   return;
2148 }
2149 //-----------------------------------------------------------------------------
2150 void AliAnalysisVertexingHF::SetD0fromDstarCuts(Double_t cut0,Double_t cut1,
2151                                    Double_t cut2,Double_t cut3,Double_t cut4,
2152                                    Double_t cut5,Double_t cut6,
2153                                    Double_t cut7,Double_t cut8)
2154 {
2155   // Set the cuts for D0 from D* selection
2156   fD0fromDstarCuts[0] = cut0;
2157   fD0fromDstarCuts[1] = cut1;
2158   fD0fromDstarCuts[2] = cut2;
2159   fD0fromDstarCuts[3] = cut3;
2160   fD0fromDstarCuts[4] = cut4;
2161   fD0fromDstarCuts[5] = cut5;
2162   fD0fromDstarCuts[6] = cut6;
2163   fD0fromDstarCuts[7] = cut7;
2164   fD0fromDstarCuts[8] = cut8;
2165
2166   return;
2167 }
2168 //-----------------------------------------------------------------------------
2169 void AliAnalysisVertexingHF::SetD0fromDstarCuts(const Double_t cuts[9])
2170 {
2171   // Set the cuts for D0 from D* selection
2172
2173   for(Int_t i=0; i<9; i++) fD0fromDstarCuts[i] = cuts[i];
2174
2175   return;
2176 }
2177 //-----------------------------------------------------------------------------
2178 void AliAnalysisVertexingHF::SetDstarCuts(Double_t cut0,Double_t cut1,
2179                                           Double_t cut2,Double_t cut3,
2180                                           Double_t cut4)
2181 {
2182   // Set the cuts for D* selection
2183   fDstarCuts[0] = cut0;
2184   fDstarCuts[1] = cut1;
2185   fDstarCuts[2] = cut2;
2186   fDstarCuts[3] = cut3;
2187   fDstarCuts[4] = cut4;
2188
2189   return;
2190 }
2191 //-----------------------------------------------------------------------------
2192 void AliAnalysisVertexingHF::SetDstarCuts(const Double_t cuts[5])
2193 {
2194   // Set the cuts for D* selection
2195
2196   for(Int_t i=0; i<5; i++) fDstarCuts[i] = cuts[i];
2197
2198   return;
2199 }
2200 //-----------------------------------------------------------------------------
2201 void AliAnalysisVertexingHF::SetBtoJPSICuts(Double_t cut0,Double_t cut1,
2202                                    Double_t cut2,Double_t cut3,Double_t cut4,
2203                                    Double_t cut5,Double_t cut6,
2204                                    Double_t cut7,Double_t cut8)
2205 {
2206   // Set the cuts for J/psi from B selection
2207   fBtoJPSICuts[0] = cut0;
2208   fBtoJPSICuts[1] = cut1;
2209   fBtoJPSICuts[2] = cut2;
2210   fBtoJPSICuts[3] = cut3;
2211   fBtoJPSICuts[4] = cut4;
2212   fBtoJPSICuts[5] = cut5;
2213   fBtoJPSICuts[6] = cut6;
2214   fBtoJPSICuts[7] = cut7;
2215   fBtoJPSICuts[8] = cut8;
2216
2217   return;
2218 }
2219 //-----------------------------------------------------------------------------
2220 void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9])
2221 {
2222   // Set the cuts for J/psi from B selection
2223
2224   for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i];
2225
2226   return;
2227 }
2228 //-----------------------------------------------------------------------------
2229 void AliAnalysisVertexingHF::SetDplusCuts(Double_t cut0,Double_t cut1,
2230                                    Double_t cut2,Double_t cut3,Double_t cut4,
2231                                    Double_t cut5,Double_t cut6,
2232                                    Double_t cut7,Double_t cut8,
2233                                    Double_t cut9,Double_t cut10,Double_t cut11)
2234 {
2235   // Set the cuts for Dplus->Kpipi selection
2236   fDplusCuts[0] = cut0;
2237   fDplusCuts[1] = cut1;
2238   fDplusCuts[2] = cut2;
2239   fDplusCuts[3] = cut3;
2240   fDplusCuts[4] = cut4;
2241   fDplusCuts[5] = cut5;
2242   fDplusCuts[6] = cut6;
2243   fDplusCuts[7] = cut7;
2244   fDplusCuts[8] = cut8;
2245   fDplusCuts[9] = cut9;
2246   fDplusCuts[10] = cut10;
2247   fDplusCuts[11] = cut11;
2248
2249   return;
2250 }
2251 //-----------------------------------------------------------------------------
2252 void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12])
2253 {
2254   // Set the cuts for Dplus->Kpipi selection
2255
2256   for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i];
2257
2258   return;
2259 }
2260 //-----------------------------------------------------------------------------
2261 void AliAnalysisVertexingHF::SetDsCuts(Double_t cut0,Double_t cut1,
2262                                        Double_t cut2,Double_t cut3,
2263                                        Double_t cut4,Double_t cut5,
2264                                        Double_t cut6,Double_t cut7,
2265                                        Double_t cut8,Double_t cut9,
2266                                        Double_t cut10,Double_t cut11,
2267                                        Double_t cut12,Double_t cut13)
2268 {
2269   // Set the cuts for Ds->KKpi selection
2270   fDsCuts[0] = cut0;
2271   fDsCuts[1] = cut1;
2272   fDsCuts[2] = cut2;
2273   fDsCuts[3] = cut3;
2274   fDsCuts[4] = cut4;
2275   fDsCuts[5] = cut5;
2276   fDsCuts[6] = cut6;
2277   fDsCuts[7] = cut7;
2278   fDsCuts[8] = cut8;
2279   fDsCuts[9] = cut9;
2280   fDsCuts[10] = cut10;
2281   fDsCuts[11] = cut11;
2282   fDsCuts[12] = cut12;
2283   fDsCuts[13] = cut13;
2284
2285   return;
2286 }
2287 //-----------------------------------------------------------------------------
2288 void AliAnalysisVertexingHF::SetDsCuts(const Double_t cuts[13])
2289 {
2290   // Set the cuts for Ds->KKpi selection
2291
2292   for(Int_t i=0; i<14; i++) fDsCuts[i] = cuts[i];
2293
2294   return;
2295 }
2296 //-----------------------------------------------------------------------------
2297 void AliAnalysisVertexingHF::SetLcCuts(Double_t cut0,Double_t cut1,
2298                                    Double_t cut2,Double_t cut3,Double_t cut4,
2299                                    Double_t cut5,Double_t cut6,
2300                                    Double_t cut7,Double_t cut8,
2301                                    Double_t cut9,Double_t cut10,Double_t cut11)
2302 {
2303   // Set the cuts for Lc->pKpi selection
2304   fLcCuts[0] = cut0;
2305   fLcCuts[1] = cut1;
2306   fLcCuts[2] = cut2;
2307   fLcCuts[3] = cut3;
2308   fLcCuts[4] = cut4;
2309   fLcCuts[5] = cut5;
2310   fLcCuts[6] = cut6;
2311   fLcCuts[7] = cut7;
2312   fLcCuts[8] = cut8;
2313   fLcCuts[9] = cut9;
2314   fLcCuts[10] = cut10;
2315   fLcCuts[11] = cut11;
2316
2317   return;
2318 }
2319 //-----------------------------------------------------------------------------
2320 void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[12])
2321 {
2322   // Set the cuts for Lc->pKpi selection
2323
2324   for(Int_t i=0; i<12; i++) fLcCuts[i] = cuts[i];
2325
2326   return;
2327 }
2328
2329 //-----------------------------------------------------------------------------
2330 void AliAnalysisVertexingHF::SetLctoV0Cuts(Double_t cut0,Double_t cut1,
2331                                            Double_t cut2,Double_t cut3,Double_t cut4,
2332                                            Double_t cut5,Double_t cut6,
2333                                            Double_t cut7,Double_t cut8)
2334 {
2335   // Set the cuts for Lc->V0+bachelor selection
2336   fLctoV0Cuts[0] = cut0;
2337   fLctoV0Cuts[1] = cut1;
2338   fLctoV0Cuts[2] = cut2;
2339   fLctoV0Cuts[3] = cut3;
2340   fLctoV0Cuts[4] = cut4;
2341   fLctoV0Cuts[5] = cut5;
2342   fLctoV0Cuts[6] = cut6;
2343   fLctoV0Cuts[7] = cut7;
2344   fLctoV0Cuts[8] = cut8;
2345
2346   return;
2347 }
2348
2349 //-----------------------------------------------------------------------------
2350 void AliAnalysisVertexingHF::SetLctoV0Cuts(const Double_t cuts[8])
2351 {
2352   // Set the cuts for Lc-> V0 + bachelor selection
2353
2354   for(Int_t i=0; i<8; i++) fLctoV0Cuts[i] = cuts[i];
2355
2356   return;
2357 }
2358
2359 //-----------------------------------------------------------------------------
2360 void AliAnalysisVertexingHF::SetD0to4ProngsCuts(Double_t cut0,Double_t cut1,
2361                                    Double_t cut2,Double_t cut3,Double_t cut4,
2362                                    Double_t cut5,Double_t cut6,
2363                                    Double_t cut7,Double_t cut8)
2364 {
2365   // Set the cuts for D0->Kpipipi selection
2366
2367   fD0to4ProngsCuts[0] = cut0;
2368   fD0to4ProngsCuts[1] = cut1;
2369   fD0to4ProngsCuts[2] = cut2;
2370   fD0to4ProngsCuts[3] = cut3;
2371   fD0to4ProngsCuts[4] = cut4;
2372   fD0to4ProngsCuts[5] = cut5;
2373   fD0to4ProngsCuts[6] = cut6;
2374   fD0to4ProngsCuts[7] = cut7;
2375   fD0to4ProngsCuts[8] = cut8;
2376
2377   return;
2378 }
2379 //-----------------------------------------------------------------------------
2380 void AliAnalysisVertexingHF::SetD0to4ProngsCuts(const Double_t cuts[9])
2381 {
2382   // Set the cuts for D0->Kpipipi selection
2383
2384   for(Int_t i=0; i<9; i++) fD0to4ProngsCuts[i] = cuts[i];
2385
2386   return;
2387 }
2388 //-----------------------------------------------------------------------------
2389 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2390                                              Bool_t &okDisplaced,Bool_t &okSoftPi) const
2391 {
2392   // Check if track passes some kinematical cuts
2393
2394   // this is needed to store the impact parameters
2395   trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2396
2397   UInt_t selectInfo;
2398   //
2399   // Track selection, displaced tracks
2400   selectInfo = 0;
2401   if(fTrackFilter) {
2402     selectInfo = fTrackFilter->IsSelected(trk);
2403   }
2404   if(selectInfo) okDisplaced=kTRUE;
2405   // Track selection, soft pions
2406   selectInfo = 0;
2407   if(fDstar && fTrackFilterSoftPi) {
2408     selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2409   }
2410   if(selectInfo) okSoftPi=kTRUE;
2411
2412   if(okDisplaced || okSoftPi) return kTRUE;
2413
2414   return kFALSE;
2415 }
2416
2417
2418 //-----------------------------------------------------------------------------
2419 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2420   //
2421   // Transform ESDv0 to AODv0
2422   //
2423   //  this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2424   //  and creates an AODv0 out of them
2425   //
2426   double vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
2427   AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2428
2429   // create the v0 neutral track to compute the DCA to the primary vertex
2430   Double_t xyz[3], pxpypz[3];
2431   esdV0->XvYvZv(xyz);
2432   esdV0->PxPyPz(pxpypz);
2433   Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2434   AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2435   if(!trackesdV0) return 0;
2436   Double_t d0z0[2],covd0z0[3];
2437   trackesdV0->PropagateToDCA(PrimaryVertex(),fBzkG,kVeryBig,d0z0,covd0z0);
2438   Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2439   // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2440   Double_t dcaV0DaughterToPrimVertex[2];
2441   AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2442   AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
2443   if( !posV0track || !negV0track) return 0;
2444   posV0track->PropagateToDCA(PrimaryVertex(),fBzkG,kVeryBig,d0z0,covd0z0);
2445   //  if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2446   //  else
2447   dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
2448   negV0track->PropagateToDCA(PrimaryVertex(),fBzkG,kVeryBig,d0z0,covd0z0);
2449   //  if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2450   //  else
2451   dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
2452   double dcaV0Daughters = esdV0->GetDcaV0Daughters();
2453   double pmom[3];  double nmom[3];
2454   esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2455   esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2456
2457   AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2458   aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2459
2460   if(trackesdV0) delete trackesdV0;
2461
2462   return aodV0;
2463 }
2464 //-----------------------------------------------------------------------------