PWG3/vertexingHF/AliAnalysisVertexingHF.cxx

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