PWG3/vertexingHF/AliAnalysisVertexingHF.cxx

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