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