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,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 exit
 609     if(trkEntries<2) {
 610       AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
 611       return;
 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 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1208                                    TObjArray *twoTrackArray,AliVEvent *event,
1209                                    AliAODVertex *secVert,
1210                                    AliAODRecoDecayHF2Prong *rd2Prong,
1211                                    Double_t dca,
1212                                    Bool_t &okDstar) const
1213 {
1214   // Make the cascade as a 2Prong decay and check if it passes Dstar
1215   // reconstruction cuts
1216
1217   okDstar = kFALSE;
1218
1219   Bool_t dummy1,dummy2,dummy3;
1220
1221   // We use Make2Prong to construct the AliAODRecoCascadeHF
1222   // (which inherits from AliAODRecoDecayHF2Prong)
1223   AliAODRecoCascadeHF *theCascade =
1224     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1225                                      dummy1,dummy2,dummy3);
1226   if(!theCascade) return 0x0;
1227
1228   // charge
1229   AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1230   theCascade->SetCharge(trackPi->Charge());
1231
1232   //--- selection cuts
1233   //
1234   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1235   tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
1236   tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
1237   AliAODVertex *primVertexAOD=0;
1238   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1239     // take event primary vertex
1240     primVertexAOD = PrimaryVertex();
1241     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1242     rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
1243   }
1244   // select D*->D0pi
1245   if(fDstar) {
1246     okDstar = (Bool_t)fCutsDStartoKpipi->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1247   }
1248   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1249   tmpCascade->UnsetOwnPrimaryVtx();
1250   delete tmpCascade; tmpCascade=NULL;
1251   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1252     rd2Prong->UnsetOwnPrimaryVtx();
1253   }
1254   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1255   //---
1256
1257   return theCascade;
1258 }
1259
1260
1261 //----------------------------------------------------------------------------
1262 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1263                                    TObjArray *twoTrackArray,AliVEvent *event,
1264                                    AliAODVertex *secVert,
1265                                    AliAODv0 *v0,
1266                                    Double_t dca,
1267                                    Bool_t &okCascades) const
1268 {
1269   //
1270   // Make the cascade as a 2Prong decay and check if it passes
1271   // cascades reconstruction cuts
1272
1273   //  AliDebug(2,Form("         building the cascade"));
1274   okCascades= kFALSE;
1275   Bool_t dummy1,dummy2,dummy3;
1276
1277   // We use Make2Prong to construct the AliAODRecoCascadeHF
1278   // (which inherits from AliAODRecoDecayHF2Prong)
1279   AliAODRecoCascadeHF *theCascade =
1280     (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1281                                      dummy1,dummy2,dummy3);
1282   if(!theCascade) return 0x0;
1283
1284   // bachelor track and charge
1285   AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1286   theCascade->SetCharge(trackBachelor->Charge());
1287
1288   //--- selection cuts
1289   //
1290   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1291   tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
1292   tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
1293   AliAODVertex *primVertexAOD=0;
1294   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1295     // take event primary vertex
1296     primVertexAOD = PrimaryVertex();
1297     if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
1298     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1299   }
1300
1301   // select Cascades
1302   Bool_t okLcksp=0, okLcLpi=0;
1303   if(fCascades && fInputAOD){
1304     okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1305     if(okCascades==1) okLcksp=1;
1306     if(okCascades==2) okLcLpi=1;
1307     if(okCascades==3) { okLcksp=1; okLcLpi=1;}
1308   }
1309   else { AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied")); okCascades=kTRUE; }// no cuts implemented from ESDs
1310   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1311   tmpCascade->UnsetOwnPrimaryVtx();
1312   delete tmpCascade; tmpCascade=NULL;
1313   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1314   //---
1315
1316   return theCascade;
1317 }
1318
1319 //-----------------------------------------------------------------------------
1320 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
1321                                    TObjArray *twoTrackArray,AliVEvent *event,
1322                                    AliAODVertex *secVert,Double_t dca,
1323                                    Bool_t &okD0,Bool_t &okJPSI,
1324                                    Bool_t &okD0fromDstar) const
1325 {
1326   // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
1327   // reconstruction cuts
1328   // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
1329
1330   okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
1331
1332   Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1333   AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1334   AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
1335
1336   // propagate tracks to secondary vertex, to compute inv. mass
1337   postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1338   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1339
1340   Double_t momentum[3];
1341   postrack->GetPxPyPz(momentum);
1342   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1343   negtrack->GetPxPyPz(momentum);
1344   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1345
1346
1347   // invariant mass cut (try to improve coding here..)
1348   Bool_t okMassCut=kFALSE;
1349   if(!okMassCut && fD0toKpi)   if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
1350   if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
1351   if(!okMassCut && fDstar)     if(SelectInvMass(3,2,px,py,pz)) okMassCut=kTRUE;
1352   if(!okMassCut) {
1353     AliDebug(2," candidate didn't pass mass cut");
1354     return 0x0;
1355   }
1356   // primary vertex to be used by this candidate
1357   AliAODVertex *primVertexAOD  = PrimaryVertex(twoTrackArray,event);
1358   if(!primVertexAOD) return 0x0;
1359
1360
1361   Double_t d0z0[2],covd0z0[3];
1362   postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1363   d0[0] = d0z0[0];
1364   d0err[0] = TMath::Sqrt(covd0z0[0]);
1365   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1366   d0[1] = d0z0[0];
1367   d0err[1] = TMath::Sqrt(covd0z0[0]);
1368
1369   // create the object AliAODRecoDecayHF2Prong
1370   AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
1371   the2Prong->SetOwnPrimaryVtx(primVertexAOD);
1372   UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1373   the2Prong->SetProngIDs(2,id);
1374   delete primVertexAOD; primVertexAOD=NULL;
1375
1376
1377   if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
1378     // select D0->Kpi
1379     if(fD0toKpi)   okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1380     //if(fDebug && fD0toKpi) printf("   %d\n",(Int_t)okD0);
1381     // select J/psi from B
1382     if(fJPSItoEle)   okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1383     //if(fDebug && fJPSItoEle) printf("   %d\n",(Int_t)okJPSI);
1384     // select D0->Kpi from Dstar
1385     if(fDstar)   okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
1386     //if(fDebug && fDstar) printf("   %d\n",(Int_t)okD0fromDstar);
1387   }
1388
1389
1390   // remove the primary vertex (was used only for selection)
1391   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1392     the2Prong->UnsetOwnPrimaryVtx();
1393   }
1394
1395   // get PID info from ESD
1396   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1397   if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1398   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1399   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1400   Double_t esdpid[10];
1401   for(Int_t i=0;i<5;i++) {
1402     esdpid[i]   = esdpid0[i];
1403     esdpid[5+i] = esdpid1[i];
1404   }
1405   the2Prong->SetPID(2,esdpid);
1406
1407   return the2Prong;
1408 }
1409 //----------------------------------------------------------------------------
1410 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
1411                              TObjArray *threeTrackArray,AliVEvent *event,
1412                              AliAODVertex *secVert,Double_t dispersion,
1413                              const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
1414                              Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1415                              Bool_t &ok3Prong) const
1416 {
1417   // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1418   // reconstruction cuts
1419   // E.Bruna, F.Prino
1420
1421
1422   ok3Prong=kFALSE;
1423   if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
1424
1425   Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1426   Double_t momentum[3];
1427
1428
1429   AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1430   AliESDtrack *negtrack  = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1431   AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1432
1433   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1434   negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1435   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1436   postrack1->GetPxPyPz(momentum);
1437   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1438   negtrack->GetPxPyPz(momentum);
1439   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1440   postrack2->GetPxPyPz(momentum);
1441   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1442
1443   // invariant mass cut for D+, Ds, Lc
1444   Bool_t okMassCut=kFALSE;
1445   if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1446   if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
1447   if(!okMassCut) {
1448     AliDebug(2," candidate didn't pass mass cut");
1449     return 0x0;
1450   }
1451
1452   // primary vertex to be used by this candidate
1453   AliAODVertex *primVertexAOD  = PrimaryVertex(threeTrackArray,event);
1454   if(!primVertexAOD) return 0x0;
1455
1456   Double_t d0z0[2],covd0z0[3];
1457   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1458   d0[0]=d0z0[0];
1459   d0err[0] = TMath::Sqrt(covd0z0[0]);
1460   negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1461   d0[1]=d0z0[0];
1462   d0err[1] = TMath::Sqrt(covd0z0[0]);
1463   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1464   d0[2]=d0z0[0];
1465   d0err[2] = TMath::Sqrt(covd0z0[0]);
1466
1467
1468   // create the object AliAODRecoDecayHF3Prong
1469   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1470   Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1471   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]));
1472   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]));
1473   Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
1474   AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1475   the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1476   UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1477   the3Prong->SetProngIDs(3,id);
1478
1479   delete primVertexAOD; primVertexAOD=NULL;
1480
1481   // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1482   if(f3Prong) {
1483     ok3Prong = kFALSE;
1484
1485     if(fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1486     if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1487     if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1488
1489   }
1490   //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1491
1492   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1493     the3Prong->UnsetOwnPrimaryVtx();
1494   }
1495
1496   // get PID info from ESD
1497   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1498   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1499   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1500   if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1501   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1502   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1503
1504   Double_t esdpid[15];
1505   for(Int_t i=0;i<5;i++) {
1506     esdpid[i]    = esdpid0[i];
1507     esdpid[5+i]  = esdpid1[i];
1508     esdpid[10+i] = esdpid2[i];
1509   }
1510   the3Prong->SetPID(3,esdpid);
1511
1512   return the3Prong;
1513 }
1514 //----------------------------------------------------------------------------
1515 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1516                              TObjArray *fourTrackArray,AliVEvent *event,
1517                              AliAODVertex *secVert,
1518                              const AliAODVertex *vertexp1n1,
1519                              const AliAODVertex *vertexp1n1p2,
1520                              Double_t dcap1n1,Double_t dcap1n2,
1521                              Double_t dcap2n1,Double_t dcap2n2,
1522                              Bool_t &ok4Prong) const
1523 {
1524   // Make 4Prong candidates and check if they pass D0toKpipipi
1525   // reconstruction cuts
1526   // G.E.Bruno, R.Romita
1527
1528   ok4Prong=kFALSE;
1529   if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
1530
1531   Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1532
1533   AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1534   AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1535   AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1536   AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1537
1538   postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1539   negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1540   postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1541   negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1542
1543   Double_t momentum[3];
1544   postrack1->GetPxPyPz(momentum);
1545   px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1546   negtrack1->GetPxPyPz(momentum);
1547   px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1548   postrack2->GetPxPyPz(momentum);
1549   px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1550   negtrack2->GetPxPyPz(momentum);
1551   px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1552
1553   // invariant mass cut for rho or D0 (try to improve coding here..)
1554   Bool_t okMassCut=kFALSE;
1555   if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1556   if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) {      //no PID, to be implemented with PID
1557     if(SelectInvMass(4,4,px,py,pz)) okMassCut=kTRUE;
1558   }
1559   if(!okMassCut) {
1560     //if(fDebug) printf(" candidate didn't pass mass cut\n");
1561     //printf(" candidate didn't pass mass cut\n");
1562     return 0x0;
1563   }
1564
1565   // primary vertex to be used by this candidate
1566   AliAODVertex *primVertexAOD  = PrimaryVertex(fourTrackArray,event);
1567   if(!primVertexAOD) return 0x0;
1568
1569   Double_t d0z0[2],covd0z0[3];
1570   postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1571   d0[0]=d0z0[0];
1572   d0err[0] = TMath::Sqrt(covd0z0[0]);
1573   negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1574   d0[1]=d0z0[0];
1575   d0err[1] = TMath::Sqrt(covd0z0[0]);
1576   postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1577   d0[2]=d0z0[0];
1578   d0err[2] = TMath::Sqrt(covd0z0[0]);
1579   negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1580   d0[3]=d0z0[0];
1581   d0err[3] = TMath::Sqrt(covd0z0[0]);
1582
1583
1584   // create the object AliAODRecoDecayHF4Prong
1585   Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1586   Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
1587   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]));
1588   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]));
1589   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]));
1590   Short_t charge=0;
1591   AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
1592   the4Prong->SetOwnPrimaryVtx(primVertexAOD);
1593   UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
1594   the4Prong->SetProngIDs(4,id);
1595
1596   delete primVertexAOD; primVertexAOD=NULL;
1597
1598   ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
1599
1600
1601   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1602     the4Prong->UnsetOwnPrimaryVtx();
1603   }
1604
1605
1606   // get PID info from ESD
1607   Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1608   if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1609   Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1610   if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
1611   Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1612   if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1613   Double_t esdpid3[5]={0.,0.,0.,0.,0.};
1614   if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
1615
1616   Double_t esdpid[20];
1617   for(Int_t i=0;i<5;i++) {
1618     esdpid[i]    = esdpid0[i];
1619     esdpid[5+i]  = esdpid1[i];
1620     esdpid[10+i] = esdpid2[i];
1621     esdpid[15+i] = esdpid3[i];
1622   }
1623   the4Prong->SetPID(4,esdpid);
1624
1625   return the4Prong;
1626 }
1627 //-----------------------------------------------------------------------------
1628 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
1629                                                     AliVEvent *event) const
1630 {
1631   // Returns primary vertex to be used for this candidate
1632
1633   AliESDVertex *vertexESD = 0;
1634   AliAODVertex *vertexAOD = 0;
1635
1636
1637   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1638     // primary vertex from the input event
1639
1640     vertexESD = new AliESDVertex(*fV1);
1641
1642   } else {
1643     // primary vertex specific to this candidate
1644
1645     Int_t nTrks = trkArray->GetEntriesFast();
1646     AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1647
1648     if(fRecoPrimVtxSkippingTrks) {
1649       // recalculating the vertex
1650
1651       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1652         Float_t diamondcovxy[3];
1653         event->GetDiamondCovXY(diamondcovxy);
1654         Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1655         Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1656         AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1657         vertexer->SetVtxStart(diamond);
1658         delete diamond; diamond=NULL;
1659         if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1660           vertexer->SetOnlyFitter();
1661       }
1662       Int_t skipped[1000];
1663       Int_t nTrksToSkip=0,id;
1664       AliExternalTrackParam *t = 0;
1665       for(Int_t i=0; i<nTrks; i++) {
1666         t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1667         id = (Int_t)t->GetID();
1668         if(id<0) continue;
1669         skipped[nTrksToSkip++] = id;
1670       }
1671       // TEMPORARY FIX
1672       // For AOD, skip also tracks without covariance matrix
1673       if(fInputAOD) {
1674         Double_t covtest[21];
1675         for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1676           AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1677           if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1678             id = (Int_t)vtrack->GetID();
1679             if(id<0) continue;
1680             skipped[nTrksToSkip++] = id;
1681           }
1682         }
1683       }
1684       //
1685       vertexer->SetSkipTracks(nTrksToSkip,skipped);
1686       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1687
1688     } else if(fRmTrksFromPrimVtx) {
1689       // removing the prongs tracks
1690
1691       TObjArray rmArray(nTrks);
1692       UShort_t *rmId = new UShort_t[nTrks];
1693       AliESDtrack *esdTrack = 0;
1694       AliESDtrack *t = 0;
1695       for(Int_t i=0; i<nTrks; i++) {
1696         t = (AliESDtrack*)trkArray->UncheckedAt(i);
1697         esdTrack = new AliESDtrack(*t);
1698         rmArray.AddLast(esdTrack);
1699         if(esdTrack->GetID()>=0) {
1700           rmId[i]=(UShort_t)esdTrack->GetID();
1701         } else {
1702           rmId[i]=9999;
1703         }
1704       }
1705       Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
1706       vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1707       delete [] rmId; rmId=NULL;
1708       rmArray.Delete();
1709
1710     }
1711
1712     if(!vertexESD) return vertexAOD;
1713     if(vertexESD->GetNContributors()<=0) {
1714       AliDebug(2,"vertexing failed");
1715       delete vertexESD; vertexESD=NULL;
1716       return vertexAOD;
1717     }
1718
1719     delete vertexer; vertexer=NULL;
1720
1721   }
1722
1723   // convert to AliAODVertex
1724   Double_t pos[3],cov[6],chi2perNDF;
1725   vertexESD->GetXYZ(pos); // position
1726   vertexESD->GetCovMatrix(cov); //covariance matrix
1727   chi2perNDF = vertexESD->GetChi2toNDF();
1728   delete vertexESD; vertexESD=NULL;
1729
1730   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1731
1732   return vertexAOD;
1733 }
1734 //-----------------------------------------------------------------------------
1735 void AliAnalysisVertexingHF::PrintStatus() const {
1736   // Print parameters being used
1737
1738   //printf("Preselections:\n");
1739   //   fTrackFilter->Dump();
1740   if(fSecVtxWithKF) {
1741     printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1742   } else {
1743     printf("Secondary vertex with AliVertexerTracks\n");
1744   }
1745   if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1746   if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1747   if(fD0toKpi) {
1748     printf("Reconstruct D0->Kpi candidates with cuts:\n");
1749     if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
1750   }
1751   if(fDstar) {
1752     printf("Reconstruct D*->D0pi candidates with cuts:\n");
1753     if(fFindVertexForDstar) {
1754       printf("    Reconstruct a secondary vertex for the D*\n");
1755     } else {
1756       printf("    Assume the D* comes from the primary vertex\n");
1757     }
1758     if(fCutsDStartoKpipi) fCutsDStartoKpipi->PrintAll();
1759   }
1760   if(fJPSItoEle) {
1761     printf("Reconstruct J/psi from B candidates with cuts:\n");
1762     if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
1763   }
1764   if(f3Prong) {
1765     printf("Reconstruct 3 prong candidates.\n");
1766     printf("  D+->Kpipi cuts:\n");
1767     if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
1768     printf("  Ds->KKpi cuts:\n");
1769     if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
1770     printf("  Lc->pKpi cuts:\n");
1771     if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
1772   }
1773   if(f4Prong) {
1774     printf("Reconstruct 4 prong candidates.\n");
1775     printf("  D0->Kpipipi cuts:\n");
1776     if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
1777   }
1778   if(fCascades) {
1779     printf("Reconstruct cascades candidates formed with v0s.\n");
1780     printf("  Lc -> k0s P & Lc -> L Pi cuts:\n");
1781     if(fCutsLctoV0) fCutsLctoV0->PrintAll();
1782   }
1783
1784   return;
1785 }
1786 //-----------------------------------------------------------------------------
1787 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
1788                                                                  Double_t &dispersion,Bool_t useTRefArray) const
1789 {
1790   // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1791
1792   AliESDVertex *vertexESD = 0;
1793   AliAODVertex *vertexAOD = 0;
1794
1795   if(!fSecVtxWithKF) { // AliVertexerTracks
1796
1797     AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1798     vertexer->SetVtxStart(fV1);
1799     vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1800     delete vertexer; vertexer=NULL;
1801
1802     if(!vertexESD) return vertexAOD;
1803
1804     if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
1805       AliDebug(2,"vertexing failed");
1806       delete vertexESD; vertexESD=NULL;
1807       return vertexAOD;
1808     }
1809
1810   } else { // Kalman Filter vertexer (AliKFParticle)
1811
1812     AliKFParticle::SetField(fBzkG);
1813
1814     AliKFVertex vertexKF;
1815
1816     Int_t nTrks = trkArray->GetEntriesFast();
1817     for(Int_t i=0; i<nTrks; i++) {
1818       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1819       AliKFParticle daughterKF(*esdTrack,211);
1820       vertexKF.AddDaughter(daughterKF);
1821     }
1822     vertexESD = new AliESDVertex(vertexKF.Parameters(),
1823                                  vertexKF.CovarianceMatrix(),
1824                                  vertexKF.GetChi2(),
1825                                  vertexKF.GetNContributors());
1826
1827   }
1828
1829   // convert to AliAODVertex
1830   Double_t pos[3],cov[6],chi2perNDF;
1831   vertexESD->GetXYZ(pos); // position
1832   vertexESD->GetCovMatrix(cov); //covariance matrix
1833   chi2perNDF = vertexESD->GetChi2toNDF();
1834   dispersion = vertexESD->GetDispersion();
1835   delete vertexESD; vertexESD=NULL;
1836
1837   Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
1838   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
1839
1840   return vertexAOD;
1841 }
1842 //-----------------------------------------------------------------------------
1843 Bool_t AliAnalysisVertexingHF::SelectInvMass(TObjArray *trkArray) const {
1844   // Invariant mass cut on tracks
1845
1846   Int_t nProngs=trkArray->GetEntriesFast();
1847   Int_t retval=kFALSE;
1848
1849   Double_t momentum[3];
1850   Double_t px3[3],py3[3],pz3[3];
1851   Double_t px4[4],py4[4],pz4[4];
1852   AliESDtrack *postrack1=0;
1853   AliESDtrack *postrack2=0;
1854   AliESDtrack *negtrack1=0;
1855   AliESDtrack *negtrack2=0;
1856
1857   switch(nProngs) {
1858   case 3:
1859     // invariant mass cut for D+, Ds, Lc
1860     postrack1 = (AliESDtrack*)trkArray->UncheckedAt(0);
1861     negtrack1  = (AliESDtrack*)trkArray->UncheckedAt(1);
1862     postrack2 = (AliESDtrack*)trkArray->UncheckedAt(2);
1863     postrack1->GetPxPyPz(momentum);
1864     px3[0] = momentum[0]; py3[0] = momentum[1]; pz3[0] = momentum[2];
1865     negtrack1->GetPxPyPz(momentum);
1866     px3[1] = momentum[0]; py3[1] = momentum[1]; pz3[1] = momentum[2];
1867     postrack2->GetPxPyPz(momentum);
1868     px3[2] = momentum[0]; py3[2] = momentum[1]; pz3[2] = momentum[2];
1869     retval = SelectInvMass(2,3,px3,py3,pz3);
1870     break;
1871   case 4:
1872     // invariant mass cut for D0->4prong
1873     postrack1 = (AliESDtrack*)trkArray->UncheckedAt(0);
1874     negtrack1 = (AliESDtrack*)trkArray->UncheckedAt(1);
1875     postrack2 = (AliESDtrack*)trkArray->UncheckedAt(2);
1876     negtrack2 = (AliESDtrack*)trkArray->UncheckedAt(3);
1877     postrack1->GetPxPyPz(momentum);
1878     px4[0] = momentum[0]; py4[0] = momentum[1]; pz4[0] = momentum[2];
1879     negtrack1->GetPxPyPz(momentum);
1880     px4[1] = momentum[0]; py4[1] = momentum[1]; pz4[1] = momentum[2];
1881     postrack2->GetPxPyPz(momentum);
1882     px4[2] = momentum[0]; py4[2] = momentum[1]; pz4[2] = momentum[2];
1883     negtrack2->GetPxPyPz(momentum);
1884     px4[3] = momentum[0]; py4[3] = momentum[1]; pz4[3] = momentum[2];
1885     retval = SelectInvMass(4,4,px4,py4,pz4);
1886     break;
1887   default:
1888     printf("SelectInvMass(): wrong decay selection\n");
1889     break;
1890   }
1891
1892   return retval;
1893 }
1894 //-----------------------------------------------------------------------------
1895 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1896                                              Int_t nprongs,
1897                                              Double_t *px,
1898                                              Double_t *py,
1899                                              Double_t *pz) const {
1900   // Check invariant mass cut
1901
1902   UInt_t pdg2[2],pdg3[3],pdg4[4];
1903   Double_t mPDG,minv;
1904
1905   Bool_t retval=kFALSE;
1906   switch (decay)
1907     {
1908     case 0:                  // D0->Kpi
1909       fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
1910       pdg2[0]=211; pdg2[1]=321;
1911       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1912       minv = fMassCalc2->InvMass(nprongs,pdg2);
1913       if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
1914       pdg2[0]=321; pdg2[1]=211;
1915       minv = fMassCalc2->InvMass(nprongs,pdg2);
1916       if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
1917       break;
1918     case 1:                  // JPSI->ee
1919       fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
1920       pdg2[0]=11; pdg2[1]=11;
1921       mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
1922       minv = fMassCalc2->InvMass(nprongs,pdg2);
1923       if(TMath::Abs(minv-mPDG)<fCutsJpsitoee->GetMassCut()) retval=kTRUE;
1924       break;
1925     case 2:                  // D+->Kpipi
1926       fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
1927       pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
1928       mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
1929       minv = fMassCalc3->InvMass(nprongs,pdg3);
1930       if(TMath::Abs(minv-mPDG)<fCutsDplustoKpipi->GetMassCut()) retval=kTRUE;
1931                             // Ds+->KKpi
1932       pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
1933       mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1934       minv = fMassCalc3->InvMass(nprongs,pdg3);
1935       if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
1936       pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
1937       minv = fMassCalc3->InvMass(nprongs,pdg3);
1938       if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
1939                             // Lc->pKpi
1940       pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
1941       mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1942       minv = fMassCalc3->InvMass(nprongs,pdg3);
1943       if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
1944       pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
1945       minv = fMassCalc3->InvMass(nprongs,pdg3);
1946       if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
1947       break;
1948     case 3:                  // D*->D0pi
1949       fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
1950       pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
1951       mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
1952       minv = fMassCalc2->InvMass(nprongs,pdg2);
1953       if(TMath::Abs(minv-mPDG)<fCutsDStartoKpipi->GetMassCut()) retval=kTRUE;
1954       break;
1955     case 4:                 // D0->Kpipipi without PID
1956       fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
1957       pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
1958       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1959       minv = fMassCalc4->InvMass(nprongs,pdg4);
1960       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1961       pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
1962       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1963       minv = fMassCalc4->InvMass(nprongs,pdg4);
1964       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1965       pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
1966       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1967       minv = fMassCalc4->InvMass(nprongs,pdg4);
1968       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1969       pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
1970       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1971       minv = fMassCalc4->InvMass(nprongs,pdg4);
1972       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1973       break;
1974     default:
1975       printf("SelectInvMass(): wrong decay selection\n");
1976       break;
1977     }
1978
1979   return retval;
1980 }
1981 //-----------------------------------------------------------------------------
1982 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
1983                                    TObjArray &seleTrksArray,Int_t &nSeleTrks,
1984                                    UChar_t *seleFlags,Int_t *evtNumber)
1985 {
1986   // Apply single-track preselection.
1987   // Fill a TObjArray with selected tracks (for displaced vertices or
1988   // soft pion from D*). Selection flag stored in seleFlags.
1989   // Create the AliESDVertex object (convert from AliAODVertex if necessary)
1990   // In case of AOD input, also fill fAODMap for track index<->ID
1991
1992   const AliVVertex *vprimary = event->GetPrimaryVertex();
1993
1994   if(fV1) { delete fV1; fV1=NULL; }
1995   if(fAODMap) { delete fAODMap; fAODMap=NULL; }
1996
1997   Int_t nindices=0;
1998   UShort_t *indices = 0;
1999   Double_t pos[3],cov[6];
2000
2001   if(!fInputAOD) { // ESD
2002     fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
2003   } else {         // AOD
2004     vprimary->GetXYZ(pos);
2005     vprimary->GetCovarianceMatrix(cov);
2006     fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2007     indices = new UShort_t[event->GetNumberOfTracks()];
2008     fAODMapSize = 100000;
2009     fAODMap = new Int_t[fAODMapSize];
2010   }
2011
2012
2013   Int_t entries = (Int_t)event->GetNumberOfTracks();
2014   Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
2015   nSeleTrks=0;
2016
2017   // transfer ITS tracks from event to arrays
2018   for(Int_t i=0; i<entries; i++) {
2019     AliVTrack *track;
2020     track = (AliVTrack*)event->GetTrack(i);
2021
2022     // skip pure ITS SA tracks
2023     if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2024
2025     // TEMPORARY: check that the cov matrix is there
2026     Double_t covtest[21];
2027     if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2028     //
2029
2030     if(fInputAOD) {
2031       AliAODTrack *aodt = (AliAODTrack*)track;
2032       if(aodt->GetUsedForPrimVtxFit()) {
2033         indices[nindices]=aodt->GetID(); nindices++;
2034       }
2035       fAODMap[(Int_t)aodt->GetID()] = i;
2036     }
2037
2038     AliESDtrack *esdt = 0;
2039
2040     if(!fInputAOD) {
2041       esdt = (AliESDtrack*)track;
2042     } else {
2043       esdt = new AliESDtrack(track);
2044     }
2045
2046     // single track selection
2047     okDisplaced=kFALSE; okSoftPi=kFALSE;
2048     if(fMixEvent){
2049       evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
2050       const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
2051       Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
2052       eventVtx->GetXYZ(vtxPos);
2053       vprimary->GetXYZ(primPos);
2054       eventVtx->GetCovarianceMatrix(primCov);
2055       for(Int_t ind=0;ind<3;ind++){
2056         trasl[ind]=vtxPos[ind]-primPos[ind];
2057       }
2058
2059       Bool_t isTransl=esdt->Translate(trasl,primCov);
2060       if(!isTransl) {
2061         delete esdt;
2062         esdt = NULL;
2063         continue;
2064       }
2065     }
2066
2067     if(SingleTrkCuts(esdt,okDisplaced,okSoftPi)) {
2068       seleTrksArray.AddLast(esdt);
2069       seleFlags[nSeleTrks]=0;
2070       if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2071       if(okSoftPi)    SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2072       nSeleTrks++;
2073     } else {
2074       if(fInputAOD) delete esdt;
2075       esdt = NULL;
2076       continue;
2077     }
2078
2079   } // end loop on tracks
2080
2081   // primary vertex from AOD
2082   if(fInputAOD) {
2083     delete fV1; fV1=NULL;
2084     vprimary->GetXYZ(pos);
2085     vprimary->GetCovarianceMatrix(cov);
2086     Double_t chi2toNDF = vprimary->GetChi2perNDF();
2087     Int_t ncontr=nindices;
2088     if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
2089     Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
2090     fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2091     fV1->SetTitle(vprimary->GetTitle());
2092     fV1->SetIndices(nindices,indices);
2093   }
2094   if(indices) { delete [] indices; indices=NULL; }
2095
2096
2097   return;
2098 }
2099 //-----------------------------------------------------------------------------
2100 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2101                                              Bool_t &okDisplaced,Bool_t &okSoftPi) const
2102 {
2103   // Check if track passes some kinematical cuts
2104
2105   // this is needed to store the impact parameters
2106   trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2107
2108   UInt_t selectInfo;
2109   //
2110   // Track selection, displaced tracks
2111   selectInfo = 0;
2112   if(fTrackFilter) {
2113     selectInfo = fTrackFilter->IsSelected(trk);
2114   }
2115   if(selectInfo) okDisplaced=kTRUE;
2116   // Track selection, soft pions
2117   selectInfo = 0;
2118   if(fDstar && fTrackFilterSoftPi) {
2119     selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2120   }
2121   if(selectInfo) okSoftPi=kTRUE;
2122
2123   if(okDisplaced || okSoftPi) return kTRUE;
2124
2125   return kFALSE;
2126 }
2127
2128
2129 //-----------------------------------------------------------------------------
2130 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2131   //
2132   // Transform ESDv0 to AODv0
2133   //
2134   //  this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2135   //  and creates an AODv0 out of them
2136   //
2137   Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
2138   AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2139
2140   // create the v0 neutral track to compute the DCA to the primary vertex
2141   Double_t xyz[3], pxpypz[3];
2142   esdV0->XvYvZv(xyz);
2143   esdV0->PxPyPz(pxpypz);
2144   Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2145   AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2146   if(!trackesdV0) return 0;
2147   Double_t d0z0[2],covd0z0[3];
2148   AliAODVertex *primVertexAOD = PrimaryVertex();
2149   trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2150   Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2151   // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2152   Double_t dcaV0DaughterToPrimVertex[2];
2153   AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2154   AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
2155   if( !posV0track || !negV0track) {
2156     if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2157     return 0;
2158   }
2159   posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2160   //  if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2161   //  else
2162   dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
2163   negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2164   //  if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2165   //  else
2166   dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
2167   Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
2168   Double_t pmom[3],nmom[3];
2169   esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2170   esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2171
2172   AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2173   aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2174
2175   if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2176   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2177
2178   return aodV0;
2179 }
2180 //-----------------------------------------------------------------------------