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