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