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