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