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