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