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