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