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           // Set selection bit for PID
 694           if(okD0) SetSelectionBitForPID(fCutsD0toKpi,rd);
 695         }
 696         // D* candidates
 697         if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
 698           // write references in io2Prong
 699           if(fInputAOD) {
 700             AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
 701           } else {
 702             vertexp1n1->AddDaughter(postrack1);
 703             vertexp1n1->AddDaughter(negtrack1);
 704           }
 705           io2Prong->SetSecondaryVtx(vertexp1n1);
 706           //printf("--->  %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
 707           // create a track from the D0
 708           AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
 709
 710           // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
 711           for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
 712
 713             if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
 714
 715             if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
 716
 717             if(fMixEvent) {
 718               if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] ||
 719                  evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] ||
 720                  evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
 721             }
 722
 723             if(iTrkSoftPi%1==0) AliDebug(1,Form("    1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));
 724
 725             trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
 726             if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
 727
 728             // get track from tracks array
 729             trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
 730             trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
 731
 732             twoTrackArrayCasc->AddAt(trackPi,0);
 733             twoTrackArrayCasc->AddAt(trackD0,1);
 734
 735             AliAODVertex *vertexCasc = 0;
 736
 737             if(fFindVertexForDstar) {
 738               // DCA between the two tracks
 739               dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
 740               // Vertexing
 741               vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
 742             } else {
 743               // assume Dstar decays at the primary vertex
 744               Double_t pos[3],cov[6],chi2perNDF;
 745               fV1->GetXYZ(pos);
 746               fV1->GetCovMatrix(cov);
 747               chi2perNDF = fV1->GetChi2toNDF();
 748               vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
 749               dcaCasc = 0.;
 750             }
 751             if(!vertexCasc) {
 752               twoTrackArrayCasc->Clear();
 753               trackPi=0;
 754               continue;
 755             }
 756
 757             ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
 758             if(okDstar) {
 759               // add the D0 to the AOD (if not already done)
 760               if(!okD0) {
 761                 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
 762                 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
 763                 rd->SetSecondaryVtx(v2Prong);
 764                 v2Prong->SetParent(rd);
 765                 AddRefs(v2Prong,rd,event,twoTrackArray1);
 766                 okD0=kTRUE; // this is done to add it only once
 767                 // Set selection bit for PID
 768                 SetSelectionBitForPID(fCutsD0toKpi,rd);
 769               }
 770               // add the vertex and the cascade to the AOD
 771               AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
 772               rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
 773               rc->SetSecondaryVtx(vCasc);
 774               vCasc->SetParent(rc);
 775               if(!fInputAOD) vCasc->AddDaughter(rd); // just to fill ref #0
 776               AddRefs(vCasc,rc,event,twoTrackArrayCasc);
 777               vCasc->AddDaughter(rd); // add the D0 (in ref #1)
 778               // Set selection bit for PID
 779               SetSelectionBitForPID(fCutsDStartoKpipi,rc);
 780             }
 781             twoTrackArrayCasc->Clear();
 782             trackPi=0;
 783             if(ioCascade) {delete ioCascade; ioCascade=NULL;}
 784             delete vertexCasc; vertexCasc=NULL;
 785           } // end loop on soft pi tracks
 786
 787           if(trackD0) {delete trackD0; trackD0=NULL;}
 788
 789         }
 790         if(io2Prong) {delete io2Prong; io2Prong=NULL;}
 791       }
 792
 793       twoTrackArray1->Clear();
 794       if( (!f3Prong && !f4Prong) ||
 795           (isLikeSign2Prong && !f3Prong) ) {
 796         negtrack1=0;
 797         delete vertexp1n1;
 798         continue;
 799       }
 800
 801
 802       // 2nd LOOP  ON  POSITIVE  TRACKS
 803       for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
 804
 805         if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
 806
 807         if(iTrkP2%1==0) AliDebug(1,Form("    2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));
 808
 809         // get track from tracks array
 810         postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
 811
 812         if(postrack2->Charge()<0) continue;
 813
 814         if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
 815
 816         if(fMixEvent) {
 817           if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
 818              evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
 819              evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
 820         }
 821
 822         if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
 823           if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
 824             isLikeSign3Prong=kTRUE;
 825           } else { // not ok
 826             continue;
 827           }
 828         } else { // normal triplet (+-+)
 829           isLikeSign3Prong=kFALSE;
 830           if(fMixEvent) {
 831             if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
 832                evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
 833                evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
 834           }
 835         }
 836
 837         // back to primary vertex
 838         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 839         postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 840         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 841         //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
 842
 843         dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
 844         if(dcap2n1>dcaMax) { postrack2=0; continue; }
 845         dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
 846         if(dcap1p2>dcaMax) { postrack2=0; continue; }
 847
 848         // check invariant mass cuts for D+,Ds,Lc
 849         massCutOK=kTRUE;
 850         if(f3Prong) {
 851           if(postrack2->Charge()>0) {
 852             threeTrackArray->AddAt(postrack1,0);
 853             threeTrackArray->AddAt(negtrack1,1);
 854             threeTrackArray->AddAt(postrack2,2);
 855           } else {
 856             threeTrackArray->AddAt(negtrack1,0);
 857             threeTrackArray->AddAt(postrack1,1);
 858             threeTrackArray->AddAt(postrack2,2);
 859           }
 860           if(fMassCutBeforeVertexing)
 861             massCutOK = SelectInvMass(threeTrackArray);
 862         }
 863
 864         if(f3Prong && !massCutOK) {
 865           threeTrackArray->Clear();
 866           if(!f4Prong) {
 867             postrack2=0;
 868             continue;
 869           }
 870         }
 871
 872         // Vertexing
 873         twoTrackArray2->AddAt(postrack2,0);
 874         twoTrackArray2->AddAt(negtrack1,1);
 875         AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
 876         if(!vertexp2n1) {
 877           twoTrackArray2->Clear();
 878           postrack2=0;
 879           continue;
 880         }
 881
 882         // 3 prong candidates
 883         if(f3Prong && massCutOK) {
 884
 885           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
 886           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
 887           if(ok3Prong) {
 888             AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
 889             if(!isLikeSign3Prong) {
 890               rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
 891               rd->SetSecondaryVtx(v3Prong);
 892               v3Prong->SetParent(rd);
 893               AddRefs(v3Prong,rd,event,threeTrackArray);
 894             } else { // isLikeSign3Prong
 895               rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
 896               rd->SetSecondaryVtx(v3Prong);
 897               v3Prong->SetParent(rd);
 898               AddRefs(v3Prong,rd,event,threeTrackArray);
 899             }
 900             // Set selection bit for PID
 901             SetSelectionBitForPID(fCutsDplustoKpipi,rd);
 902             SetSelectionBitForPID(fCutsDstoKKpi,rd);
 903             SetSelectionBitForPID(fCutsLctopKpi,rd);
 904           }
 905           if(io3Prong) {delete io3Prong; io3Prong=NULL;}
 906           if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
 907         }
 908
 909         // 4 prong candidates
 910         if(f4Prong
 911            // don't make 4 prong with like-sign pairs and triplets
 912            && !isLikeSign2Prong && !isLikeSign3Prong
 913            // track-to-track dca cuts already now
 914            && dcap1n1 < fCutsD0toKpipipi->GetDCACut()
 915            && dcap2n1 < fCutsD0toKpipipi->GetDCACut()) {
 916
 917           // back to primary vertex
 918           postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 919           postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 920           negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 921           // Vertexing for these 3 (can be taken from above?)
 922           threeTrackArray->AddAt(postrack1,0);
 923           threeTrackArray->AddAt(negtrack1,1);
 924           threeTrackArray->AddAt(postrack2,2);
 925           AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
 926
 927           // 3rd LOOP  ON  NEGATIVE  TRACKS (for 4 prong)
 928           for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
 929
 930             if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
 931
 932             if(iTrkN2%1==0) AliDebug(1,Form("    3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
 933
 934             // get track from tracks array
 935             negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
 936
 937             if(negtrack2->Charge()>0) continue;
 938
 939             if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
 940             if(fMixEvent){
 941               if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
 942                  evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
 943                  evtNumber[iTrkP2]==evtNumber[iTrkN2] ||
 944                  evtNumber[iTrkP1]==evtNumber[iTrkN1] ||
 945                  evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
 946                  evtNumber[iTrkN1]==evtNumber[iTrkP2]) continue;
 947             }
 948
 949             // back to primary vertex
 950             postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 951             postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 952             negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
 953             negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
 954             dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 955             if(dcap1n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
 956             dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
 957             if(dcap2n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
 958
 959
 960             fourTrackArray->AddAt(postrack1,0);
 961             fourTrackArray->AddAt(negtrack1,1);
 962             fourTrackArray->AddAt(postrack2,2);
 963             fourTrackArray->AddAt(negtrack2,3);
 964
 965             // check invariant mass cuts for D0
 966             massCutOK=kTRUE;
 967             if(fMassCutBeforeVertexing)
 968               massCutOK = SelectInvMass(fourTrackArray);
 969
 970             if(!massCutOK) {
 971               fourTrackArray->Clear();
 972               negtrack2=0;
 973               continue;
 974             }
 975
 976             // Vertexing
 977             AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
 978             io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
 979             if(ok4Prong) {
 980               AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
 981               rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
 982               rd->SetSecondaryVtx(v4Prong);
 983               v4Prong->SetParent(rd);
 984               AddRefs(v4Prong,rd,event,fourTrackArray);
 985             }
 986
 987             if(io4Prong) {delete io4Prong; io4Prong=NULL;}
 988             if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;}
 989             fourTrackArray->Clear();
 990             negtrack2 = 0;
 991
 992           } // end loop on negative tracks
 993
 994           threeTrackArray->Clear();
 995           delete vertexp1n1p2;
 996
 997         }
 998
 999         postrack2 = 0;
1000         delete vertexp2n1;
1001
1002       } // end 2nd loop on positive tracks
1003
1004       twoTrackArray2->Clear();
1005
1006       // 2nd LOOP  ON  NEGATIVE  TRACKS (for 3 prong -+-)
1007       for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
1008
1009         if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1010
1011         if(iTrkN2%1==0) AliDebug(1,Form("    2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
1012
1013         // get track from tracks array
1014         negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1015
1016         if(negtrack2->Charge()>0) continue;
1017
1018         if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1019
1020         if(fMixEvent) {
1021           if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
1022              evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
1023              evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1024         }
1025
1026         if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
1027           if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
1028             isLikeSign3Prong=kTRUE;
1029           } else { // not ok
1030             continue;
1031           }
1032         } else { // normal triplet (-+-)
1033           isLikeSign3Prong=kFALSE;
1034         }
1035
1036         // back to primary vertex
1037         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1038         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1039         negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1040         //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
1041
1042         dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1043         if(dcap1n2>dcaMax) { negtrack2=0; continue; }
1044         dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1045         if(dcan1n2>dcaMax) { negtrack2=0; continue; }
1046
1047         threeTrackArray->AddAt(negtrack1,0);
1048         threeTrackArray->AddAt(postrack1,1);
1049         threeTrackArray->AddAt(negtrack2,2);
1050
1051         // check invariant mass cuts for D+,Ds,Lc
1052         massCutOK=kTRUE;
1053         if(fMassCutBeforeVertexing && f3Prong)
1054           massCutOK = SelectInvMass(threeTrackArray);
1055
1056         if(!massCutOK) {
1057           threeTrackArray->Clear();
1058           negtrack2=0;
1059           continue;
1060         }
1061
1062         // Vertexing
1063         twoTrackArray2->AddAt(postrack1,0);
1064         twoTrackArray2->AddAt(negtrack2,1);
1065
1066         AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1067         if(!vertexp1n2) {
1068           twoTrackArray2->Clear();
1069           negtrack2=0;
1070           continue;
1071         }
1072
1073         if(f3Prong) {
1074           AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1075           io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
1076           if(ok3Prong) {
1077             AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1078             if(!isLikeSign3Prong) {
1079               rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1080               rd->SetSecondaryVtx(v3Prong);
1081               v3Prong->SetParent(rd);
1082               AddRefs(v3Prong,rd,event,threeTrackArray);
1083             } else { // isLikeSign3Prong
1084               rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1085               rd->SetSecondaryVtx(v3Prong);
1086               v3Prong->SetParent(rd);
1087               AddRefs(v3Prong,rd,event,threeTrackArray);
1088             }
1089             // Set selection bit for PID
1090             SetSelectionBitForPID(fCutsDplustoKpipi,rd);
1091             SetSelectionBitForPID(fCutsDstoKKpi,rd);
1092             SetSelectionBitForPID(fCutsLctopKpi,rd);
1093           }
1094           if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1095           if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1096         }
1097         threeTrackArray->Clear();
1098         negtrack2 = 0;
1099         delete vertexp1n2;
1100
1101       } // end 2nd loop on negative tracks
1102
1103       twoTrackArray2->Clear();
1104
1105       negtrack1 = 0;
1106       delete vertexp1n1;
1107     } // end 1st loop on negative tracks
1108
1109     postrack1 = 0;
1110   }  // end 1st loop on positive tracks
1111
1112
1113   AliDebug(1,Form(" Total HF vertices in event = %d;",
1114                   (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
1115   if(fD0toKpi) {
1116     AliDebug(1,Form(" D0->Kpi in event = %d;",
1117                     (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
1118   }
1119   if(fJPSItoEle) {
1120     AliDebug(1,Form(" JPSI->ee in event = %d;",
1121                     (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
1122   }
1123   if(f3Prong) {
1124     AliDebug(1,Form(" Charm->3Prong in event = %d;",
1125                     (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
1126   }
1127   if(f4Prong) {
1128     AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
1129                     (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
1130   }
1131   if(fDstar) {
1132     AliDebug(1,Form(" D*->D0pi in event = %d;\n",
1133                     (Int_t)aodDstarTClArr->GetEntriesFast()));
1134   }
1135   if(fCascades){
1136     AliDebug(1,Form(" cascades -> v0 + track in event = %d;\n",
1137                     (Int_t)aodCascadesTClArr->GetEntriesFast()));
1138   }
1139   if(fLikeSign) {
1140     AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
1141                     (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
1142   }
1143   if(fLikeSign && f3Prong) {
1144     AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
1145                     (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
1146   }
1147
1148
1149   twoTrackArray1->Delete();  delete twoTrackArray1;
1150   twoTrackArray2->Delete();  delete twoTrackArray2;
1151   twoTrackArrayCasc->Delete();  delete twoTrackArrayCasc;
1152   twoTrackArrayV0->Delete();  delete twoTrackArrayV0;
1153   threeTrackArray->Clear();
1154   threeTrackArray->Delete(); delete threeTrackArray;
1155   fourTrackArray->Delete();  delete fourTrackArray;
1156   delete [] seleFlags; seleFlags=NULL;
1157   if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
1158
1159   if(fInputAOD) {
1160     seleTrksArray.Delete();
1161     if(fAODMap) { delete fAODMap; fAODMap=NULL; }
1162   }
1163
1164   return;
1165 }
1166 //----------------------------------------------------------------------------
1167 void AliAnalysisVertexingHF::AddRefs(AliAODVertex *v,AliAODRecoDecayHF *rd,
1168                                      const AliVEvent *event,
1169                                      const TObjArray *trkArray) const
1170 {
1171   // Add the AOD tracks as daughters of the vertex (TRef)
1172   // Also add the references to the primary vertex and to the cuts
1173
1174   if(fInputAOD) {
1175     AddDaughterRefs(v,event,trkArray);
1176     rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1177   }
1178
1179   /*
1180   rd->SetListOfCutsRef((TList*)fListOfCuts);
1181   //fListOfCuts->Print();
1182   cout<<fListOfCuts<<endl;
1183   TList *l=(TList*)rd->GetListOfCuts();
1184   cout<<l<<endl;
1185   if(l) {l->Print(); }else{printf("error\n");}
1186   */
1187
1188   return;
1189 }
1190 //----------------------------------------------------------------------------
1191 void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,
1192                                              const AliVEvent *event,
1193                                              const TObjArray *trkArray) const
1194 {
1195   // Add the AOD tracks as daughters of the vertex (TRef)
1196
1197   Int_t nDg = v->GetNDaughters();
1198   TObject *dg = 0;
1199   if(nDg) dg = v->GetDaughter(0);
1200
1201   if(dg) return; // daughters already added
1202
1203   Int_t nTrks = trkArray->GetEntriesFast();
1204
1205   AliExternalTrackParam *track = 0;
1206   AliAODTrack *aodTrack = 0;
1207   Int_t id;
1208
1209   for(Int_t i=0; i<nTrks; i++) {
1210     track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1211     id = (Int_t)track->GetID();
1212     //printf("---> %d\n",id);
1213     if(id<0) continue; // this track is a AliAODRecoDecay
1214     aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
1215     v->AddDaughter(aodTrack);
1216   }
1217
1218   return;
1219 }
1220 //---------------------------------------------------------------------------
1221 void AliAnalysisVertexingHF::FixReferences(AliAODEvent *aod)
1222 {
1223   // Checks that the references to the daughter tracks are properly
1224   // assigned and reassigns them if needed
1225   //
1226
1227
1228   TClonesArray *inputArray=(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
1229   if(!inputArray) return;
1230
1231   AliAODTrack *track = 0;
1232   AliAODVertex *vertex = 0;
1233
1234   Bool_t needtofix=kFALSE;
1235   for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1236     vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1237     for(Int_t id=0; id<vertex->GetNDaughters(); id++) {
1238       track = (AliAODTrack*)vertex->GetDaughter(id);
1239       if(!track->GetStatus()) needtofix=kTRUE;
1240     }
1241     if(needtofix) break;
1242   }
1243
1244   if(!needtofix) return;
1245
1246
1247   printf("Fixing references\n");
1248
1249   fAODMapSize = 100000;
1250   fAODMap = new Int_t[fAODMapSize];
1251
1252   for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
1253     track = aod->GetTrack(i);
1254
1255     // skip pure ITS SA tracks
1256     if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
1257
1258     // skip tracks without ITS
1259     if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
1260
1261     // TEMPORARY: check that the cov matrix is there
1262     Double_t covtest[21];
1263     if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1264     //
1265
1266     fAODMap[(Int_t)track->GetID()] = i;
1267   }
1268
1269
1270   Int_t ids[4]={-1,-1,-1,-1};
1271   for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1272     Bool_t cascade=kFALSE;
1273     vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1274     Int_t id=0;
1275     Int_t nDgs = vertex->GetNDaughters();
1276     for(id=0; id<nDgs; id++) {
1277       track = (AliAODTrack*)vertex->GetDaughter(id);
1278       if(track->Charge()==0) {cascade=kTRUE; continue;} // cascade
1279       ids[id]=(Int_t)track->GetID();
1280       vertex->RemoveDaughter(track);
1281     }
1282     if(cascade) continue;
1283     for(id=0; id<nDgs; id++) {
1284       track = aod->GetTrack(fAODMap[ids[id]]);
1285       vertex->AddDaughter(track);
1286     }
1287
1288   }
1289
1290   return;
1291 }
1292 //----------------------------------------------------------------------------
1293 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1294                                    TObjArray *twoTrackArray,AliVEvent *event,
1295                                    AliAODVertex *secVert,
1296                                    AliAODRecoDecayHF2Prong *rd2Prong,
1297                                    Double_t dca,
1298                                    Bool_t &okDstar) const
1299 {
1300   // Make the cascade as a 2Prong decay and check if it passes Dstar
1301   // reconstruction cuts
1302
1303   okDstar = kFALSE;
1304
1305   Bool_t dummy1,dummy2,dummy3;
1306
1307   // We use Make2Prong to construct the AliAODRecoCascadeHF
1308   // (which inherits from AliAODRecoDecayHF2Prong)
1309   AliAODRecoCascadeHF *theCascade =
1310     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1311                                      dummy1,dummy2,dummy3);
1312   if(!theCascade) return 0x0;
1313
1314   // charge
1315   AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1316   theCascade->SetCharge(trackPi->Charge());
1317
1318   //--- selection cuts
1319   //
1320   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1321   tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
1322   tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
1323   AliAODVertex *primVertexAOD=0;
1324   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1325     // take event primary vertex
1326     primVertexAOD = PrimaryVertex();
1327     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1328     rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
1329   }
1330   // select D*->D0pi
1331   if(fDstar) {
1332     okDstar = (Bool_t)fCutsDStartoKpipi->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1333     if(okDstar) theCascade->SetSelectionBit(AliRDHFCuts::kDstarCuts);
1334   }
1335   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1336   tmpCascade->UnsetOwnPrimaryVtx();
1337   delete tmpCascade; tmpCascade=NULL;
1338   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1339     rd2Prong->UnsetOwnPrimaryVtx();
1340   }
1341   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1342   //---
1343
1344
1345   return theCascade;
1346 }
1347
1348
1349 //----------------------------------------------------------------------------
1350 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1351                                    TObjArray *twoTrackArray,AliVEvent *event,
1352                                    AliAODVertex *secVert,
1353                                    AliAODv0 *v0,
1354                                    Double_t dca,
1355                                    Bool_t &okCascades) const
1356 {
1357   //
1358   // Make the cascade as a 2Prong decay and check if it passes
1359   // cascades reconstruction cuts
1360
1361   //  AliDebug(2,Form("         building the cascade"));
1362   okCascades= kFALSE;
1363   Bool_t dummy1,dummy2,dummy3;
1364
1365   // We use Make2Prong to construct the AliAODRecoCascadeHF
1366   // (which inherits from AliAODRecoDecayHF2Prong)
1367   AliAODRecoCascadeHF *theCascade =
1368     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1369                                      dummy1,dummy2,dummy3);
1370   if(!theCascade) return 0x0;
1371
1372   // bachelor track and charge
1373   AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1374   theCascade->SetCharge(trackBachelor->Charge());
1375
1376   //--- selection cuts
1377   //
1378   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1379   tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
1380   tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
1381   AliAODVertex *primVertexAOD=0;
1382   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1383     // take event primary vertex
1384     primVertexAOD = PrimaryVertex();
1385     if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
1386     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1387   }
1388
1389   // select Cascades
1390   if(fCascades && fInputAOD){
1391     okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1392   }
1393   else { AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied")); okCascades=kTRUE; }// no cuts implemented from ESDs
1394   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1395   tmpCascade->UnsetOwnPrimaryVtx();
1396   delete tmpCascade; tmpCascade=NULL;
1397   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1398   //---
1399
1400   return theCascade;
1401 }
1402
1403 //-----------------------------------------------------------------------------
1404 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
1405                                    TObjArray *twoTrackArray,AliVEvent *event,
1406                                    AliAODVertex *secVert,Double_t dca,
1407                                    Bool_t &okD0,Bool_t &okJPSI,
1408                                    Bool_t &okD0fromDstar) const
1409 {
1410   // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
1411   // reconstruction cuts
1412   // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
1413
1414   okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
1415
1416   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1417   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1418   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
1419
1420   // propagate tracks to secondary vertex, to compute inv. mass
1421   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1422   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1423
1424   Double_t momentum[3];
1425   postrack->GetPxPyPz(momentum);
1426   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1427   negtrack->GetPxPyPz(momentum);
1428   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1429
1430
1431   // invariant mass cut (try to improve coding here..)
1432   Bool_t okMassCut=kFALSE;
1433   if(!okMassCut && fD0toKpi)   if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
1434   if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
1435   if(!okMassCut && fDstar)     if(SelectInvMass(3,2,px,py,pz)) okMassCut=kTRUE;
1436   if(!okMassCut) {
1437     AliDebug(2," candidate didn't pass mass cut");
1438     return 0x0;
1439   }
1440   // primary vertex to be used by this candidate
1441   AliAODVertex *primVertexAOD  = PrimaryVertex(twoTrackArray,event);
1442   if(!primVertexAOD) return 0x0;
1443
1444
1445   Double_t d0z0[2],covd0z0[3];
1446   postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1447   d0[0] = d0z0[0];
1448   d0err[0] = TMath::Sqrt(covd0z0[0]);
1449   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1450   d0[1] = d0z0[0];
1451   d0err[1] = TMath::Sqrt(covd0z0[0]);
1452
1453   // create the object AliAODRecoDecayHF2Prong
1454   AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
1455   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
1456   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1457   the2Prong->SetProngIDs(2,id);
1458   delete primVertexAOD; primVertexAOD=NULL;
1459
1460
1461   if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
1462     // select D0->Kpi
1463     if(fD0toKpi)   {
1464       okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1465       if(okD0) the2Prong->SetSelectionBit(AliRDHFCuts::kD0toKpiCuts);
1466     }
1467     //if(fDebug && fD0toKpi) printf("   %d\n",(Int_t)okD0);
1468     // select J/psi from B
1469     if(fJPSItoEle)   {
1470       okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1471     }
1472     //if(fDebug && fJPSItoEle) printf("   %d\n",(Int_t)okJPSI);
1473     // select D0->Kpi from Dstar
1474     if(fDstar)   {
1475       okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
1476       if(okD0fromDstar) the2Prong->SetSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
1477     }
1478     //if(fDebug && fDstar) printf("   %d\n",(Int_t)okD0fromDstar);
1479   }
1480
1481
1482   // remove the primary vertex (was used only for selection)
1483   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1484     the2Prong->UnsetOwnPrimaryVtx();
1485   }
1486
1487   // get PID info from ESD
1488   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1489   if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1490   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1491   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1492   Double_t esdpid[10];
1493   for(Int_t i=0;i<5;i++) {
1494     esdpid[i]   = esdpid0[i];
1495     esdpid[5+i] = esdpid1[i];
1496   }
1497   the2Prong->SetPID(2,esdpid);
1498
1499   return the2Prong;
1500 }
1501 //----------------------------------------------------------------------------
1502 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
1503                              TObjArray *threeTrackArray,AliVEvent *event,
1504                              AliAODVertex *secVert,Double_t dispersion,
1505                              const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
1506                              Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1507                              Bool_t &ok3Prong) const
1508 {
1509   // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1510   // reconstruction cuts
1511   // E.Bruna, F.Prino
1512
1513
1514   ok3Prong=kFALSE;
1515   if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
1516
1517   Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1518   Double_t momentum[3];
1519
1520
1521   AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1522   AliESDtrack *negtrack  = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1523   AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1524
1525   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1526   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1527   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1528   postrack1->GetPxPyPz(momentum);
1529   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1530   negtrack->GetPxPyPz(momentum);
1531   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1532   postrack2->GetPxPyPz(momentum);
1533   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1534
1535   // invariant mass cut for D+, Ds, Lc
1536   Bool_t okMassCut=kFALSE;
1537   if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1538   if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
1539   if(!okMassCut) {
1540     AliDebug(2," candidate didn't pass mass cut");
1541     return 0x0;
1542   }
1543
1544   // primary vertex to be used by this candidate
1545   AliAODVertex *primVertexAOD  = PrimaryVertex(threeTrackArray,event);
1546   if(!primVertexAOD) return 0x0;
1547
1548   Double_t d0z0[2],covd0z0[3];
1549   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1550   d0[0]=d0z0[0];
1551   d0err[0] = TMath::Sqrt(covd0z0[0]);
1552   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1553   d0[1]=d0z0[0];
1554   d0err[1] = TMath::Sqrt(covd0z0[0]);
1555   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1556   d0[2]=d0z0[0];
1557   d0err[2] = TMath::Sqrt(covd0z0[0]);
1558
1559
1560   // create the object AliAODRecoDecayHF3Prong
1561   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1562   Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1563   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]));
1564   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]));
1565   Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
1566   AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1567   the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1568   UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1569   the3Prong->SetProngIDs(3,id);
1570
1571   delete primVertexAOD; primVertexAOD=NULL;
1572
1573   // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1574   if(f3Prong) {
1575     ok3Prong = kFALSE;
1576
1577     if(fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) {
1578       ok3Prong = kTRUE;
1579       the3Prong->SetSelectionBit(AliRDHFCuts::kDplusCuts);
1580     }
1581     if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) {
1582       ok3Prong = kTRUE;
1583       the3Prong->SetSelectionBit(AliRDHFCuts::kDsCuts);
1584     }
1585     if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) {
1586       ok3Prong = kTRUE;
1587       the3Prong->SetSelectionBit(AliRDHFCuts::kLcCuts);
1588     }
1589   }
1590   //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1591
1592   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1593     the3Prong->UnsetOwnPrimaryVtx();
1594   }
1595
1596   // get PID info from ESD
1597   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1598   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1599   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1600   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1601   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1602   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1603
1604   Double_t esdpid[15];
1605   for(Int_t i=0;i<5;i++) {
1606     esdpid[i]    = esdpid0[i];
1607     esdpid[5+i]  = esdpid1[i];
1608     esdpid[10+i] = esdpid2[i];
1609   }
1610   the3Prong->SetPID(3,esdpid);
1611
1612   return the3Prong;
1613 }
1614 //----------------------------------------------------------------------------
1615 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1616                              TObjArray *fourTrackArray,AliVEvent *event,
1617                              AliAODVertex *secVert,
1618                              const AliAODVertex *vertexp1n1,
1619                              const AliAODVertex *vertexp1n1p2,
1620                              Double_t dcap1n1,Double_t dcap1n2,
1621                              Double_t dcap2n1,Double_t dcap2n2,
1622                              Bool_t &ok4Prong) const
1623 {
1624   // Make 4Prong candidates and check if they pass D0toKpipipi
1625   // reconstruction cuts
1626   // G.E.Bruno, R.Romita
1627
1628   ok4Prong=kFALSE;
1629   if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
1630
1631   Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1632
1633   AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1634   AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1635   AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1636   AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1637
1638   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1639   negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1640   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1641   negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1642
1643   Double_t momentum[3];
1644   postrack1->GetPxPyPz(momentum);
1645   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1646   negtrack1->GetPxPyPz(momentum);
1647   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1648   postrack2->GetPxPyPz(momentum);
1649   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1650   negtrack2->GetPxPyPz(momentum);
1651   px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1652
1653   // invariant mass cut for rho or D0 (try to improve coding here..)
1654   Bool_t okMassCut=kFALSE;
1655   if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1656   if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) {      //no PID, to be implemented with PID
1657     if(SelectInvMass(4,4,px,py,pz)) okMassCut=kTRUE;
1658   }
1659   if(!okMassCut) {
1660     //if(fDebug) printf(" candidate didn't pass mass cut\n");
1661     //printf(" candidate didn't pass mass cut\n");
1662     return 0x0;
1663   }
1664
1665   // primary vertex to be used by this candidate
1666   AliAODVertex *primVertexAOD  = PrimaryVertex(fourTrackArray,event);
1667   if(!primVertexAOD) return 0x0;
1668
1669   Double_t d0z0[2],covd0z0[3];
1670   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1671   d0[0]=d0z0[0];
1672   d0err[0] = TMath::Sqrt(covd0z0[0]);
1673   negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1674   d0[1]=d0z0[0];
1675   d0err[1] = TMath::Sqrt(covd0z0[0]);
1676   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1677   d0[2]=d0z0[0];
1678   d0err[2] = TMath::Sqrt(covd0z0[0]);
1679   negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1680   d0[3]=d0z0[0];
1681   d0err[3] = TMath::Sqrt(covd0z0[0]);
1682
1683
1684   // create the object AliAODRecoDecayHF4Prong
1685   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1686   Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
1687   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]));
1688   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]));
1689   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]));
1690   Short_t charge=0;
1691   AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
1692   the4Prong->SetOwnPrimaryVtx(primVertexAOD);
1693   UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
1694   the4Prong->SetProngIDs(4,id);
1695
1696   delete primVertexAOD; primVertexAOD=NULL;
1697
1698   ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
1699
1700
1701   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1702     the4Prong->UnsetOwnPrimaryVtx();
1703   }
1704
1705
1706   // get PID info from ESD
1707   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1708   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1709   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1710   if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
1711   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1712   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1713   Double_t esdpid3[5]={0.,0.,0.,0.,0.};
1714   if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
1715
1716   Double_t esdpid[20];
1717   for(Int_t i=0;i<5;i++) {
1718     esdpid[i]    = esdpid0[i];
1719     esdpid[5+i]  = esdpid1[i];
1720     esdpid[10+i] = esdpid2[i];
1721     esdpid[15+i] = esdpid3[i];
1722   }
1723   the4Prong->SetPID(4,esdpid);
1724
1725   return the4Prong;
1726 }
1727 //-----------------------------------------------------------------------------
1728 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
1729                                                     AliVEvent *event) const
1730 {
1731   // Returns primary vertex to be used for this candidate
1732
1733   AliESDVertex *vertexESD = 0;
1734   AliAODVertex *vertexAOD = 0;
1735
1736
1737   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1738     // primary vertex from the input event
1739
1740     vertexESD = new AliESDVertex(*fV1);
1741
1742   } else {
1743     // primary vertex specific to this candidate
1744
1745     Int_t nTrks = trkArray->GetEntriesFast();
1746     AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1747
1748     if(fRecoPrimVtxSkippingTrks) {
1749       // recalculating the vertex
1750
1751       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1752         Float_t diamondcovxy[3];
1753         event->GetDiamondCovXY(diamondcovxy);
1754         Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1755         Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1756         AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1757         vertexer->SetVtxStart(diamond);
1758         delete diamond; diamond=NULL;
1759         if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1760           vertexer->SetOnlyFitter();
1761       }
1762       Int_t skipped[1000];
1763       Int_t nTrksToSkip=0,id;
1764       AliExternalTrackParam *t = 0;
1765       for(Int_t i=0; i<nTrks; i++) {
1766         t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1767         id = (Int_t)t->GetID();
1768         if(id<0) continue;
1769         skipped[nTrksToSkip++] = id;
1770       }
1771       // TEMPORARY FIX
1772       // For AOD, skip also tracks without covariance matrix
1773       if(fInputAOD) {
1774         Double_t covtest[21];
1775         for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1776           AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1777           if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1778             id = (Int_t)vtrack->GetID();
1779             if(id<0) continue;
1780             skipped[nTrksToSkip++] = id;
1781           }
1782         }
1783       }
1784       //
1785       vertexer->SetSkipTracks(nTrksToSkip,skipped);
1786       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1787
1788     } else if(fRmTrksFromPrimVtx) {
1789       // removing the prongs tracks
1790
1791       TObjArray rmArray(nTrks);
1792       UShort_t *rmId = new UShort_t[nTrks];
1793       AliESDtrack *esdTrack = 0;
1794       AliESDtrack *t = 0;
1795       for(Int_t i=0; i<nTrks; i++) {
1796         t = (AliESDtrack*)trkArray->UncheckedAt(i);
1797         esdTrack = new AliESDtrack(*t);
1798         rmArray.AddLast(esdTrack);
1799         if(esdTrack->GetID()>=0) {
1800           rmId[i]=(UShort_t)esdTrack->GetID();
1801         } else {
1802           rmId[i]=9999;
1803         }
1804       }
1805       Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
1806       vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1807       delete [] rmId; rmId=NULL;
1808       rmArray.Delete();
1809
1810     }
1811
1812     if(!vertexESD) return vertexAOD;
1813     if(vertexESD->GetNContributors()<=0) {
1814       AliDebug(2,"vertexing failed");
1815       delete vertexESD; vertexESD=NULL;
1816       return vertexAOD;
1817     }
1818
1819     delete vertexer; vertexer=NULL;
1820
1821   }
1822
1823   // convert to AliAODVertex
1824   Double_t pos[3],cov[6],chi2perNDF;
1825   vertexESD->GetXYZ(pos); // position
1826   vertexESD->GetCovMatrix(cov); //covariance matrix
1827   chi2perNDF = vertexESD->GetChi2toNDF();
1828   delete vertexESD; vertexESD=NULL;
1829
1830   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1831
1832   return vertexAOD;
1833 }
1834 //-----------------------------------------------------------------------------
1835 void AliAnalysisVertexingHF::PrintStatus() const {
1836   // Print parameters being used
1837
1838   //printf("Preselections:\n");
1839   //   fTrackFilter->Dump();
1840   if(fSecVtxWithKF) {
1841     printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1842   } else {
1843     printf("Secondary vertex with AliVertexerTracks\n");
1844   }
1845   if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1846   if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1847   if(fD0toKpi) {
1848     printf("Reconstruct D0->Kpi candidates with cuts:\n");
1849     if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
1850   }
1851   if(fDstar) {
1852     printf("Reconstruct D*->D0pi candidates with cuts:\n");
1853     if(fFindVertexForDstar) {
1854       printf("    Reconstruct a secondary vertex for the D*\n");
1855     } else {
1856       printf("    Assume the D* comes from the primary vertex\n");
1857     }
1858     if(fCutsDStartoKpipi) fCutsDStartoKpipi->PrintAll();
1859   }
1860   if(fJPSItoEle) {
1861     printf("Reconstruct J/psi from B candidates with cuts:\n");
1862     if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
1863   }
1864   if(f3Prong) {
1865     printf("Reconstruct 3 prong candidates.\n");
1866     printf("  D+->Kpipi cuts:\n");
1867     if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
1868     printf("  Ds->KKpi cuts:\n");
1869     if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
1870     printf("  Lc->pKpi cuts:\n");
1871     if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
1872   }
1873   if(f4Prong) {
1874     printf("Reconstruct 4 prong candidates.\n");
1875     printf("  D0->Kpipipi cuts:\n");
1876     if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
1877   }
1878   if(fCascades) {
1879     printf("Reconstruct cascades candidates formed with v0s.\n");
1880     printf("  Lc -> k0s P & Lc -> L Pi cuts:\n");
1881     if(fCutsLctoV0) fCutsLctoV0->PrintAll();
1882   }
1883
1884   return;
1885 }
1886 //-----------------------------------------------------------------------------
1887 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
1888                                                                  Double_t &dispersion,Bool_t useTRefArray) const
1889 {
1890   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1891
1892   AliESDVertex *vertexESD = 0;
1893   AliAODVertex *vertexAOD = 0;
1894
1895   if(!fSecVtxWithKF) { // AliVertexerTracks
1896
1897     AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1898     vertexer->SetVtxStart(fV1);
1899     vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1900     delete vertexer; vertexer=NULL;
1901
1902     if(!vertexESD) return vertexAOD;
1903
1904     if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
1905       AliDebug(2,"vertexing failed");
1906       delete vertexESD; vertexESD=NULL;
1907       return vertexAOD;
1908     }
1909
1910   } else { // Kalman Filter vertexer (AliKFParticle)
1911
1912     AliKFParticle::SetField(fBzkG);
1913
1914     AliKFVertex vertexKF;
1915
1916     Int_t nTrks = trkArray->GetEntriesFast();
1917     for(Int_t i=0; i<nTrks; i++) {
1918       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1919       AliKFParticle daughterKF(*esdTrack,211);
1920       vertexKF.AddDaughter(daughterKF);
1921     }
1922     vertexESD = new AliESDVertex(vertexKF.Parameters(),
1923                                  vertexKF.CovarianceMatrix(),
1924                                  vertexKF.GetChi2(),
1925                                  vertexKF.GetNContributors());
1926
1927   }
1928
1929   // convert to AliAODVertex
1930   Double_t pos[3],cov[6],chi2perNDF;
1931   vertexESD->GetXYZ(pos); // position
1932   vertexESD->GetCovMatrix(cov); //covariance matrix
1933   chi2perNDF = vertexESD->GetChi2toNDF();
1934   dispersion = vertexESD->GetDispersion();
1935   delete vertexESD; vertexESD=NULL;
1936
1937   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
1938   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
1939
1940   return vertexAOD;
1941 }
1942 //-----------------------------------------------------------------------------
1943 Bool_t AliAnalysisVertexingHF::SelectInvMass(TObjArray *trkArray) const {
1944   // Invariant mass cut on tracks
1945
1946   Int_t nProngs=trkArray->GetEntriesFast();
1947   Int_t retval=kFALSE;
1948
1949   Double_t momentum[3];
1950   Double_t px3[3],py3[3],pz3[3];
1951   Double_t px4[4],py4[4],pz4[4];
1952   AliESDtrack *postrack1=0;
1953   AliESDtrack *postrack2=0;
1954   AliESDtrack *negtrack1=0;
1955   AliESDtrack *negtrack2=0;
1956
1957   switch(nProngs) {
1958   case 3:
1959     // invariant mass cut for D+, Ds, Lc
1960     postrack1 = (AliESDtrack*)trkArray->UncheckedAt(0);
1961     negtrack1  = (AliESDtrack*)trkArray->UncheckedAt(1);
1962     postrack2 = (AliESDtrack*)trkArray->UncheckedAt(2);
1963     postrack1->GetPxPyPz(momentum);
1964     px3[0] = momentum[0]; py3[0] = momentum[1]; pz3[0] = momentum[2];
1965     negtrack1->GetPxPyPz(momentum);
1966     px3[1] = momentum[0]; py3[1] = momentum[1]; pz3[1] = momentum[2];
1967     postrack2->GetPxPyPz(momentum);
1968     px3[2] = momentum[0]; py3[2] = momentum[1]; pz3[2] = momentum[2];
1969     retval = SelectInvMass(2,3,px3,py3,pz3);
1970     break;
1971   case 4:
1972     // invariant mass cut for D0->4prong
1973     postrack1 = (AliESDtrack*)trkArray->UncheckedAt(0);
1974     negtrack1 = (AliESDtrack*)trkArray->UncheckedAt(1);
1975     postrack2 = (AliESDtrack*)trkArray->UncheckedAt(2);
1976     negtrack2 = (AliESDtrack*)trkArray->UncheckedAt(3);
1977     postrack1->GetPxPyPz(momentum);
1978     px4[0] = momentum[0]; py4[0] = momentum[1]; pz4[0] = momentum[2];
1979     negtrack1->GetPxPyPz(momentum);
1980     px4[1] = momentum[0]; py4[1] = momentum[1]; pz4[1] = momentum[2];
1981     postrack2->GetPxPyPz(momentum);
1982     px4[2] = momentum[0]; py4[2] = momentum[1]; pz4[2] = momentum[2];
1983     negtrack2->GetPxPyPz(momentum);
1984     px4[3] = momentum[0]; py4[3] = momentum[1]; pz4[3] = momentum[2];
1985     retval = SelectInvMass(4,4,px4,py4,pz4);
1986     break;
1987   default:
1988     printf("SelectInvMass(): wrong decay selection\n");
1989     break;
1990   }
1991
1992   return retval;
1993 }
1994 //-----------------------------------------------------------------------------
1995 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1996                                              Int_t nprongs,
1997                                              Double_t *px,
1998                                              Double_t *py,
1999                                              Double_t *pz) const {
2000   // Check invariant mass cut
2001
2002   UInt_t pdg2[2],pdg3[3],pdg4[4];
2003   Double_t mPDG,minv;
2004
2005   Bool_t retval=kFALSE;
2006   switch (decay)
2007     {
2008     case 0:                  // D0->Kpi
2009       fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2010       pdg2[0]=211; pdg2[1]=321;
2011       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2012       minv = fMassCalc2->InvMass(nprongs,pdg2);
2013       if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
2014       pdg2[0]=321; pdg2[1]=211;
2015       minv = fMassCalc2->InvMass(nprongs,pdg2);
2016       if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
2017       break;
2018     case 1:                  // JPSI->ee
2019       fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2020       pdg2[0]=11; pdg2[1]=11;
2021       mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
2022       minv = fMassCalc2->InvMass(nprongs,pdg2);
2023       if(TMath::Abs(minv-mPDG)<fCutsJpsitoee->GetMassCut()) retval=kTRUE;
2024       break;
2025     case 2:                  // D+->Kpipi
2026       fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
2027       pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
2028       mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
2029       minv = fMassCalc3->InvMass(nprongs,pdg3);
2030       if(TMath::Abs(minv-mPDG)<fCutsDplustoKpipi->GetMassCut()) retval=kTRUE;
2031                             // Ds+->KKpi
2032       pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
2033       mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
2034       minv = fMassCalc3->InvMass(nprongs,pdg3);
2035       if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
2036       pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
2037       minv = fMassCalc3->InvMass(nprongs,pdg3);
2038       if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
2039                             // Lc->pKpi
2040       pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
2041       mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
2042       minv = fMassCalc3->InvMass(nprongs,pdg3);
2043       if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
2044       pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
2045       minv = fMassCalc3->InvMass(nprongs,pdg3);
2046       if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
2047       break;
2048     case 3:                  // D*->D0pi
2049       fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2050       pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
2051       mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
2052       minv = fMassCalc2->InvMass(nprongs,pdg2);
2053       if(TMath::Abs(minv-mPDG)<fCutsDStartoKpipi->GetMassCut()) retval=kTRUE;
2054       break;
2055     case 4:                 // D0->Kpipipi without PID
2056       fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
2057       pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
2058       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2059       minv = fMassCalc4->InvMass(nprongs,pdg4);
2060       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
2061       pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
2062       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2063       minv = fMassCalc4->InvMass(nprongs,pdg4);
2064       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
2065       pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
2066       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2067       minv = fMassCalc4->InvMass(nprongs,pdg4);
2068       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
2069       pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
2070       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2071       minv = fMassCalc4->InvMass(nprongs,pdg4);
2072       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
2073       break;
2074     default:
2075       printf("SelectInvMass(): wrong decay selection\n");
2076       break;
2077     }
2078
2079   return retval;
2080 }
2081 //-----------------------------------------------------------------------------
2082 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
2083                                                        Int_t trkEntries,
2084                                    TObjArray &seleTrksArray,Int_t &nSeleTrks,
2085                                    UChar_t *seleFlags,Int_t *evtNumber)
2086 {
2087   // Apply single-track preselection.
2088   // Fill a TObjArray with selected tracks (for displaced vertices or
2089   // soft pion from D*). Selection flag stored in seleFlags.
2090   // Create the AliESDVertex object (convert from AliAODVertex if necessary)
2091   // In case of AOD input, also fill fAODMap for track index<->ID
2092
2093   const AliVVertex *vprimary = event->GetPrimaryVertex();
2094
2095   if(fV1) { delete fV1; fV1=NULL; }
2096   if(fAODMap) { delete fAODMap; fAODMap=NULL; }
2097
2098   Int_t nindices=0;
2099   UShort_t *indices = 0;
2100   Double_t pos[3],cov[6];
2101
2102   if(!fInputAOD) { // ESD
2103     fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
2104   } else {         // AOD
2105     vprimary->GetXYZ(pos);
2106     vprimary->GetCovarianceMatrix(cov);
2107     fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2108     indices = new UShort_t[event->GetNumberOfTracks()];
2109     fAODMapSize = 100000;
2110     fAODMap = new Int_t[fAODMapSize];
2111   }
2112
2113
2114   Int_t entries = (Int_t)event->GetNumberOfTracks();
2115   Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
2116   nSeleTrks=0;
2117
2118   // transfer ITS tracks from event to arrays
2119   for(Int_t i=0; i<entries; i++) {
2120     AliVTrack *track;
2121     track = (AliVTrack*)event->GetTrack(i);
2122
2123     // skip pure ITS SA tracks
2124     if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2125
2126     // skip tracks without ITS
2127     if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
2128
2129     // TEMPORARY: check that the cov matrix is there
2130     Double_t covtest[21];
2131     if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2132     //
2133
2134     if(fInputAOD) {
2135       AliAODTrack *aodt = (AliAODTrack*)track;
2136       if(aodt->GetUsedForPrimVtxFit()) {
2137         indices[nindices]=aodt->GetID(); nindices++;
2138       }
2139       fAODMap[(Int_t)aodt->GetID()] = i;
2140     }
2141
2142     AliESDtrack *esdt = 0;
2143
2144     if(!fInputAOD) {
2145       esdt = (AliESDtrack*)track;
2146     } else {
2147       esdt = new AliESDtrack(track);
2148     }
2149
2150     // single track selection
2151     okDisplaced=kFALSE; okSoftPi=kFALSE;
2152     if(fMixEvent && i<trkEntries){
2153       evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
2154       const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
2155       Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
2156       eventVtx->GetXYZ(vtxPos);
2157       vprimary->GetXYZ(primPos);
2158       eventVtx->GetCovarianceMatrix(primCov);
2159       for(Int_t ind=0;ind<3;ind++){
2160         trasl[ind]=vtxPos[ind]-primPos[ind];
2161       }
2162
2163       Bool_t isTransl=esdt->Translate(trasl,primCov);
2164       if(!isTransl) {
2165         delete esdt;
2166         esdt = NULL;
2167         continue;
2168       }
2169     }
2170
2171     if(SingleTrkCuts(esdt,okDisplaced,okSoftPi) && nSeleTrks<trkEntries) {
2172       seleTrksArray.AddLast(esdt);
2173       seleFlags[nSeleTrks]=0;
2174       if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2175       if(okSoftPi)    SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2176       nSeleTrks++;
2177     } else {
2178       if(fInputAOD) delete esdt;
2179       esdt = NULL;
2180       continue;
2181     }
2182
2183   } // end loop on tracks
2184
2185   // primary vertex from AOD
2186   if(fInputAOD) {
2187     delete fV1; fV1=NULL;
2188     vprimary->GetXYZ(pos);
2189     vprimary->GetCovarianceMatrix(cov);
2190     Double_t chi2toNDF = vprimary->GetChi2perNDF();
2191     Int_t ncontr=nindices;
2192     if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
2193     Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
2194     fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2195     fV1->SetTitle(vprimary->GetTitle());
2196     fV1->SetIndices(nindices,indices);
2197   }
2198   if(indices) { delete [] indices; indices=NULL; }
2199
2200
2201   return;
2202 }
2203 //-----------------------------------------------------------------------------
2204 void AliAnalysisVertexingHF::SetSelectionBitForPID(AliRDHFCuts *cuts,AliAODRecoDecayHF *rd) {
2205   //
2206   // Set the selection bit for PID
2207   //
2208   if(cuts->GetPidHF()) {
2209     Bool_t usepid=cuts->GetIsUsePID();
2210     cuts->SetUsePID(kTRUE);
2211     if(cuts->IsSelectedPID(rd))
2212       rd->SetSelectionBit(AliRDHFCuts::kDstarPID);
2213     cuts->SetUsePID(usepid);
2214   }
2215   return;
2216 }
2217 //-----------------------------------------------------------------------------
2218 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2219                                              Bool_t &okDisplaced,Bool_t &okSoftPi) const
2220 {
2221   // Check if track passes some kinematical cuts
2222
2223   // this is needed to store the impact parameters
2224   trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2225
2226   UInt_t selectInfo;
2227   //
2228   // Track selection, displaced tracks
2229   selectInfo = 0;
2230   if(fTrackFilter) {
2231     selectInfo = fTrackFilter->IsSelected(trk);
2232   }
2233   if(selectInfo) okDisplaced=kTRUE;
2234
2235   // Track selection, soft pions
2236   selectInfo = 0;
2237   if(fDstar && fTrackFilterSoftPi) {
2238     selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2239   }
2240   if(selectInfo) okSoftPi=kTRUE;
2241
2242   if(okDisplaced || okSoftPi) return kTRUE;
2243
2244   return kFALSE;
2245 }
2246
2247
2248 //-----------------------------------------------------------------------------
2249 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2250   //
2251   // Transform ESDv0 to AODv0
2252   //
2253   //  this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2254   //  and creates an AODv0 out of them
2255   //
2256   Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
2257   AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2258
2259   // create the v0 neutral track to compute the DCA to the primary vertex
2260   Double_t xyz[3], pxpypz[3];
2261   esdV0->XvYvZv(xyz);
2262   esdV0->PxPyPz(pxpypz);
2263   Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2264   AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2265   if(!trackesdV0) {
2266     delete vertexV0;
2267     return 0;
2268   }
2269   Double_t d0z0[2],covd0z0[3];
2270   AliAODVertex *primVertexAOD = PrimaryVertex();
2271   trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2272   Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2273   // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2274   Double_t dcaV0DaughterToPrimVertex[2];
2275   AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2276   AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
2277   if( !posV0track || !negV0track) {
2278     if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2279     delete vertexV0;
2280     delete primVertexAOD;
2281     return 0;
2282   }
2283   posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2284   //  if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2285   //  else
2286   dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
2287   negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2288   //  if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2289   //  else
2290   dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
2291   Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
2292   Double_t pmom[3],nmom[3];
2293   esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2294   esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2295
2296   AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2297   aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2298
2299   delete trackesdV0;
2300   delete primVertexAOD;
2301
2302   return aodV0;
2303 }
2304 //-----------------------------------------------------------------------------