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