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 = (AliAODTrack*)event->GetTrack(fAODMap[id]);
1503     v->AddDaughter(aodTrack);
1504   }
1505
1506   return;
1507 }
1508 //---------------------------------------------------------------------------
1509 void AliAnalysisVertexingHF::FixReferences(AliAODEvent *aod)
1510 {
1511   // Checks that the references to the daughter tracks are properly
1512   // assigned and reassigns them if needed
1513   //
1514   //AliCodeTimerAuto("",0);
1515
1516
1517   TClonesArray *inputArray=(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
1518   if(!inputArray) return;
1519
1520   AliAODTrack *track = 0;
1521   AliAODVertex *vertex = 0;
1522
1523   Bool_t needtofix=kFALSE;
1524   for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1525     vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1526     for(Int_t id=0; id<vertex->GetNDaughters(); id++) {
1527       track = (AliAODTrack*)vertex->GetDaughter(id);
1528       if(!track->GetStatus()) needtofix=kTRUE;
1529     }
1530     if(needtofix) break;
1531   }
1532
1533   if(!needtofix) return;
1534
1535
1536   printf("Fixing references\n");
1537
1538   fAODMapSize = 100000;
1539   fAODMap = new Int_t[fAODMapSize];
1540   memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
1541
1542   for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
1543     track = aod->GetTrack(i);
1544
1545     // skip pure ITS SA tracks
1546     if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
1547
1548     // skip tracks without ITS
1549     if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
1550
1551     // TEMPORARY: check that the cov matrix is there
1552     Double_t covtest[21];
1553     if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1554     //
1555
1556     Int_t ind = (Int_t)track->GetID();
1557     if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
1558   }
1559
1560
1561   Int_t ids[4]={-1,-1,-1,-1};
1562   for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1563     Bool_t cascade=kFALSE;
1564     vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1565     Int_t id=0;
1566     Int_t nDgs = vertex->GetNDaughters();
1567     for(id=0; id<nDgs; id++) {
1568       track = (AliAODTrack*)vertex->GetDaughter(id);
1569       if(track->Charge()==0) {cascade=kTRUE; continue;} // cascade
1570       ids[id]=(Int_t)track->GetID();
1571       vertex->RemoveDaughter(track);
1572     }
1573     if(cascade) continue;
1574     for(id=0; id<nDgs; id++) {
1575       if (ids[id]>-1 && ids[id] < fAODMapSize) {
1576         track = aod->GetTrack(fAODMap[ids[id]]);
1577         vertex->AddDaughter(track);
1578       }
1579     }
1580
1581   }
1582
1583   return;
1584 }
1585 //----------------------------------------------------------------------------
1586 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1587                                    TObjArray *twoTrackArray,AliVEvent *event,
1588                                    AliAODVertex *secVert,
1589                                    AliAODRecoDecayHF2Prong *rd2Prong,
1590                                    Double_t dca,
1591                                    Bool_t &okDstar)
1592 {
1593   // Make the cascade as a 2Prong decay and check if it passes Dstar
1594   // reconstruction cuts
1595   //AliCodeTimerAuto("",0);
1596
1597   okDstar = kFALSE;
1598
1599   Bool_t dummy1,dummy2,dummy3;
1600
1601   // We use Make2Prong to construct the AliAODRecoCascadeHF
1602   // (which inherits from AliAODRecoDecayHF2Prong)
1603   AliAODRecoCascadeHF *theCascade =
1604     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1605                                      dummy1,dummy2,dummy3);
1606   if(!theCascade) return 0x0;
1607
1608   // charge
1609   AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1610   theCascade->SetCharge(trackPi->Charge());
1611
1612   //--- selection cuts
1613   //
1614   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1615   if(fInputAOD){
1616     Int_t idSoftPi=(Int_t)trackPi->GetID();
1617     if (idSoftPi > -1 && idSoftPi < fAODMapSize) {
1618       AliAODTrack* trackPiAOD=(AliAODTrack*)event->GetTrack(fAODMap[idSoftPi]);
1619       tmpCascade->GetSecondaryVtx()->AddDaughter(trackPiAOD);
1620     }
1621   }else{
1622     tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
1623   }
1624   tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
1625
1626   AliAODVertex *primVertexAOD=0;
1627   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1628     // take event primary vertex
1629     primVertexAOD = PrimaryVertex();
1630     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1631     rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
1632   }
1633   // select D*->D0pi
1634   if(fDstar) {
1635     okDstar = (Bool_t)fCutsDStartoKpipi->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1636     if(okDstar) theCascade->SetSelectionBit(AliRDHFCuts::kDstarCuts);
1637   }
1638   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1639   tmpCascade->UnsetOwnPrimaryVtx();
1640   delete tmpCascade; tmpCascade=NULL;
1641   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1642     rd2Prong->UnsetOwnPrimaryVtx();
1643   }
1644   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1645   //---
1646
1647
1648   return theCascade;
1649 }
1650
1651
1652 //----------------------------------------------------------------------------
1653 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1654                                    TObjArray *twoTrackArray,AliVEvent *event,
1655                                    AliAODVertex *secVert,
1656                                    AliAODv0 *v0,
1657                                    Double_t dca,
1658                                    Bool_t &okCascades)
1659 {
1660   //
1661   // Make the cascade as a 2Prong decay and check if it passes
1662   // cascades reconstruction cuts
1663   //AliCodeTimerAuto("",0);
1664
1665   //  AliDebug(2,Form("         building the cascade"));
1666   okCascades= kFALSE;
1667   Bool_t dummy1,dummy2,dummy3;
1668
1669   // We use Make2Prong to construct the AliAODRecoCascadeHF
1670   // (which inherits from AliAODRecoDecayHF2Prong)
1671   AliAODRecoCascadeHF *theCascade =
1672     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1673                                      dummy1,dummy2,dummy3);
1674   if(!theCascade) return 0x0;
1675
1676   // bachelor track and charge
1677   AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1678   theCascade->SetCharge(trackBachelor->Charge());
1679
1680   //--- selection cuts
1681   //
1682
1683   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1684   if(fInputAOD){
1685     Int_t idBachelor=(Int_t)trackBachelor->GetID();
1686     if (idBachelor > -1 && idBachelor < fAODMapSize) {
1687       AliAODTrack* trackBachelorAOD=(AliAODTrack*)event->GetTrack(fAODMap[idBachelor]);
1688       tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelorAOD);
1689     }
1690   }else{
1691     tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
1692   }
1693   tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
1694
1695   AliAODVertex *primVertexAOD=0;
1696   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1697     // take event primary vertex
1698     primVertexAOD = PrimaryVertex();
1699     if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
1700     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1701   }
1702
1703   // select Cascades
1704   if(fCascades && fInputAOD){
1705     okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1706   }
1707   else {
1708     //AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied"));
1709     okCascades=kTRUE;
1710   }// no cuts implemented from ESDs
1711   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1712   tmpCascade->UnsetOwnPrimaryVtx();
1713   delete tmpCascade; tmpCascade=NULL;
1714   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1715   //---
1716
1717   return theCascade;
1718 }
1719
1720 //-----------------------------------------------------------------------------
1721 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
1722                                    TObjArray *twoTrackArray,AliVEvent *event,
1723                                    AliAODVertex *secVert,Double_t dca,
1724                                    Bool_t &okD0,Bool_t &okJPSI,
1725                                    Bool_t &okD0fromDstar)
1726 {
1727   // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
1728   // reconstruction cuts
1729   // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
1730   //AliCodeTimerAuto("",0);
1731
1732   okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
1733
1734   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1735   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1736   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
1737
1738   // propagate tracks to secondary vertex, to compute inv. mass
1739   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1740   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1741
1742   Double_t momentum[3];
1743   postrack->GetPxPyPz(momentum);
1744   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1745   negtrack->GetPxPyPz(momentum);
1746   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1747
1748
1749   // invariant mass cut (try to improve coding here..)
1750   Bool_t okMassCut=kFALSE;
1751   if(!okMassCut && fD0toKpi)   if(SelectInvMassAndPtD0Kpi(px,py,pz)) okMassCut=kTRUE;
1752   if(!okMassCut && fJPSItoEle) if(SelectInvMassAndPtJpsiee(px,py,pz)) okMassCut=kTRUE;
1753   if(!okMassCut && fDstar)     if(SelectInvMassAndPtDstarD0pi(px,py,pz)) okMassCut=kTRUE;
1754   if(!okMassCut && fCascades)  if(SelectInvMassAndPtCascade(px,py,pz)) okMassCut=kTRUE;
1755   if(!okMassCut) {
1756     //AliDebug(2," candidate didn't pass mass cut");
1757     return 0x0;
1758   }
1759   // primary vertex to be used by this candidate
1760   AliAODVertex *primVertexAOD  = PrimaryVertex(twoTrackArray,event);
1761   if(!primVertexAOD) return 0x0;
1762
1763
1764   Double_t d0z0[2],covd0z0[3];
1765   postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1766   d0[0] = d0z0[0];
1767   d0err[0] = TMath::Sqrt(covd0z0[0]);
1768   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1769   d0[1] = d0z0[0];
1770   d0err[1] = TMath::Sqrt(covd0z0[0]);
1771
1772   // create the object AliAODRecoDecayHF2Prong
1773   AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
1774   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
1775   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1776   the2Prong->SetProngIDs(2,id);
1777   delete primVertexAOD; primVertexAOD=NULL;
1778
1779
1780   if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
1781
1782     // Add daughter references already here
1783     if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
1784
1785     // select D0->Kpi
1786     if(fD0toKpi)   {
1787       okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event);
1788       if(okD0) the2Prong->SetSelectionBit(AliRDHFCuts::kD0toKpiCuts);
1789     }
1790     //if(fDebug && fD0toKpi) printf("   %d\n",(Int_t)okD0);
1791     // select J/psi from B
1792     if(fJPSItoEle)   {
1793       okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1794     }
1795     //if(fDebug && fJPSItoEle) printf("   %d\n",(Int_t)okJPSI);
1796     // select D0->Kpi from Dstar
1797     if(fDstar)   {
1798       okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
1799       if(okD0fromDstar) the2Prong->SetSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
1800     }
1801     //if(fDebug && fDstar) printf("   %d\n",(Int_t)okD0fromDstar);
1802   }
1803
1804
1805   // remove the primary vertex (was used only for selection)
1806   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1807     the2Prong->UnsetOwnPrimaryVtx();
1808   }
1809
1810   // get PID info from ESD
1811   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1812   if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1813   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1814   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1815   Double_t esdpid[10];
1816   for(Int_t i=0;i<5;i++) {
1817     esdpid[i]   = esdpid0[i];
1818     esdpid[5+i] = esdpid1[i];
1819   }
1820   the2Prong->SetPID(2,esdpid);
1821
1822   return the2Prong;
1823 }
1824 //----------------------------------------------------------------------------
1825 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
1826                              TObjArray *threeTrackArray,AliVEvent *event,
1827                              AliAODVertex *secVert,Double_t dispersion,
1828                              const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
1829                              Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1830                              Bool_t useForLc, Bool_t useForDs, Bool_t &ok3Prong)
1831 {
1832   // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1833   // reconstruction cuts
1834   // E.Bruna, F.Prino
1835
1836   //AliCodeTimerAuto("",0);
1837
1838   ok3Prong=kFALSE;
1839   if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
1840
1841   Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1842   Double_t momentum[3];
1843
1844
1845   AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1846   AliESDtrack *negtrack  = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1847   AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1848
1849   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1850   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1851   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1852   postrack1->GetPxPyPz(momentum);
1853   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1854   negtrack->GetPxPyPz(momentum);
1855   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1856   postrack2->GetPxPyPz(momentum);
1857   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1858
1859   // invariant mass cut for D+, Ds, Lc
1860   Bool_t okMassCut=kFALSE;
1861   if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1862   if(!okMassCut && f3Prong) if(SelectInvMassAndPt3prong(px,py,pz)) okMassCut=kTRUE;
1863   if(!okMassCut) {
1864     //AliDebug(2," candidate didn't pass mass cut");
1865     return 0x0;
1866   }
1867
1868   // primary vertex to be used by this candidate
1869   AliAODVertex *primVertexAOD  = PrimaryVertex(threeTrackArray,event);
1870   if(!primVertexAOD) return 0x0;
1871
1872   Double_t d0z0[2],covd0z0[3];
1873   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1874   d0[0]=d0z0[0];
1875   d0err[0] = TMath::Sqrt(covd0z0[0]);
1876   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1877   d0[1]=d0z0[0];
1878   d0err[1] = TMath::Sqrt(covd0z0[0]);
1879   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1880   d0[2]=d0z0[0];
1881   d0err[2] = TMath::Sqrt(covd0z0[0]);
1882
1883
1884   // create the object AliAODRecoDecayHF3Prong
1885   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1886   Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1887   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]));
1888   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]));
1889   Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
1890   AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1891   the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1892   UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1893   the3Prong->SetProngIDs(3,id);
1894
1895   delete primVertexAOD; primVertexAOD=NULL;
1896
1897   // Add daughter references already here
1898   if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
1899
1900   // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1901   if(f3Prong) {
1902     ok3Prong = kFALSE;
1903
1904     if(fOKInvMassDplus && fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1905       ok3Prong = kTRUE;
1906       the3Prong->SetSelectionBit(AliRDHFCuts::kDplusCuts);
1907     }
1908     if(useForDs && fOKInvMassDs){
1909       if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1910         ok3Prong = kTRUE;
1911         the3Prong->SetSelectionBit(AliRDHFCuts::kDsCuts);
1912       }
1913     }
1914     if(useForLc && fOKInvMassLc){
1915       if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1916         ok3Prong = kTRUE;
1917         the3Prong->SetSelectionBit(AliRDHFCuts::kLcCuts);
1918       }
1919     }
1920   }
1921   //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1922
1923   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1924     the3Prong->UnsetOwnPrimaryVtx();
1925   }
1926
1927   // get PID info from ESD
1928   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1929   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1930   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1931   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1932   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1933   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1934
1935   Double_t esdpid[15];
1936   for(Int_t i=0;i<5;i++) {
1937     esdpid[i]    = esdpid0[i];
1938     esdpid[5+i]  = esdpid1[i];
1939     esdpid[10+i] = esdpid2[i];
1940   }
1941   the3Prong->SetPID(3,esdpid);
1942
1943   return the3Prong;
1944 }
1945 //----------------------------------------------------------------------------
1946 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1947                              TObjArray *fourTrackArray,AliVEvent *event,
1948                              AliAODVertex *secVert,
1949                              const AliAODVertex *vertexp1n1,
1950                              const AliAODVertex *vertexp1n1p2,
1951                              Double_t dcap1n1,Double_t dcap1n2,
1952                              Double_t dcap2n1,Double_t dcap2n2,
1953                              Bool_t &ok4Prong)
1954 {
1955   // Make 4Prong candidates and check if they pass D0toKpipipi
1956   // reconstruction cuts
1957   // G.E.Bruno, R.Romita
1958   //AliCodeTimerAuto("",0);
1959
1960   ok4Prong=kFALSE;
1961   if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
1962
1963   Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1964
1965   AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1966   AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1967   AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1968   AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1969
1970   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1971   negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1972   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1973   negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1974
1975   Double_t momentum[3];
1976   postrack1->GetPxPyPz(momentum);
1977   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1978   negtrack1->GetPxPyPz(momentum);
1979   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1980   postrack2->GetPxPyPz(momentum);
1981   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1982   negtrack2->GetPxPyPz(momentum);
1983   px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1984
1985   // invariant mass cut for rho or D0 (try to improve coding here..)
1986   Bool_t okMassCut=kFALSE;
1987   if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1988   if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) {      //no PID, to be implemented with PID
1989     if(SelectInvMassAndPt4prong(px,py,pz)) okMassCut=kTRUE;
1990   }
1991   if(!okMassCut) {
1992     //if(fDebug) printf(" candidate didn't pass mass cut\n");
1993     //printf(" candidate didn't pass mass cut\n");
1994     return 0x0;
1995   }
1996
1997   // primary vertex to be used by this candidate
1998   AliAODVertex *primVertexAOD  = PrimaryVertex(fourTrackArray,event);
1999   if(!primVertexAOD) return 0x0;
2000
2001   Double_t d0z0[2],covd0z0[3];
2002   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2003   d0[0]=d0z0[0];
2004   d0err[0] = TMath::Sqrt(covd0z0[0]);
2005   negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2006   d0[1]=d0z0[0];
2007   d0err[1] = TMath::Sqrt(covd0z0[0]);
2008   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2009   d0[2]=d0z0[0];
2010   d0err[2] = TMath::Sqrt(covd0z0[0]);
2011   negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2012   d0[3]=d0z0[0];
2013   d0err[3] = TMath::Sqrt(covd0z0[0]);
2014
2015
2016   // create the object AliAODRecoDecayHF4Prong
2017   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2018   Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
2019   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]));
2020   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]));
2021   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]));
2022   Short_t charge=0;
2023   AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
2024   the4Prong->SetOwnPrimaryVtx(primVertexAOD);
2025   UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
2026   the4Prong->SetProngIDs(4,id);
2027
2028   delete primVertexAOD; primVertexAOD=NULL;
2029
2030   ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
2031
2032
2033   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
2034     the4Prong->UnsetOwnPrimaryVtx();
2035   }
2036
2037
2038   // get PID info from ESD
2039   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2040   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2041   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2042   if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
2043   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2044   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2045   Double_t esdpid3[5]={0.,0.,0.,0.,0.};
2046   if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
2047
2048   Double_t esdpid[20];
2049   for(Int_t i=0;i<5;i++) {
2050     esdpid[i]    = esdpid0[i];
2051     esdpid[5+i]  = esdpid1[i];
2052     esdpid[10+i] = esdpid2[i];
2053     esdpid[15+i] = esdpid3[i];
2054   }
2055   the4Prong->SetPID(4,esdpid);
2056
2057   return the4Prong;
2058 }
2059 //-----------------------------------------------------------------------------
2060 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
2061                                                     AliVEvent *event) const
2062 {
2063   // Returns primary vertex to be used for this candidate
2064   //AliCodeTimerAuto("",0);
2065
2066   AliESDVertex *vertexESD = 0;
2067   AliAODVertex *vertexAOD = 0;
2068
2069
2070   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
2071     // primary vertex from the input event
2072
2073     vertexESD = new AliESDVertex(*fV1);
2074
2075   } else {
2076     // primary vertex specific to this candidate
2077
2078     Int_t nTrks = trkArray->GetEntriesFast();
2079     AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
2080
2081     if(fRecoPrimVtxSkippingTrks) {
2082       // recalculating the vertex
2083
2084       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
2085         Float_t diamondcovxy[3];
2086         event->GetDiamondCovXY(diamondcovxy);
2087         Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
2088         Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
2089         AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
2090         vertexer->SetVtxStart(diamond);
2091         delete diamond; diamond=NULL;
2092         if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
2093           vertexer->SetOnlyFitter();
2094       }
2095       Int_t skipped[1000];
2096       Int_t nTrksToSkip=0,id;
2097       AliExternalTrackParam *t = 0;
2098       for(Int_t i=0; i<nTrks; i++) {
2099         t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
2100         id = (Int_t)t->GetID();
2101         if(id<0) continue;
2102         skipped[nTrksToSkip++] = id;
2103       }
2104       // TEMPORARY FIX
2105       // For AOD, skip also tracks without covariance matrix
2106       if(fInputAOD) {
2107         Double_t covtest[21];
2108         for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
2109           AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
2110           if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
2111             id = (Int_t)vtrack->GetID();
2112             if(id<0) continue;
2113             skipped[nTrksToSkip++] = id;
2114           }
2115         }
2116       }
2117       for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
2118       //
2119       vertexer->SetSkipTracks(nTrksToSkip,skipped);
2120       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
2121
2122     } else if(fRmTrksFromPrimVtx && nTrks>0) {
2123       // removing the prongs tracks
2124
2125       TObjArray rmArray(nTrks);
2126       UShort_t *rmId = new UShort_t[nTrks];
2127       AliESDtrack *esdTrack = 0;
2128       AliESDtrack *t = 0;
2129       for(Int_t i=0; i<nTrks; i++) {
2130         t = (AliESDtrack*)trkArray->UncheckedAt(i);
2131         esdTrack = new AliESDtrack(*t);
2132         rmArray.AddLast(esdTrack);
2133         if(esdTrack->GetID()>=0) {
2134           rmId[i]=(UShort_t)esdTrack->GetID();
2135         } else {
2136           rmId[i]=9999;
2137         }
2138       }
2139       Float_t diamondxy[2]={static_cast<Float_t>(event->GetDiamondX()),static_cast<Float_t>(event->GetDiamondY())};
2140       vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
2141       delete [] rmId; rmId=NULL;
2142       rmArray.Delete();
2143
2144     }
2145
2146     if(!vertexESD) return vertexAOD;
2147     if(vertexESD->GetNContributors()<=0) {
2148       //AliDebug(2,"vertexing failed");
2149       delete vertexESD; vertexESD=NULL;
2150       return vertexAOD;
2151     }
2152
2153     delete vertexer; vertexer=NULL;
2154
2155   }
2156
2157   // convert to AliAODVertex
2158   Double_t pos[3],cov[6],chi2perNDF;
2159   vertexESD->GetXYZ(pos); // position
2160   vertexESD->GetCovMatrix(cov); //covariance matrix
2161   chi2perNDF = vertexESD->GetChi2toNDF();
2162   delete vertexESD; vertexESD=NULL;
2163
2164   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
2165
2166   return vertexAOD;
2167 }
2168 //-----------------------------------------------------------------------------
2169 void AliAnalysisVertexingHF::PrintStatus() const {
2170   // Print parameters being used
2171
2172   //printf("Preselections:\n");
2173   //   fTrackFilter->Dump();
2174   if(fSecVtxWithKF) {
2175     printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
2176   } else {
2177     printf("Secondary vertex with AliVertexerTracks\n");
2178   }
2179   if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
2180   if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
2181   if(fD0toKpi) {
2182     printf("Reconstruct D0->Kpi candidates with cuts:\n");
2183     if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
2184   }
2185   if(fDstar) {
2186     printf("Reconstruct D*->D0pi candidates with cuts:\n");
2187     if(fFindVertexForDstar) {
2188       printf("    Reconstruct a secondary vertex for the D*\n");
2189     } else {
2190       printf("    Assume the D* comes from the primary vertex\n");
2191     }
2192     if(fCutsDStartoKpipi) fCutsDStartoKpipi->PrintAll();
2193   }
2194   if(fJPSItoEle) {
2195     printf("Reconstruct J/psi from B candidates with cuts:\n");
2196     if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
2197   }
2198   if(f3Prong) {
2199     printf("Reconstruct 3 prong candidates.\n");
2200     printf("  D+->Kpipi cuts:\n");
2201     if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
2202     printf("  Ds->KKpi cuts:\n");
2203     if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
2204     printf("  Lc->pKpi cuts:\n");
2205     if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
2206   }
2207   if(f4Prong) {
2208     printf("Reconstruct 4 prong candidates.\n");
2209     printf("  D0->Kpipipi cuts:\n");
2210     if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
2211   }
2212   if(fCascades) {
2213     printf("Reconstruct cascades candidates formed with v0s.\n");
2214     printf("  Lc -> k0s P & Lc -> L Pi cuts:\n");
2215     if(fCutsLctoV0) fCutsLctoV0->PrintAll();
2216   }
2217
2218   return;
2219 }
2220 //-----------------------------------------------------------------------------
2221 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
2222                                                                  Double_t &dispersion,Bool_t useTRefArray) const
2223 {
2224   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
2225   //AliCodeTimerAuto("",0);
2226
2227   AliESDVertex *vertexESD = 0;
2228   AliAODVertex *vertexAOD = 0;
2229
2230   if(!fSecVtxWithKF) { // AliVertexerTracks
2231
2232     fVertexerTracks->SetVtxStart(fV1);
2233     vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
2234
2235     if(!vertexESD) return vertexAOD;
2236
2237     if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
2238       //AliDebug(2,"vertexing failed");
2239       delete vertexESD; vertexESD=NULL;
2240       return vertexAOD;
2241     }
2242
2243     Double_t vertRadius2=vertexESD->GetXv()*vertexESD->GetXv()+vertexESD->GetYv()*vertexESD->GetYv();
2244     if(vertRadius2>8.){
2245       // vertex outside beam pipe, reject candidate to avoid propagation through material
2246       delete vertexESD; vertexESD=NULL;
2247       return vertexAOD;
2248     }
2249
2250   } else { // Kalman Filter vertexer (AliKFParticle)
2251
2252     AliKFParticle::SetField(fBzkG);
2253
2254     AliKFVertex vertexKF;
2255
2256     Int_t nTrks = trkArray->GetEntriesFast();
2257     for(Int_t i=0; i<nTrks; i++) {
2258       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
2259       AliKFParticle daughterKF(*esdTrack,211);
2260       vertexKF.AddDaughter(daughterKF);
2261     }
2262     vertexESD = new AliESDVertex(vertexKF.Parameters(),
2263                                  vertexKF.CovarianceMatrix(),
2264                                  vertexKF.GetChi2(),
2265                                  vertexKF.GetNContributors());
2266
2267   }
2268
2269   // convert to AliAODVertex
2270   Double_t pos[3],cov[6],chi2perNDF;
2271   vertexESD->GetXYZ(pos); // position
2272   vertexESD->GetCovMatrix(cov); //covariance matrix
2273   chi2perNDF = vertexESD->GetChi2toNDF();
2274   dispersion = vertexESD->GetDispersion();
2275   delete vertexESD; vertexESD=NULL;
2276
2277   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
2278   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
2279
2280   return vertexAOD;
2281 }
2282 //-----------------------------------------------------------------------------
2283 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt3prong(TObjArray *trkArray){
2284   // Invariant mass cut on tracks
2285   //AliCodeTimerAuto("",0);
2286
2287   Int_t retval=kFALSE;
2288   Double_t momentum[3];
2289   Double_t px[3],py[3],pz[3];
2290   for(Int_t iTrack=0; iTrack<3; iTrack++){
2291     AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
2292     track->GetPxPyPz(momentum);
2293     px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
2294   }
2295   retval = SelectInvMassAndPt3prong(px,py,pz);
2296
2297   return retval;
2298 }
2299
2300 //-----------------------------------------------------------------------------
2301 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt4prong(TObjArray *trkArray){
2302   // Invariant mass cut on tracks
2303   //AliCodeTimerAuto("",0);
2304
2305   Int_t retval=kFALSE;
2306   Double_t momentum[3];
2307   Double_t px[4],py[4],pz[4];
2308
2309   for(Int_t iTrack=0; iTrack<4; iTrack++){
2310     AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
2311     track->GetPxPyPz(momentum);
2312     px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
2313   }
2314
2315   retval = SelectInvMassAndPt4prong(px,py,pz);
2316
2317   return retval;
2318 }
2319 //-----------------------------------------------------------------------------
2320 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtDstarD0pi(TObjArray *trkArray){
2321   // Invariant mass cut on tracks
2322   //AliCodeTimerAuto("",0);
2323
2324   Int_t retval=kFALSE;
2325   Double_t momentum[3];
2326   Double_t px[2],py[2],pz[2];
2327
2328   for(Int_t iTrack=0; iTrack<2; iTrack++){
2329     AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
2330     track->GetPxPyPz(momentum);
2331     px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
2332   }
2333   retval = SelectInvMassAndPtDstarD0pi(px,py,pz);
2334
2335   return retval;
2336 }
2337 //-----------------------------------------------------------------------------
2338 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtD0Kpi(Double_t *px,
2339                                                        Double_t *py,
2340                                                        Double_t *pz){
2341   // Check invariant mass cut and pt candidate cut
2342   //AliCodeTimerAuto("",0);
2343
2344   UInt_t pdg2[2];
2345   Int_t nprongs=2;
2346   Double_t minv2,mrange;
2347   Double_t lolim,hilim;
2348   Double_t minPt=0;
2349   Bool_t retval=kFALSE;
2350
2351   fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2352   fOKInvMassD0=kFALSE;
2353   // pt cut
2354   minPt=fCutsD0toKpi->GetMinPtCandidate();
2355   if(minPt>0.1)
2356     if(fMassCalc2->Pt2() < minPt*minPt) return retval;
2357   // mass cut
2358   mrange=fCutsD0toKpi->GetMassCut();
2359   lolim=fMassDzero-mrange;
2360   hilim=fMassDzero+mrange;
2361   pdg2[0]=211; pdg2[1]=321;
2362   minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2363   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2364     retval=kTRUE;
2365     fOKInvMassD0=kTRUE;
2366   }
2367   pdg2[0]=321; pdg2[1]=211;
2368   minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2369   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2370     retval=kTRUE;
2371     fOKInvMassD0=kTRUE;
2372   }
2373   return retval;
2374 }
2375
2376 //-----------------------------------------------------------------------------
2377 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtJpsiee(Double_t *px,
2378                                                         Double_t *py,
2379                                                         Double_t *pz){
2380   // Check invariant mass cut and pt candidate cut
2381   //AliCodeTimerAuto("",0);
2382
2383   UInt_t pdg2[2];
2384   Int_t nprongs=2;
2385   Double_t minv2,mrange;
2386   Double_t lolim,hilim;
2387   Double_t minPt=0;
2388   Bool_t retval=kFALSE;
2389
2390   fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2391   fOKInvMassJpsi=kFALSE;
2392   // pt cut
2393   minPt=fCutsJpsitoee->GetMinPtCandidate();
2394   if(minPt>0.1)
2395     if(fMassCalc2->Pt2() < minPt*minPt) return retval;
2396   // mass cut
2397   mrange=fCutsJpsitoee->GetMassCut();
2398   lolim=fMassJpsi-mrange;
2399   hilim=fMassJpsi+mrange;
2400
2401   pdg2[0]=11; pdg2[1]=11;
2402   minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2403   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2404     retval=kTRUE;
2405     fOKInvMassJpsi=kTRUE;
2406   }
2407
2408   return retval;
2409 }
2410 //-----------------------------------------------------------------------------
2411 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt3prong(Double_t *px,
2412                                                         Double_t *py,
2413                                                         Double_t *pz,
2414                                                         Int_t pidLcStatus){
2415   // Check invariant mass cut and pt candidate cut
2416   //AliCodeTimerAuto("",0);
2417
2418   UInt_t pdg3[3];
2419   Int_t nprongs=3;
2420   Double_t minv2,mrange;
2421   Double_t lolim,hilim;
2422   Double_t minPt=0;
2423   Bool_t retval=kFALSE;
2424
2425
2426   fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
2427   fOKInvMassDplus=kFALSE;
2428   fOKInvMassDs=kFALSE;
2429   fOKInvMassLc=kFALSE;
2430   // pt cut
2431   minPt=TMath::Min(fCutsDplustoKpipi->GetMinPtCandidate(),fCutsDstoKKpi->GetMinPtCandidate());
2432   minPt=TMath::Min(minPt,fCutsLctopKpi->GetMinPtCandidate());
2433   if(minPt>0.1)
2434     if(fMassCalc3->Pt2() < minPt*minPt) return retval;
2435   // D+->Kpipi
2436   mrange=fCutsDplustoKpipi->GetMassCut();
2437   lolim=fMassDplus-mrange;
2438   hilim=fMassDplus+mrange;
2439   pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
2440   minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2441   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2442     retval=kTRUE;
2443     fOKInvMassDplus=kTRUE;
2444   }
2445   // Ds+->KKpi
2446   mrange=fCutsDstoKKpi->GetMassCut();
2447   lolim=fMassDs-mrange;
2448   hilim=fMassDs+mrange;
2449   pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
2450   minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2451   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2452     retval=kTRUE;
2453     fOKInvMassDs=kTRUE;
2454   }
2455   pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
2456   minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2457   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2458     retval=kTRUE;
2459     fOKInvMassDs=kTRUE;
2460   }
2461   // Lc->pKpi
2462   mrange=fCutsLctopKpi->GetMassCut();
2463   lolim=fMassLambdaC-mrange;
2464   hilim=fMassLambdaC+mrange;
2465   if(pidLcStatus&1){
2466     pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
2467     minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2468     if(minv2>lolim*lolim && minv2<hilim*hilim ){
2469       retval=kTRUE;
2470       fOKInvMassLc=kTRUE;
2471     }
2472   }
2473   if(pidLcStatus&2){
2474     pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
2475     minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
2476     if(minv2>lolim*lolim && minv2<hilim*hilim ){
2477       retval=kTRUE;
2478       fOKInvMassLc=kTRUE;
2479     }
2480   }
2481
2482   return retval;
2483 }
2484
2485 //-----------------------------------------------------------------------------
2486 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtDstarD0pi(Double_t *px,
2487                                                            Double_t *py,
2488                                                            Double_t *pz){
2489   // Check invariant mass cut and pt candidate cut
2490   //AliCodeTimerAuto("",0);
2491
2492   UInt_t pdg2[2];
2493   Int_t nprongs=2;
2494   Double_t minv2,mrange;
2495   Double_t lolim,hilim;
2496   Double_t minPt=0;
2497   Bool_t retval=kFALSE;
2498
2499   fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2500   fOKInvMassDstar=kFALSE;
2501   // pt cut
2502   minPt=fCutsDStartoKpipi->GetMinPtCandidate();
2503   if(minPt>0.1)
2504     if(fMassCalc2->Pt2() < minPt*minPt) return retval;
2505   // mass cut
2506   pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
2507   mrange=fCutsDStartoKpipi->GetMassCut();
2508   lolim=fMassDstar-mrange;
2509   hilim=fMassDstar+mrange;
2510   minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
2511   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2512     retval=kTRUE;
2513     fOKInvMassDstar=kTRUE;
2514   }
2515
2516   return retval;
2517 }
2518
2519 //-----------------------------------------------------------------------------
2520 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt4prong(Double_t *px,
2521                                                         Double_t *py,
2522                                                         Double_t *pz){
2523   // Check invariant mass cut and pt candidate cut
2524   //AliCodeTimerAuto("",0);
2525
2526   UInt_t pdg4[4];
2527   Int_t nprongs=4;
2528   Double_t minv2,mrange;
2529   Double_t lolim,hilim;
2530   Double_t minPt=0;
2531   Bool_t retval=kFALSE;
2532
2533   // D0->Kpipipi without PID
2534   fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
2535   fOKInvMassD0to4p=kFALSE;
2536   // pt cut
2537   minPt=fCutsD0toKpipipi->GetMinPtCandidate();
2538   if(minPt>0.1)
2539     if(fMassCalc4->Pt2() < minPt*minPt) return retval;
2540   // mass cut
2541   mrange=fCutsD0toKpipipi->GetMassCut();
2542   lolim=fMassDzero-mrange;
2543   hilim=fMassDzero+mrange;
2544
2545   pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
2546   minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2547   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2548     retval=kTRUE;
2549     fOKInvMassD0to4p=kTRUE;
2550   }
2551
2552   pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
2553   minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2554   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2555     retval=kTRUE;
2556     fOKInvMassD0to4p=kTRUE;
2557   }
2558
2559   pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
2560   minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2561   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2562     retval=kTRUE;
2563     fOKInvMassD0to4p=kTRUE;
2564   }
2565
2566   pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
2567   minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
2568   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2569     retval=kTRUE;
2570     fOKInvMassD0to4p=kTRUE;
2571   }
2572
2573   return retval;
2574 }
2575 //-----------------------------------------------------------------------------
2576 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPtCascade(Double_t *px,
2577                                                          Double_t *py,
2578                                                          Double_t *pz){
2579   // Check invariant mass cut and pt candidate cut
2580   //AliCodeTimerAuto("",0);
2581
2582   UInt_t pdg2[2];
2583   Int_t nprongs=2;
2584   Double_t minv2,mrange;
2585   Double_t lolim,hilim;
2586   //  Double_t minPt=0;
2587   Bool_t retval=kFALSE;
2588
2589   fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2590   //  minPt=fCutsLctoV0->GetMinPtCandidate();
2591   fOKInvMassLctoV0=kFALSE;
2592   mrange=fCutsLctoV0->GetMassCut();
2593   lolim=fMassLambdaC-mrange;
2594   hilim=fMassLambdaC+mrange;
2595   pdg2[0]=2212;pdg2[1]=310;
2596   minv2=fMassCalc2->InvMass2(2,pdg2);
2597   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2598     retval=kTRUE;
2599     fOKInvMassLctoV0=kTRUE;
2600   }
2601   pdg2[0]=211;pdg2[1]=3122;
2602   minv2=fMassCalc2->InvMass2(2,pdg2);
2603   if(minv2>lolim*lolim && minv2<hilim*hilim ){
2604     retval=kTRUE;
2605     fOKInvMassLctoV0=kTRUE;
2606   }
2607
2608   return retval;
2609 }
2610 //-----------------------------------------------------------------------------
2611 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
2612                                                        Int_t trkEntries,
2613                                                        TObjArray &seleTrksArray,
2614                                                        TObjArray &tracksAtVertex,
2615                                                        Int_t &nSeleTrks,
2616                                                        UChar_t *seleFlags,Int_t *evtNumber)
2617 {
2618   // Apply single-track preselection.
2619   // Fill a TObjArray with selected tracks (for displaced vertices or
2620   // soft pion from D*). Selection flag stored in seleFlags.
2621   // Create the AliESDVertex object (convert from AliAODVertex if necessary)
2622   // In case of AOD input, also fill fAODMap for track index<->ID
2623   //AliCodeTimerAuto("",0);
2624
2625   const AliVVertex *vprimary = event->GetPrimaryVertex();
2626
2627   if(fV1) { delete fV1; fV1=NULL; }
2628   if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
2629
2630   Int_t nindices=0;
2631   UShort_t *indices = 0;
2632   Double_t pos[3],cov[6];
2633   const Int_t entries = event->GetNumberOfTracks();
2634   AliCentrality* cent;
2635
2636   if(!fInputAOD) { // ESD
2637     fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
2638     cent=((AliESDEvent*)event)->GetCentrality();
2639  } else {         // AOD
2640     vprimary->GetXYZ(pos);
2641     vprimary->GetCovarianceMatrix(cov);
2642     fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2643     if(entries<=0) return;
2644     indices = new UShort_t[entries];
2645     memset(indices,0,sizeof(UShort_t)*entries);
2646     fAODMapSize = 100000;
2647     fAODMap = new Int_t[fAODMapSize];
2648     memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
2649     cent=((AliAODEvent*)event)->GetCentrality();
2650   }
2651   Float_t centperc=cent->GetCentralityPercentile("V0M");
2652
2653   Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE,okFor3Prong=kFALSE;
2654   nSeleTrks=0;
2655
2656   // transfer ITS tracks from event to arrays
2657   for(Int_t i=0; i<entries; i++) {
2658     AliVTrack *track;
2659     track = (AliVTrack*)event->GetTrack(i);
2660
2661     // skip pure ITS SA tracks
2662     if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2663
2664     // skip tracks without ITS
2665     if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
2666
2667     // skip tracks with negative ID
2668     // (these are duplicated TPC-only AOD tracks, for jet analysis...)
2669     if(track->GetID()<0) continue;
2670
2671     // TEMPORARY: check that the cov matrix is there
2672     Double_t covtest[21];
2673     if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2674     //
2675
2676     if(fInputAOD) {
2677       AliAODTrack *aodt = (AliAODTrack*)track;
2678       if(aodt->GetUsedForPrimVtxFit()) {
2679         indices[nindices]=aodt->GetID(); nindices++;
2680       }
2681       Int_t ind = (Int_t)aodt->GetID();
2682       if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
2683     }
2684
2685     AliESDtrack *esdt = 0;
2686
2687     if(!fInputAOD) {
2688       esdt = (AliESDtrack*)track;
2689     } else {
2690       esdt = new AliESDtrack(track);
2691     }
2692
2693     // single track selection
2694     okDisplaced=kFALSE; okSoftPi=kFALSE; okFor3Prong=kFALSE;
2695     if(fMixEvent && i<trkEntries){
2696       evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
2697       const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
2698       Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
2699       eventVtx->GetXYZ(vtxPos);
2700       vprimary->GetXYZ(primPos);
2701       eventVtx->GetCovarianceMatrix(primCov);
2702       for(Int_t ind=0;ind<3;ind++){
2703         trasl[ind]=vtxPos[ind]-primPos[ind];
2704       }
2705
2706       Bool_t isTransl=esdt->Translate(trasl,primCov);
2707       if(!isTransl) {
2708         delete esdt;
2709         esdt = NULL;
2710         continue;
2711       }
2712     }
2713
2714     if(SingleTrkCuts(esdt,centperc,okDisplaced,okSoftPi,okFor3Prong) && nSeleTrks<trkEntries) {
2715       esdt->PropagateToDCA(fV1,fBzkG,kVeryBig);
2716       seleTrksArray.AddLast(esdt);
2717       tracksAtVertex.AddLast(new AliExternalTrackParam(*esdt));
2718       seleFlags[nSeleTrks]=0;
2719       if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2720       if(okFor3Prong)    SETBIT(seleFlags[nSeleTrks],kBit3Prong);
2721       if(okSoftPi)    SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2722
2723       // Check the PID
2724       SETBIT(seleFlags[nSeleTrks],kBitPionCompat);
2725       SETBIT(seleFlags[nSeleTrks],kBitKaonCompat);
2726       SETBIT(seleFlags[nSeleTrks],kBitProtonCompat);
2727       Bool_t useTPC=kTRUE;
2728       if(fUseTOFPID){
2729         Double_t nsigmatofPi= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kPion);
2730         if(nsigmatofPi>-990. && (nsigmatofPi<-fnSigmaTOFPionLow || nsigmatofPi>fnSigmaTOFPionHi)){
2731           CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
2732         }
2733         Double_t nsigmatofK= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kKaon);
2734         if(nsigmatofK>-990. && (nsigmatofK<-fnSigmaTOFKaonLow || nsigmatofK>fnSigmaTOFKaonHi)){
2735           CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
2736         }
2737         Double_t nsigmatofP= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kProton);
2738         if(nsigmatofP>-990. && (nsigmatofP<-fnSigmaTOFProtonLow || nsigmatofP>fnSigmaTOFProtonHi)){
2739           CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
2740         }
2741         if(fUseTPCPIDOnlyIfNoTOF && nsigmatofPi>-990.) useTPC=kFALSE;
2742       }
2743       if(useTPC && fUseTPCPID && esdt->P()<fMaxMomForTPCPid){
2744         Double_t nsigmatpcPi= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kPion);
2745         if(nsigmatpcPi>-990. && (nsigmatpcPi<-fnSigmaTPCPionLow || nsigmatpcPi>fnSigmaTPCPionHi)){
2746           CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
2747         }
2748         Double_t nsigmatpcK= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kKaon);
2749         if(nsigmatpcK>-990. && (nsigmatpcK<-fnSigmaTPCKaonLow || nsigmatpcK>fnSigmaTPCKaonHi)){
2750           CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
2751         }
2752         Double_t nsigmatpcP= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kProton);
2753         if(nsigmatpcP>-990. && (nsigmatpcP<-fnSigmaTPCProtonLow || nsigmatpcP>fnSigmaTPCProtonHi)){
2754           CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
2755         }
2756       }
2757       nSeleTrks++;
2758     } else {
2759       if(fInputAOD) delete esdt;
2760       esdt = NULL;
2761       continue;
2762     }
2763
2764   } // end loop on tracks
2765
2766   // primary vertex from AOD
2767   if(fInputAOD) {
2768     delete fV1; fV1=NULL;
2769     vprimary->GetXYZ(pos);
2770     vprimary->GetCovarianceMatrix(cov);
2771     Double_t chi2toNDF = vprimary->GetChi2perNDF();
2772     Int_t ncontr=nindices;
2773     if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
2774     Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
2775     fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2776     fV1->SetTitle(vprimary->GetTitle());
2777     fV1->SetIndices(nindices,indices);
2778   }
2779   if(indices) { delete [] indices; indices=NULL; }
2780
2781
2782   return;
2783 }
2784 //-----------------------------------------------------------------------------
2785 void AliAnalysisVertexingHF::SetSelectionBitForPID(AliRDHFCuts *cuts,AliAODRecoDecayHF *rd,Int_t bit) {
2786   //
2787   // Set the selection bit for PID
2788   //
2789   //AliCodeTimerAuto("",0);
2790   if(cuts->GetPidHF()) {
2791     Bool_t usepid=cuts->GetIsUsePID();
2792     cuts->SetUsePID(kTRUE);
2793     if(cuts->IsSelectedPID(rd))
2794       rd->SetSelectionBit(bit);
2795     cuts->SetUsePID(usepid);
2796   }
2797   return;
2798 }
2799 //-----------------------------------------------------------------------------
2800 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2801                                              Float_t centralityperc,
2802                                              Bool_t &okDisplaced,
2803                                              Bool_t &okSoftPi,
2804                                              Bool_t &okFor3Prong) const
2805 {
2806   // Check if track passes some kinematical cuts
2807
2808   // this is needed to store the impact parameters
2809   //AliCodeTimerAuto("",0);
2810
2811   trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2812
2813   UInt_t selectInfo;
2814   //
2815   // Track selection, displaced tracks -- 2 prong
2816   selectInfo = 0;
2817   if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
2818      && fTrackFilter2prongCentral && centralityperc<fMaxCentPercentileForTightCuts){
2819     // central PbPb events, tighter cuts
2820     selectInfo = fTrackFilter2prongCentral->IsSelected(trk);
2821   }else{
2822     // standard cuts
2823     if(fTrackFilter) {
2824       selectInfo = fTrackFilter->IsSelected(trk);
2825     }
2826   }
2827   if(selectInfo) okDisplaced=kTRUE;
2828
2829   // Track selection, displaced tracks -- 3 prong
2830   selectInfo = 0;
2831   if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
2832      && fTrackFilter3prongCentral && centralityperc<fMaxCentPercentileForTightCuts){
2833     // central PbPb events, tighter cuts
2834     selectInfo = fTrackFilter3prongCentral->IsSelected(trk);
2835   }else{
2836     // standard cuts
2837     if(fTrackFilter) {
2838       selectInfo = fTrackFilter->IsSelected(trk);
2839     }
2840   }
2841   if(selectInfo) okFor3Prong=kTRUE;
2842
2843   // Track selection, soft pions
2844   selectInfo = 0;
2845   if(fDstar && fTrackFilterSoftPi) {
2846     selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2847   }
2848   if(selectInfo) okSoftPi=kTRUE;
2849
2850   if(okDisplaced || okSoftPi || okFor3Prong) return kTRUE;
2851
2852   return kFALSE;
2853 }
2854
2855
2856 //-----------------------------------------------------------------------------
2857 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2858   //
2859   // Transform ESDv0 to AODv0
2860   //
2861   //  this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2862   //  and creates an AODv0 out of them
2863   //
2864   //AliCodeTimerAuto("",0);
2865   Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
2866   AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2867
2868   // create the v0 neutral track to compute the DCA to the primary vertex
2869   Double_t xyz[3], pxpypz[3];
2870   esdV0->XvYvZv(xyz);
2871   esdV0->PxPyPz(pxpypz);
2872   Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2873   AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2874   if(!trackesdV0) {
2875     delete vertexV0;
2876     return 0;
2877   }
2878   Double_t d0z0[2],covd0z0[3];
2879   AliAODVertex *primVertexAOD = PrimaryVertex();
2880   trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2881   Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2882   // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2883   Double_t dcaV0DaughterToPrimVertex[2];
2884   AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2885   AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
2886   if( !posV0track || !negV0track) {
2887     if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2888     delete vertexV0;
2889     delete primVertexAOD;
2890     return 0;
2891   }
2892   posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2893   //  if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2894   //  else
2895   dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
2896   negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2897   //  if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2898   //  else
2899   dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
2900   Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
2901   Double_t pmom[3],nmom[3];
2902   esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2903   esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2904
2905   AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2906   aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2907
2908   delete trackesdV0;
2909   delete primVertexAOD;
2910
2911   return aodV0;
2912 }
2913 //-----------------------------------------------------------------------------
2914 void AliAnalysisVertexingHF::SetParametersAtVertex(AliESDtrack* esdt, const AliExternalTrackParam* extpar) const{
2915   // Set the stored track parameters at primary vertex into AliESDtrack
2916   //AliCodeTimerAuto("",0);
2917
2918   const Double_t *par=extpar->GetParameter();
2919   const Double_t *cov=extpar->GetCovariance();
2920   Double_t alpha=extpar->GetAlpha();
2921   Double_t x=extpar->GetX();
2922   esdt->Set(x,alpha,par,cov);
2923   return;
2924 }
2925 //-----------------------------------------------------------------------------
2926 void AliAnalysisVertexingHF::SetMasses(){
2927   // Set the hadron mass values from TDatabasePDG
2928
2929   fMassDzero=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2930   fMassDplus=TDatabasePDG::Instance()->GetParticle(411)->Mass();
2931   fMassDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
2932   fMassLambdaC=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
2933   fMassDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
2934   fMassJpsi=TDatabasePDG::Instance()->GetParticle(443)->Mass();
2935 }
2936 //-----------------------------------------------------------------------------
2937 Bool_t AliAnalysisVertexingHF::CheckCutsConsistency(){
2938   //
2939   // Check the Vertexer and the analysts task consitstecny
2940   //
2941
2942
2943   //___ Check if the V0 type from AliRDHFCutsLctoV0 is the same as the one set in the ConfigVertexingHF.C for AliAnalysisVertexingHF
2944
2945
2946   if ( fCutsLctoV0 && fV0TypeForCascadeVertex != fCutsLctoV0->GetV0Type())
2947   {
2948     printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHFCutsLctoV0 (%d)\n",fV0TypeForCascadeVertex,fCutsLctoV0->GetV0Type());
2949     return kFALSE;
2950   }
2951   return kTRUE;
2952 }
2953 //-----------------------------------------------------------------------------
2954