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