AliPhysics  c7b8e89 (c7b8e89)
AliAnalysisVertexingHF.cxx
Go to the documentation of this file.
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 /* $Id$ */
17 
18 //----------------------------------------------------------------------------
19 // Implementation of the heavy-flavour vertexing analysis class
20 // Candidates are stored in the AOD as objects deriving from AliAODRecoDecay.
21 // To be used as a task of AliAnalysisManager by means of the interface
22 // class AliAnalysisTaskSEVertexingHF.
23 // An example of usage in the macro AliAnalysisTaskSEVertexingHFTest.C.
24 //
25 // Contact: andrea.dainese@pd.infn.it
26 // Contributors: E.Bruna, G.E.Bruno, A.Dainese, C.Di Gliglio,
27 // F.Prino, R.Romita, X.M.Zhang
28 //----------------------------------------------------------------------------
29 #include <Riostream.h>
30 #include <TFile.h>
31 #include <TDatabasePDG.h>
32 #include <TString.h>
33 #include <TList.h>
34 #include <TProcessID.h>
35 #include "AliLog.h"
36 #include "AliVEvent.h"
37 #include "AliVVertex.h"
38 #include "AliVTrack.h"
39 #include "AliVertexerTracks.h"
40 #include "AliKFVertex.h"
41 #include "AliESDEvent.h"
42 #include "AliESDVertex.h"
43 #include "AliExternalTrackParam.h"
44 #include "AliNeutralTrackParam.h"
45 #include "AliESDtrack.h"
46 #include "AliESDtrackCuts.h"
47 #include "AliAODEvent.h"
48 #include "AliPIDResponse.h"
49 #include "AliAODRecoDecay.h"
50 #include "AliAODRecoDecayHF.h"
54 #include "AliAODRecoCascadeHF.h"
55 #include "AliRDHFCutsD0toKpi.h"
56 #include "AliRDHFCutsJpsitoee.h"
59 #include "AliRDHFCutsDstoK0sK.h"
60 #include "AliRDHFCutsDstoKKpi.h"
61 #include "AliRDHFCutsLctopKpi.h"
62 #include "AliRDHFCutsLctoV0.h"
63 #include "AliRDHFCutsD0toKpipipi.h"
65 #include "AliAnalysisFilter.h"
66 #include "AliAnalysisVertexingHF.h"
67 #include "AliMixedEvent.h"
68 #include "AliESDv0.h"
69 #include "AliAODv0.h"
70 #include "AliCodeTimer.h"
71 #include "AliMultSelection.h"
72 #include <cstring>
73 
75 ClassImp(AliAnalysisVertexingHF);
77 
78 //----------------------------------------------------------------------------
80 fInputAOD(kFALSE),
81 fAODMapSize(0),
82 fAODMap(0),
83 fVertexerTracks(0x0),
84 fBzkG(0.),
85 fSecVtxWithKF(kFALSE),
86 fRecoPrimVtxSkippingTrks(kFALSE),
87 fRmTrksFromPrimVtx(kFALSE),
88 fV1(0x0),
89 fD0toKpi(kTRUE),
90 fJPSItoEle(kTRUE),
91 f3Prong(kTRUE),
92 f4Prong(kTRUE),
93 fDstar(kTRUE),
94 fCascades(kTRUE),
95 fLikeSign(kFALSE),
96 fLikeSign3prong(kFALSE),
97 fMixEvent(kFALSE),
98 fPidResponse(0x0),
99 fUseKaonPIDfor3Prong(kFALSE),
100 fUsePIDforLc(0),
101 fUsePIDforLc2V0(0),
102 fUseKaonPIDforDs(kFALSE),
103 fUseTPCPID(kFALSE),
104 fUseTOFPID(kFALSE),
105 fUseTPCPIDOnlyIfNoTOF(kFALSE),
106 fMaxMomForTPCPid(1.),
107 fUsePidTag(kFALSE),
108 fnSigmaTPCPionLow(5.),
109 fnSigmaTPCPionHi(5.),
110 fnSigmaTOFPionLow(5.),
111 fnSigmaTOFPionHi(5.),
112 fnSigmaTPCKaonLow(5.),
113 fnSigmaTPCKaonHi(5.),
114 fnSigmaTOFKaonLow(5.),
115 fnSigmaTOFKaonHi(5.),
116 fnSigmaTPCProtonLow(5.),
117 fnSigmaTPCProtonHi(5.),
118 fnSigmaTOFProtonLow(5.),
119 fnSigmaTOFProtonHi(5.),
120 fMaxCentPercentileForTightCuts(-9999),
121 fTrackFilter(0x0),
122 fTrackFilter2prongCentral(0x0),
123 fTrackFilter3prongCentral(0x0),
124 fTrackFilterSoftPi(0x0),
125 fTrackFilterBachelor(0x0),
126 fCutsD0toKpi(0x0),
127 fCutsJpsitoee(0x0),
128 fCutsDplustoK0spi(0x0),
129 fCutsDplustoKpipi(0x0),
130 fCutsDstoK0sK(0x0),
131 fCutsDstoKKpi(0x0),
132 fCutsLctopKpi(0x0),
133 fCutsLctoV0(0x0),
134 fCutsD0toKpipipi(0x0),
135 fCutsDStartoKpipi(0x0),
136 fListOfCuts(0x0),
137 fFindVertexForDstar(kTRUE),
138 fFindVertexForCascades(kTRUE),
139 fV0TypeForCascadeVertex(0),
140 fMassCutBeforeVertexing(kFALSE),
141 fMassCalc2(0),
142 fMassCalc3(0),
143 fMassCalc4(0),
144 fOKInvMassD0(kFALSE),
145 fOKInvMassJpsi(kFALSE),
146 fOKInvMassDplus(kFALSE),
147 fOKInvMassDs(kFALSE),
148 fOKInvMassLc(kFALSE),
149 fOKInvMassDstar(kFALSE),
150 fOKInvMassD0to4p(kFALSE),
151 fOKInvMassLctoV0(kFALSE),
152 fnTrksTotal(0),
153 fnSeleTrksTotal(0),
154 fMakeReducedRHF(kFALSE),
155 fMassDzero(0.),
156 fMassDplus(0.),
157 fMassDs(0.),
158 fMassLambdaC(0.),
159 fMassDstar(0.),
160 fMassJpsi(0.)
161 {
163 
164  Double_t d02[2]={0.,0.};
165  Double_t d03[3]={0.,0.,0.};
166  Double_t d04[4]={0.,0.,0.,0.};
167  fMassCalc2 = new AliAODRecoDecay(0x0,2,0,d02);
168  fMassCalc3 = new AliAODRecoDecay(0x0,3,1,d03);
169  fMassCalc4 = new AliAODRecoDecay(0x0,4,0,d04);
170  SetMasses();
171 }
172 //--------------------------------------------------------------------------
174 TNamed(source),
175 fInputAOD(source.fInputAOD),
176 fAODMapSize(source.fAODMapSize),
177 fAODMap(source.fAODMap),
179 fBzkG(source.fBzkG),
183 fV1(source.fV1),
184 fD0toKpi(source.fD0toKpi),
185 fJPSItoEle(source.fJPSItoEle),
186 f3Prong(source.f3Prong),
187 f4Prong(source.f4Prong),
188 fDstar(source.fDstar),
189 fCascades(source.fCascades),
190 fLikeSign(source.fLikeSign),
192 fMixEvent(source.fMixEvent),
193 fPidResponse(source.fPidResponse),
195 fUsePIDforLc(source.fUsePIDforLc),
198 fUseTPCPID(source.fUseTPCPID),
199 fUseTOFPID(source.fUseTOFPID),
202 fUsePidTag(source.fUsePidTag),
216 fTrackFilter(source.fTrackFilter),
221 fCutsD0toKpi(source.fCutsD0toKpi),
228 fCutsLctoV0(source.fCutsLctoV0),
231 fListOfCuts(source.fListOfCuts),
236 fMassCalc2(source.fMassCalc2),
237 fMassCalc3(source.fMassCalc3),
238 fMassCalc4(source.fMassCalc4),
239 fOKInvMassD0(source.fOKInvMassD0),
242 fOKInvMassDs(source.fOKInvMassDs),
243 fOKInvMassLc(source.fOKInvMassLc),
247 fnTrksTotal(0),
248 fnSeleTrksTotal(0),
249 fMakeReducedRHF(kFALSE),
250 fMassDzero(source.fMassDzero),
251 fMassDplus(source.fMassDplus),
252 fMassDs(source.fMassDs),
253 fMassLambdaC(source.fMassLambdaC),
254 fMassDstar(source.fMassDstar),
255 fMassJpsi(source.fMassJpsi)
256 {
260 }
261 //--------------------------------------------------------------------------
263 {
267  if(&source == this) return *this;
268  fInputAOD = source.fInputAOD;
269  fAODMapSize = source.fAODMapSize;
271  fBzkG = source.fBzkG;
272  fSecVtxWithKF = source.fSecVtxWithKF;
275  fV1 = source.fV1;
276  fD0toKpi = source.fD0toKpi;
277  fJPSItoEle = source.fJPSItoEle;
278  f3Prong = source.f3Prong;
279  f4Prong = source.f4Prong;
280  fDstar = source.fDstar;
281  fCascades = source.fCascades;
282  fLikeSign = source.fLikeSign;
284  fMixEvent = source.fMixEvent;
285  fPidResponse = source.fPidResponse;
287  fUsePIDforLc = source.fUsePIDforLc;
290  fUseTPCPID = source.fUseTPCPID;
291  fUseTOFPID = source.fUseTOFPID;
294  fUsePidTag = source.fUsePidTag;
310  fTrackFilter = source.fTrackFilter;
315  fCutsD0toKpi = source.fCutsD0toKpi;
316  fCutsJpsitoee = source.fCutsJpsitoee;
319  fCutsDstoK0sK = source.fCutsDstoK0sK;
320  fCutsDstoKKpi = source.fCutsDstoKKpi;
321  fCutsLctopKpi = source.fCutsLctopKpi;
322  fCutsLctoV0 = source.fCutsLctoV0;
325  fListOfCuts = source.fListOfCuts;
330  fMassCalc2 = source.fMassCalc2;
331  fMassCalc3 = source.fMassCalc3;
332  fMassCalc4 = source.fMassCalc4;
333  fOKInvMassD0 = source.fOKInvMassD0;
336  fOKInvMassDs = source.fOKInvMassDs;
337  fOKInvMassLc = source.fOKInvMassLc;
341  fMassDzero = source.fMassDzero;
342  fMassDplus = source.fMassDplus;
343  fMassDs = source.fMassDs;
344  fMassLambdaC = source.fMassLambdaC;
345  fMassDstar = source.fMassDstar;
346  fMassJpsi = source.fMassJpsi;
347 
348  return *this;
349 }
350 //----------------------------------------------------------------------------
353  if(fV1) { delete fV1; fV1=0; }
354  delete fVertexerTracks;
355  if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
360  if(fCutsD0toKpi) { delete fCutsD0toKpi; fCutsD0toKpi=0; }
361  if(fCutsJpsitoee) { delete fCutsJpsitoee; fCutsJpsitoee=0; }
364  if(fCutsDstoK0sK) { delete fCutsDstoK0sK; fCutsDstoK0sK=0; }
365  if(fCutsDstoKKpi) { delete fCutsDstoKKpi; fCutsDstoKKpi=0; }
366  if(fCutsLctopKpi) { delete fCutsLctopKpi; fCutsLctopKpi=0; }
367  if(fCutsLctoV0) { delete fCutsLctoV0; fCutsLctoV0=0; }
370  if(fAODMap) { delete [] fAODMap; fAODMap=0; }
371  if(fMassCalc2) { delete fMassCalc2; fMassCalc2=0; }
372  if(fMassCalc3) { delete fMassCalc3; fMassCalc3=0; }
373  if(fMassCalc4) { delete fMassCalc4; fMassCalc4=0; }
374 }
375 //----------------------------------------------------------------------------
378 
379  TList *list = new TList();
380  list->SetOwner();
381  list->SetName("ListOfCuts");
382 
383  if(fCutsD0toKpi) {
385  list->Add(cutsD0toKpi);
386  }
387  if(fCutsJpsitoee) {
389  list->Add(cutsJpsitoee);
390  }
391  if (fCutsDplustoK0spi) {
393  list->Add(cutsDplustoK0spi);
394  }
395  if(fCutsDplustoKpipi) {
397  list->Add(cutsDplustoKpipi);
398  }
399  if (fCutsDstoK0sK) {
401  list->Add(cutsDstoK0sK);
402  }
403  if(fCutsDstoKKpi) {
405  list->Add(cutsDstoKKpi);
406  }
407  if(fCutsLctopKpi) {
409  list->Add(cutsLctopKpi);
410  }
411  if(fCutsLctoV0){
412  AliRDHFCutsLctoV0 *cutsLctoV0 = new AliRDHFCutsLctoV0(*fCutsLctoV0);
413  list->Add(cutsLctoV0);
414  }
415  if(fCutsD0toKpipipi) {
417  list->Add(cutsD0toKpipipi);
418  }
419  if(fCutsDStartoKpipi) {
421  list->Add(cutsDStartoKpipi);
422  }
423 
424  //___ Check consitstency of cuts between vertexer and analysis tasks
425  Bool_t bCutsOk = CheckCutsConsistency();
426  if (bCutsOk == kFALSE) {AliFatal("AliAnalysisVertexingHF::FillListOfCuts vertexing and the analysis task cuts are not consistent!");}
427 
428  // keep a pointer to the list
429  fListOfCuts = list;
430 
431  return list;
432 }
433 //----------------------------------------------------------------------------
435  TClonesArray *aodVerticesHFTClArr,
436  TClonesArray *aodD0toKpiTClArr,
437  TClonesArray *aodJPSItoEleTClArr,
438  TClonesArray *aodCharm3ProngTClArr,
439  TClonesArray *aodCharm4ProngTClArr,
440  TClonesArray *aodDstarTClArr,
441  TClonesArray *aodCascadesTClArr,
442  TClonesArray *aodLikeSign2ProngTClArr,
443  TClonesArray *aodLikeSign3ProngTClArr)
444 {
448  //AliCodeTimerAuto("",0);
449 
450  if(!fMixEvent){
451  TString evtype = event->IsA()->GetName();
452  fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
453  } // if we do mixing AliVEvent is a AliMixedEvent
454 
455  if(fInputAOD) {
456  AliDebug(2,"Creating HF candidates from AOD");
457  } else {
458  AliDebug(2,"Creating HF candidates from ESD");
459  }
460 
461  if(!aodVerticesHFTClArr) {
462  printf("ERROR: no aodVerticesHFTClArr");
463  return;
464  }
465  if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
466  printf("ERROR: no aodD0toKpiTClArr");
467  return;
468  }
469  if(fJPSItoEle && !aodJPSItoEleTClArr) {
470  printf("ERROR: no aodJPSItoEleTClArr");
471  return;
472  }
473  if(f3Prong && !aodCharm3ProngTClArr) {
474  printf("ERROR: no aodCharm3ProngTClArr");
475  return;
476  }
477  if(f4Prong && !aodCharm4ProngTClArr) {
478  printf("ERROR: no aodCharm4ProngTClArr");
479  return;
480  }
481  if(fDstar && !aodDstarTClArr) {
482  printf("ERROR: no aodDstarTClArr");
483  return;
484  }
485  if(fCascades && !aodCascadesTClArr){
486  printf("ERROR: no aodCascadesTClArr ");
487  return;
488  }
489  if(fLikeSign && !aodLikeSign2ProngTClArr) {
490  printf("ERROR: no aodLikeSign2ProngTClArr");
491  return;
492  }
493  if(fLikeSign3prong && f3Prong && !aodLikeSign3ProngTClArr) {
494  printf("ERROR: no aodLikeSign3ProngTClArr");
495  return;
496  }
497 
498  // delete candidates from previous event and create references
499  Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0,iCascades=0,iLikeSign2Prong=0,iLikeSign3Prong=0;
500  aodVerticesHFTClArr->Delete();
501  iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
502  TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
503  if(fD0toKpi || fDstar) {
504  aodD0toKpiTClArr->Delete();
505  iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
506  }
507  if(fJPSItoEle) {
508  aodJPSItoEleTClArr->Delete();
509  iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
510  }
511  if(f3Prong) {
512  aodCharm3ProngTClArr->Delete();
513  i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
514  }
515  if(f4Prong) {
516  aodCharm4ProngTClArr->Delete();
517  i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
518  }
519  if(fDstar) {
520  aodDstarTClArr->Delete();
521  iDstar = aodDstarTClArr->GetEntriesFast();
522  }
523  if(fCascades) {
524  aodCascadesTClArr->Delete();
525  iCascades = aodCascadesTClArr->GetEntriesFast();
526  }
527  if(fLikeSign) {
528  aodLikeSign2ProngTClArr->Delete();
529  iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast();
530  }
531  if(fLikeSign3prong && f3Prong) {
532  aodLikeSign3ProngTClArr->Delete();
533  iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast();
534  }
535 
536  TClonesArray &aodD0toKpiRef = *aodD0toKpiTClArr;
537  TClonesArray &aodJPSItoEleRef = *aodJPSItoEleTClArr;
538  TClonesArray &aodCharm3ProngRef = *aodCharm3ProngTClArr;
539  TClonesArray &aodCharm4ProngRef = *aodCharm4ProngTClArr;
540  TClonesArray &aodDstarRef = *aodDstarTClArr;
541  TClonesArray &aodCascadesRef = *aodCascadesTClArr;
542  TClonesArray &aodLikeSign2ProngRef = *aodLikeSign2ProngTClArr;
543  TClonesArray &aodLikeSign3ProngRef = *aodLikeSign3ProngTClArr;
544 
545 
546  AliAODRecoDecayHF2Prong *io2Prong = 0;
547  AliAODRecoDecayHF3Prong *io3Prong = 0;
548  AliAODRecoDecayHF4Prong *io4Prong = 0;
549  AliAODRecoCascadeHF *ioCascade = 0;
550 
551  Int_t iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries,iv0,nv0;
552  Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaV0,dcaCasc;
553  Bool_t okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
554  Bool_t okDstar=kFALSE,okD0fromDstar=kFALSE;
555  Bool_t okCascades=kFALSE;
556  AliESDtrack *postrack1 = 0;
557  AliESDtrack *postrack2 = 0;
558  AliESDtrack *negtrack1 = 0;
559  AliESDtrack *negtrack2 = 0;
560  AliESDtrack *trackPi = 0;
561  Double_t mompos1[3],mompos2[3],momneg1[3],momneg2[3];
562  // AliESDtrack *posV0track = 0;
563  // AliESDtrack *negV0track = 0;
564  Float_t dcaMax = fCutsD0toKpi->GetDCACut();
565  if(fCutsJpsitoee) dcaMax=TMath::Max(dcaMax,fCutsJpsitoee->GetDCACut());
566  if(fCutsDplustoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDplustoKpipi->GetDCACut());
567  if(fCutsDstoKKpi) dcaMax=TMath::Max(dcaMax,fCutsDstoKKpi->GetDCACut());
568  if(fCutsLctopKpi) dcaMax=TMath::Max(dcaMax,fCutsLctopKpi->GetDCACut());
569  if(fCutsD0toKpipipi) dcaMax=TMath::Max(dcaMax,fCutsD0toKpipipi->GetDCACut());
570  if(fCutsDStartoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDStartoKpipi->GetDCACut());
571 
572  AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
573 
574 
575  // get Bz
576  fBzkG = (Double_t)event->GetMagneticField();
577  if(!fVertexerTracks){
578  fVertexerTracks=new AliVertexerTracks(fBzkG);
579  }else{
580  Double_t oldField=fVertexerTracks->GetFieldkG();
581  if(oldField!=fBzkG) fVertexerTracks->SetFieldkG(fBzkG);
582  }
583 
584  trkEntries = (Int_t)event->GetNumberOfTracks();
585  AliDebug(1,Form(" Number of tracks: %d",trkEntries));
586  fnTrksTotal += trkEntries;
587 
588  nv0 = (Int_t)event->GetNumberOfV0s();
589  AliDebug(1,Form(" Number of V0s: %d",nv0));
590 
591  if( trkEntries<2 && (trkEntries<1 || nv0<1) ) {
592  AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
593  return;
594  }
595 
596  // event selection + PID configuration
597  if(!fCutsD0toKpi->IsEventSelected(event)) return;
604  if(fCutsLctoV0) fCutsLctoV0->SetupPID(event);
607 
608  // call function that applies sigle-track selection,
609  // for displaced tracks and soft pions (both charges) for D*,
610  // and retrieves primary vertex
611  TObjArray seleTrksArray(trkEntries);
612  TObjArray tracksAtVertex(trkEntries);
613  UChar_t *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi, bit 2: 3 prong, bits 3-4-5: for PID
614  Int_t nSeleTrks=0;
615  Int_t *evtNumber = new Int_t[trkEntries];
616  SelectTracksAndCopyVertex(event,trkEntries,seleTrksArray,tracksAtVertex,nSeleTrks,seleFlags,evtNumber);
617 
618  AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
619  fnSeleTrksTotal += nSeleTrks;
620 
621 
622  TObjArray *twoTrackArray1 = new TObjArray(2);
623  TObjArray *twoTrackArray2 = new TObjArray(2);
624  TObjArray *twoTrackArrayV0 = new TObjArray(2);
625  TObjArray *twoTrackArrayCasc = new TObjArray(2);
626  TObjArray *threeTrackArray = new TObjArray(3);
627  TObjArray *fourTrackArray = new TObjArray(4);
628 
629  Double_t dispersion;
630  Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
631 
632  AliAODRecoDecayHF *rd = 0;
633  AliAODRecoCascadeHF *rc = 0;
634  AliAODv0 *v0 = 0;
635  AliESDv0 *esdV0 = 0;
636 
637  Bool_t massCutOK=kTRUE;
638 
639  // LOOP ON POSITIVE TRACKS
640  for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
641 
642  //if(iTrkP1%1==0) AliDebug(1,Form(" 1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));
643  //if(iTrkP1%1==0) printf(" 1st loop on pos: track number %d of %d\n",iTrkP1,nSeleTrks);
644 
645  // get track from tracks array
646  postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
647  postrack1->GetPxPyPz(mompos1);
648 
649  // Make cascades with V0+track
650  //
651  if(fCascades) {
652  // loop on V0's
653  for(iv0=0; iv0<nv0; iv0++){
654 
655  //AliDebug(1,Form(" loop on v0s for track number %d and v0 number %d",iTrkP1,iv0));
656  if ( !TESTBIT(seleFlags[iTrkP1],kBitBachelor) ) continue;
657 
658  if ( fUsePIDforLc2V0 && !TESTBIT(seleFlags[iTrkP1],kBitProtonCompat) ) continue; //clm
659 
660  // Get the V0
661  if(fInputAOD) {
662  v0 = ((AliAODEvent*)event)->GetV0(iv0);
663  } else {
664  esdV0 = ((AliESDEvent*)event)->GetV0(iv0);
665  }
666  if ( (!v0 || !v0->IsA()->InheritsFrom("AliAODv0") ) &&
667  (!esdV0 || !esdV0->IsA()->InheritsFrom("AliESDv0") ) ) continue;
668 
669  if ( v0 && ((v0->GetOnFlyStatus() == kTRUE && fV0TypeForCascadeVertex == AliRDHFCuts::kOnlyOfflineV0s) ||
670  (v0->GetOnFlyStatus() == kFALSE && fV0TypeForCascadeVertex == AliRDHFCuts::kOnlyOnTheFlyV0s)) ) continue;
671 
672  if ( esdV0 && ((esdV0->GetOnFlyStatus() == kTRUE && fV0TypeForCascadeVertex == AliRDHFCuts::kOnlyOfflineV0s) ||
673  ( esdV0->GetOnFlyStatus() == kFALSE && fV0TypeForCascadeVertex == AliRDHFCuts::kOnlyOnTheFlyV0s)) ) continue;
674 
675  // Get the tracks that form the V0
676  // ( parameters at primary vertex )
677  // and define an AliExternalTrackParam out of them
678  AliExternalTrackParam * posV0track;
679  AliExternalTrackParam * negV0track;
680 
681  if(fInputAOD){
682  AliAODTrack *posVV0track = (AliAODTrack*)(v0->GetDaughter(0));
683  AliAODTrack *negVV0track = (AliAODTrack*)(v0->GetDaughter(1));
684  if( !posVV0track || !negVV0track ) continue;
685  //
686  // Apply some basic V0 daughter criteria
687  //
688  // bachelor must not be a v0-track
689  if (posVV0track->GetID() == postrack1->GetID() ||
690  negVV0track->GetID() == postrack1->GetID()) continue;
691  // reject like-sign v0
692  if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
693  // avoid ghost TPC tracks
694  if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
695  !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
696  // Get AliExternalTrackParam out of the AliAODTracks
697  Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign;
698  posVV0track->PxPyPz(pxpypz); posVV0track->XvYvZv(xyz);
699  posVV0track->GetCovarianceXYZPxPyPz(cv); sign=posVV0track->Charge();
700  posV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
701  negVV0track->PxPyPz(pxpypz); negVV0track->XvYvZv(xyz);
702  negVV0track->GetCovarianceXYZPxPyPz(cv); sign=negVV0track->Charge();
703  negV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
704  } else {
705  AliESDtrack *posVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetPindex() ));
706  AliESDtrack *negVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetNindex() ));
707  if( !posVV0track || !negVV0track ) continue;
708  //
709  // Apply some basic V0 daughter criteria
710  //
711  // bachelor must not be a v0-track
712  if (posVV0track->GetID() == postrack1->GetID() ||
713  negVV0track->GetID() == postrack1->GetID()) continue;
714  // reject like-sign v0
715  if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
716  // avoid ghost TPC tracks
717  if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
718  !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
719  // reject kinks (only necessary on AliESDtracks)
720  if (posVV0track->GetKinkIndex(0)>0 || negVV0track->GetKinkIndex(0)>0) continue;
721  // Get AliExternalTrackParam out of the AliESDtracks
722  posV0track = new AliExternalTrackParam(*posVV0track);
723  negV0track = new AliExternalTrackParam(*negVV0track);
724 
725  // Define the AODv0 from ESDv0 if reading ESDs
726  v0 = TransformESDv0toAODv0(esdV0,twoTrackArrayV0);
727  }
728  if( !posV0track || !negV0track ){
729  AliDebug(1,Form(" Couldn't get the V0 daughters"));
730  continue;
731  }
732 
733  // fill in the v0 two-external-track-param array
734  twoTrackArrayV0->AddAt(posV0track,0);
735  twoTrackArrayV0->AddAt(negV0track,1);
736 
737  // Get the V0 dca
738  dcaV0 = v0->DcaV0Daughters();
739 
740  // Define the V0 (neutral) track
741  AliNeutralTrackParam *trackV0=NULL;
742  if(fInputAOD) {
743  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
744  if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
745  } else {
746  Double_t xyz[3], pxpypz[3];
747  esdV0->XvYvZv(xyz);
748  esdV0->PxPyPz(pxpypz);
749  Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
750  trackV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
751  }
752 
753 
754  // Fill in the object array to create the cascade
755  twoTrackArrayCasc->AddAt(postrack1,0);
756  twoTrackArrayCasc->AddAt(trackV0,1);
757  // Compute the cascade vertex
758  AliAODVertex *vertexCasc = 0;
760  // DCA between the two tracks
761  dcaCasc = postrack1->GetDCA(trackV0,fBzkG,xdummy,ydummy);
762  // Vertexing+
763  vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
764  } else {
765  // assume Cascade decays at the primary vertex
766  Double_t pos[3],cov[6],chi2perNDF;
767  fV1->GetXYZ(pos);
768  fV1->GetCovMatrix(cov);
769  chi2perNDF = fV1->GetChi2toNDF();
770  vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
771  dcaCasc = 0.;
772  }
773  if(!vertexCasc) {
774  delete posV0track; posV0track=NULL;
775  delete negV0track; negV0track=NULL;
776  delete trackV0; trackV0=NULL;
777  if(!fInputAOD) {delete v0; v0=NULL;}
778  twoTrackArrayV0->Clear();
779  twoTrackArrayCasc->Clear();
780  continue;
781  }
782 
783  // Create and store the Cascade if passed the cuts
784  ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,v0,dcaCasc,okCascades);
785  if(okCascades && ioCascade) {
786  //AliDebug(1,Form("Storing a cascade object... "));
787  // add the vertex and the cascade to the AOD
788  rc = new(aodCascadesRef[iCascades++])AliAODRecoCascadeHF(*ioCascade);
789  if(fMakeReducedRHF){
790  UShort_t id[2]={(UShort_t)postrack1->GetID(),(UShort_t)iv0};
791  rc->SetProngIDs(2,id);
792  rc->DeleteRecoD();
793  }else{
794  AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
795  rc->SetSecondaryVtx(vCasc);
796  vCasc->SetParent(rc);
797  if(!fInputAOD) vCasc->AddDaughter(v0); // just to fill ref #0 ??
798  AddRefs(vCasc,rc,event,twoTrackArrayCasc); // add the track (proton)
799  vCasc->AddDaughter(v0); // fill the 2prong V0
800  }
801  rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
802  }
803 
804 
805  // Clean up
806  delete posV0track; posV0track=NULL;
807  delete negV0track; negV0track=NULL;
808  delete trackV0; trackV0=NULL;
809  twoTrackArrayV0->Clear();
810  twoTrackArrayCasc->Clear();
811  if(ioCascade) { delete ioCascade; ioCascade=NULL; }
812  if(vertexCasc) { delete vertexCasc; vertexCasc=NULL; }
813  if(!fInputAOD) {delete v0; v0=NULL;}
814 
815  } // end loop on V0's
816 
817  // re-set parameters at vertex
818  SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
819  } // end fCascades
820 
821  // If there is less than 2 particles continue
822  if(trkEntries<2) {
823  AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
824  continue;
825  }
826 
827  if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
828  if(postrack1->Charge()<0 && !fLikeSign) continue;
829 
830  // LOOP ON NEGATIVE TRACKS
831  for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
832 
833  //if(iTrkN1%1==0) AliDebug(1,Form(" 1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));
834  //if(iTrkN1%1==0) printf(" 1st loop on neg: track number %d of %d\n",iTrkN1,nSeleTrks);
835 
836  if(iTrkN1==iTrkP1) continue;
837 
838  // get track from tracks array
839  negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
840 
841  if(negtrack1->Charge()>0 && !fLikeSign) continue;
842 
843  if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
844 
845  if(fMixEvent) {
846  if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
847  }
848 
849  if(postrack1->Charge()==negtrack1->Charge()) { // like-sign
850  isLikeSign2Prong=kTRUE;
851  if(!fLikeSign) continue;
852  if(iTrkN1<iTrkP1) continue; // this is needed to avoid double-counting of like-sign
853  } else { // unlike-sign
854  isLikeSign2Prong=kFALSE;
855  if(postrack1->Charge()<0 || negtrack1->Charge()>0) continue; // this is needed to avoid double-counting of unlike-sign
856  if(fMixEvent) {
857  if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
858  }
859 
860  }
861 
862  // back to primary vertex
863  // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
864  // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
865  SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
866  SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
867  negtrack1->GetPxPyPz(momneg1);
868 
869  // DCA between the two tracks
870  dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
871  if(dcap1n1>dcaMax) { negtrack1=0; continue; }
872 
873  // Vertexing
874  twoTrackArray1->AddAt(postrack1,0);
875  twoTrackArray1->AddAt(negtrack1,1);
876  AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
877  if(!vertexp1n1) {
878  twoTrackArray1->Clear();
879  negtrack1=0;
880  continue;
881  }
882  // 2 prong candidate
883  if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) {
884 
885  io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
886 
887  if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI) || (isLikeSign2Prong && (okD0 || okJPSI))) {
888  // add the vertex and the decay to the AOD
889  AliAODVertex *v2Prong =0x0;
890  if(!fMakeReducedRHF)v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
891  if(!isLikeSign2Prong) {
892  if(okD0) {
893  rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
895 
896  if(fMakeReducedRHF){
897  rd->DeleteRecoD();
898  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
899  }else{
900  rd->SetSecondaryVtx(v2Prong);
901  v2Prong->SetParent(rd);
902  AddRefs(v2Prong,rd,event,twoTrackArray1);
903  }
904  }
905  if(okJPSI) {
906  rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
907  if(fMakeReducedRHF){
908  rd->DeleteRecoD();
909  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
910  }else{
911  if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
912  AddRefs(v2Prong,rd,event,twoTrackArray1);
913  }
914  }
915  } else { // isLikeSign2Prong
916  rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
917  //Set selection bit for PID
919  if(fMakeReducedRHF){
920  rd->DeleteRecoD();
921  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
922  }else{
923  rd->SetSecondaryVtx(v2Prong);
924  v2Prong->SetParent(rd);
925  AddRefs(v2Prong,rd,event,twoTrackArray1);
926  }
927  }
928  }
929  // D* candidates
930  if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
931  // write references in io2Prong
932  if(fInputAOD) {
933  AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
934  } else {
935  vertexp1n1->AddDaughter(postrack1);
936  vertexp1n1->AddDaughter(negtrack1);
937  }
938  io2Prong->SetSecondaryVtx(vertexp1n1);
939  //printf("---> %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
940  // create a track from the D0
941  AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
942 
943  // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
944  for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
945 
946  if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
947 
948  if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
949 
950  if(fMixEvent) {
951  if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] ||
952  evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] ||
953  evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
954  }
955 
956  //if(iTrkSoftPi%1==0) AliDebug(1,Form(" 1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));
957 
958  trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
959  if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
960 
961  // get track from tracks array
962  trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
963  // trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
964  SetParametersAtVertex(trackPi,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkSoftPi));
965  twoTrackArrayCasc->AddAt(trackPi,0);
966  twoTrackArrayCasc->AddAt(trackD0,1);
967  if(!SelectInvMassAndPtDstarD0pi(twoTrackArrayCasc)){
968  twoTrackArrayCasc->Clear();
969  trackPi=0;
970  continue;
971  }
972 
973  AliAODVertex *vertexCasc = 0;
974 
975  if(fFindVertexForDstar) {
976  // DCA between the two tracks
977  dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
978  // Vertexing
979  vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
980  } else {
981  // assume Dstar decays at the primary vertex
982  Double_t pos[3],cov[6],chi2perNDF;
983  fV1->GetXYZ(pos);
984  fV1->GetCovMatrix(cov);
985  chi2perNDF = fV1->GetChi2toNDF();
986  vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
987  dcaCasc = 0.;
988  }
989  if(!vertexCasc) {
990  twoTrackArrayCasc->Clear();
991  trackPi=0;
992  continue;
993  }
994 
995  ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
996  if(okDstar) {
997  // add the D0 to the AOD (if not already done)
998  if(!okD0) {
999  rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
1000  if(fMakeReducedRHF){
1001  rd->DeleteRecoD();
1002  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1003  }else{
1004  AliAODVertex *v2Prong = new (verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
1005  rd->SetSecondaryVtx(v2Prong);
1006  v2Prong->SetParent(rd);
1007  AddRefs(v2Prong,rd,event,twoTrackArray1);
1008  }
1009  okD0=kTRUE; // this is done to add it only once
1010  }
1011  // add the vertex and the cascade to the AOD
1012  rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
1013  // Set selection bit for PID
1015  AliAODVertex *vCasc = 0x0;
1016  if(fMakeReducedRHF){
1017  //assign a ID to the D0 candidate, daughter of the Cascade. ID = position in the D0toKpi array
1018  UShort_t idCasc[2]={(UShort_t)trackPi->GetID(),(UShort_t)(iD0toKpi-1)};
1019  rc->SetProngIDs(2,idCasc);
1020  rc->DeleteRecoD();
1021  rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1022  }else{
1023  AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
1024  rc->SetSecondaryVtx(vCasc);
1025  vCasc->SetParent(rc);
1026  if(!fInputAOD) vCasc->AddDaughter(rd); // just to fill ref #0
1027  AddRefs(vCasc,rc,event,twoTrackArrayCasc);
1028  vCasc->AddDaughter(rd); // add the D0 (in ref #1)
1029  }
1030  }
1031  twoTrackArrayCasc->Clear();
1032  trackPi=0;
1033  if(ioCascade) {delete ioCascade; ioCascade=NULL;}
1034  delete vertexCasc; vertexCasc=NULL;
1035  } // end loop on soft pi tracks
1036 
1037  if(trackD0) {delete trackD0; trackD0=NULL;}
1038 
1039  }
1040  if(io2Prong) {delete io2Prong; io2Prong=NULL;}
1041  }
1042 
1043  twoTrackArray1->Clear();
1044  if( (!f3Prong && !f4Prong) ||
1045  (isLikeSign2Prong && !f3Prong) ) {
1046  negtrack1=0;
1047  delete vertexp1n1;
1048  continue;
1049  }
1050 
1051 
1052  // 2nd LOOP ON POSITIVE TRACKS
1053  for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
1054 
1055  if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
1056 
1057  //if(iTrkP2%1==0) AliDebug(1,Form(" 2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));
1058 
1059  // get track from tracks array
1060  postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
1061 
1062  if(postrack2->Charge()<0) continue;
1063 
1064  if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
1065 
1066  // Check single tracks cuts specific for 3 prongs
1067  if(!TESTBIT(seleFlags[iTrkP2],kBit3Prong)) continue;
1068  if(!TESTBIT(seleFlags[iTrkP1],kBit3Prong)) continue;
1069  if(!TESTBIT(seleFlags[iTrkN1],kBit3Prong)) continue;
1070 
1071  if(fMixEvent) {
1072  if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
1073  evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
1074  evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1075  }
1076 
1077  if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
1078  if(!fLikeSign3prong) continue;
1079  if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
1080  isLikeSign3Prong=kTRUE;
1081  } else { // not ok
1082  continue;
1083  }
1084  } else { // normal triplet (+-+)
1085  isLikeSign3Prong=kFALSE;
1086  if(fMixEvent) {
1087  if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
1088  evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
1089  evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1090  }
1091  }
1092 
1094  if(!TESTBIT(seleFlags[iTrkN1],kBitKaonCompat)) continue;
1095  }
1096  Bool_t okForLcTopKpi=kTRUE;
1097  Int_t pidLcStatus=3; // 3= OK as pKpi and Kpipi
1098  if(fUsePIDforLc>0){
1099  if(!TESTBIT(seleFlags[iTrkP1],kBitProtonCompat) &&
1100  !TESTBIT(seleFlags[iTrkP2],kBitProtonCompat) ){
1101  okForLcTopKpi=kFALSE;
1102  pidLcStatus=0;
1103  }
1104  if(okForLcTopKpi && fUsePIDforLc>1){
1105  okForLcTopKpi=kFALSE;
1106  pidLcStatus=0;
1107  if(TESTBIT(seleFlags[iTrkP1],kBitProtonCompat) &&
1108  TESTBIT(seleFlags[iTrkP2],kBitPionCompat) ){
1109  okForLcTopKpi=kTRUE;
1110  pidLcStatus+=1; // 1= OK as pKpi
1111  }
1112  if(TESTBIT(seleFlags[iTrkP2],kBitProtonCompat) &&
1113  TESTBIT(seleFlags[iTrkP1],kBitPionCompat) ){
1114  okForLcTopKpi=kTRUE;
1115  pidLcStatus+=2; // 2= OK as piKp
1116  }
1117  }
1118  }
1119  Bool_t okForDsToKKpi=kTRUE;
1120  if(fUseKaonPIDforDs){
1121  if(!TESTBIT(seleFlags[iTrkP1],kBitKaonCompat) &&
1122  !TESTBIT(seleFlags[iTrkP2],kBitKaonCompat) ) okForDsToKKpi=kFALSE;
1123  }
1124  // back to primary vertex
1125  // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1126  // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1127  // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1128  SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1129  SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1130  SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
1131 
1132  //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
1133 
1134  dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
1135  if(dcap2n1>dcaMax) { postrack2=0; continue; }
1136  dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
1137  if(dcap1p2>dcaMax) { postrack2=0; continue; }
1138 
1139  // check invariant mass cuts for D+,Ds,Lc
1140  massCutOK=kTRUE;
1141  if(f3Prong) {
1142  if(postrack2->Charge()>0) {
1143  threeTrackArray->AddAt(postrack1,0);
1144  threeTrackArray->AddAt(negtrack1,1);
1145  threeTrackArray->AddAt(postrack2,2);
1146  } else {
1147  threeTrackArray->AddAt(negtrack1,0);
1148  threeTrackArray->AddAt(postrack1,1);
1149  threeTrackArray->AddAt(postrack2,2);
1150  }
1152  postrack2->GetPxPyPz(mompos2);
1153  Double_t pxDau[3]={mompos1[0],momneg1[0],mompos2[0]};
1154  Double_t pyDau[3]={mompos1[1],momneg1[1],mompos2[1]};
1155  Double_t pzDau[3]={mompos1[2],momneg1[2],mompos2[2]};
1156  // massCutOK = SelectInvMassAndPt3prong(threeTrackArray);
1157  massCutOK = SelectInvMassAndPt3prong(pxDau,pyDau,pzDau,pidLcStatus);
1158  }
1159  }
1160 
1161  if(f3Prong && !massCutOK) {
1162  threeTrackArray->Clear();
1163  if(!f4Prong) {
1164  postrack2=0;
1165  continue;
1166  }
1167  }
1168 
1169  // Vertexing
1170  twoTrackArray2->AddAt(postrack2,0);
1171  twoTrackArray2->AddAt(negtrack1,1);
1172  AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1173  if(!vertexp2n1) {
1174  twoTrackArray2->Clear();
1175  postrack2=0;
1176  continue;
1177  }
1178 
1179  // 3 prong candidates
1180  if(f3Prong && massCutOK) {
1181 
1182  AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1183  io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,okForLcTopKpi,okForDsToKKpi,ok3Prong);
1184  if(ok3Prong) {
1185  AliAODVertex *v3Prong=0x0;
1186  if(!fMakeReducedRHF)v3Prong = new (verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1187  if(!isLikeSign3Prong) {
1188  rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1189  // Set selection bit for PID
1193  if(fMakeReducedRHF){
1194  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1195  ((AliAODRecoDecayHF3Prong*)rd)->DeleteRecoD();
1196  }else{
1197  v3Prong = new (verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1198  rd->SetSecondaryVtx(v3Prong);
1199  v3Prong->SetParent(rd);
1200  AddRefs(v3Prong,rd,event,threeTrackArray);
1201  }
1202  } else { // isLikeSign3Prong
1203  if(fLikeSign3prong){
1204  rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1205  // Set selection bit for PID
1209  if(fMakeReducedRHF){
1210  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1211  ((AliAODRecoDecayHF3Prong*)rd)->DeleteRecoD();
1212  }else{
1213  rd->SetSecondaryVtx(v3Prong);
1214  v3Prong->SetParent(rd);
1215  AddRefs(v3Prong,rd,event,threeTrackArray);
1216  }
1217  }
1218  }
1219 
1220  }
1221  if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1222  if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1223  }
1224 
1225  // 4 prong candidates
1226  if(f4Prong
1227  // don't make 4 prong with like-sign pairs and triplets
1228  && !isLikeSign2Prong && !isLikeSign3Prong
1229  // track-to-track dca cuts already now
1230  && dcap1n1 < fCutsD0toKpipipi->GetDCACut()
1231  && dcap2n1 < fCutsD0toKpipipi->GetDCACut()) {
1232  // back to primary vertex
1233  // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1234  // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1235  // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1236  SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1237  SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1238  SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
1239 
1240  // Vertexing for these 3 (can be taken from above?)
1241  threeTrackArray->AddAt(postrack1,0);
1242  threeTrackArray->AddAt(negtrack1,1);
1243  threeTrackArray->AddAt(postrack2,2);
1244  AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1245 
1246  // 3rd LOOP ON NEGATIVE TRACKS (for 4 prong)
1247  for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
1248 
1249  if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1250 
1251  //if(iTrkN2%1==0) AliDebug(1,Form(" 3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
1252 
1253  // get track from tracks array
1254  negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1255 
1256  if(negtrack2->Charge()>0) continue;
1257 
1258  if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1259  if(fMixEvent){
1260  if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
1261  evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
1262  evtNumber[iTrkP2]==evtNumber[iTrkN2] ||
1263  evtNumber[iTrkP1]==evtNumber[iTrkN1] ||
1264  evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
1265  evtNumber[iTrkN1]==evtNumber[iTrkP2]) continue;
1266  }
1267 
1268  // back to primary vertex
1269  // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1270  // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1271  // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1272  // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1273  SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1274  SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1275  SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
1276  SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
1277 
1278  dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1279  if(dcap1n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
1280  dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1281  if(dcap2n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
1282 
1283 
1284  fourTrackArray->AddAt(postrack1,0);
1285  fourTrackArray->AddAt(negtrack1,1);
1286  fourTrackArray->AddAt(postrack2,2);
1287  fourTrackArray->AddAt(negtrack2,3);
1288 
1289  // check invariant mass cuts for D0
1290  massCutOK=kTRUE;
1292  massCutOK = SelectInvMassAndPt4prong(fourTrackArray);
1293 
1294  if(!massCutOK) {
1295  fourTrackArray->Clear();
1296  negtrack2=0;
1297  continue;
1298  }
1299 
1300  // Vertexing
1301  AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
1302  io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
1303  if(ok4Prong) {
1304  rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
1305  if(fMakeReducedRHF){
1306  rd->DeleteRecoD();
1307  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1308  }else{
1309  AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
1310  rd->SetSecondaryVtx(v4Prong);
1311  v4Prong->SetParent(rd);
1312  AddRefs(v4Prong,rd,event,fourTrackArray);
1313  }
1314  }
1315 
1316  if(io4Prong) {delete io4Prong; io4Prong=NULL;}
1317  if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;}
1318  fourTrackArray->Clear();
1319  negtrack2 = 0;
1320 
1321  } // end loop on negative tracks
1322 
1323  threeTrackArray->Clear();
1324  delete vertexp1n1p2;
1325 
1326  }
1327 
1328  postrack2 = 0;
1329  delete vertexp2n1;
1330 
1331  } // end 2nd loop on positive tracks
1332 
1333  twoTrackArray2->Clear();
1334 
1335  // 2nd LOOP ON NEGATIVE TRACKS (for 3 prong -+-)
1336  for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
1337 
1338  if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1339 
1340  //if(iTrkN2%1==0) AliDebug(1,Form(" 2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
1341 
1342  // get track from tracks array
1343  negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1344 
1345  if(negtrack2->Charge()>0) continue;
1346 
1347  if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1348 
1349  // Check single tracks cuts specific for 3 prongs
1350  if(!TESTBIT(seleFlags[iTrkN2],kBit3Prong)) continue;
1351  if(!TESTBIT(seleFlags[iTrkP1],kBit3Prong)) continue;
1352  if(!TESTBIT(seleFlags[iTrkN1],kBit3Prong)) continue;
1353 
1354  if(fMixEvent) {
1355  if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
1356  evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
1357  evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1358  }
1359 
1360  if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
1361  if(!fLikeSign3prong) continue;
1362  if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
1363  isLikeSign3Prong=kTRUE;
1364  } else { // not ok
1365  continue;
1366  }
1367  } else { // normal triplet (-+-)
1368  isLikeSign3Prong=kFALSE;
1369  }
1370 
1372  if(!TESTBIT(seleFlags[iTrkP1],kBitKaonCompat)) continue;
1373  }
1374  Bool_t okForLcTopKpi=kTRUE;
1375  Int_t pidLcStatus=3; // 3= OK as pKpi and Kpipi
1376  if(fUsePIDforLc>0){
1377  if(!TESTBIT(seleFlags[iTrkN1],kBitProtonCompat) &&
1378  !TESTBIT(seleFlags[iTrkN2],kBitProtonCompat) ){
1379  okForLcTopKpi=kFALSE;
1380  pidLcStatus=0;
1381  }
1382  if(okForLcTopKpi && fUsePIDforLc>1){
1383  okForLcTopKpi=kFALSE;
1384  pidLcStatus=0;
1385  if(TESTBIT(seleFlags[iTrkN1],kBitProtonCompat) &&
1386  TESTBIT(seleFlags[iTrkN2],kBitPionCompat) ){
1387  okForLcTopKpi=kTRUE;
1388  pidLcStatus+=1; // 1= OK as pKpi
1389  }
1390  if(TESTBIT(seleFlags[iTrkN2],kBitProtonCompat) &&
1391  TESTBIT(seleFlags[iTrkN1],kBitPionCompat) ){
1392  okForLcTopKpi=kTRUE;
1393  pidLcStatus+=2; // 2= OK as piKp
1394  }
1395  }
1396  }
1397  Bool_t okForDsToKKpi=kTRUE;
1398  if(fUseKaonPIDforDs){
1399  if(!TESTBIT(seleFlags[iTrkN1],kBitKaonCompat) &&
1400  !TESTBIT(seleFlags[iTrkN2],kBitKaonCompat) ) okForDsToKKpi=kFALSE;
1401  }
1402 
1403  // back to primary vertex
1404  // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1405  // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1406  // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1407  SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1408  SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1409  SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
1410  //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
1411 
1412  dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1413  if(dcap1n2>dcaMax) { negtrack2=0; continue; }
1414  dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1415  if(dcan1n2>dcaMax) { negtrack2=0; continue; }
1416 
1417  threeTrackArray->AddAt(negtrack1,0);
1418  threeTrackArray->AddAt(postrack1,1);
1419  threeTrackArray->AddAt(negtrack2,2);
1420 
1421  // check invariant mass cuts for D+,Ds,Lc
1422  massCutOK=kTRUE;
1424  negtrack2->GetPxPyPz(momneg2);
1425  Double_t pxDau[3]={momneg1[0],mompos1[0],momneg2[0]};
1426  Double_t pyDau[3]={momneg1[1],mompos1[1],momneg2[1]};
1427  Double_t pzDau[3]={momneg1[2],mompos1[2],momneg2[2]};
1428  // massCutOK = SelectInvMassAndPt3prong(threeTrackArray);
1429  massCutOK = SelectInvMassAndPt3prong(pxDau,pyDau,pzDau,pidLcStatus);
1430  }
1431  if(!massCutOK) {
1432  threeTrackArray->Clear();
1433  negtrack2=0;
1434  continue;
1435  }
1436 
1437  // Vertexing
1438  twoTrackArray2->AddAt(postrack1,0);
1439  twoTrackArray2->AddAt(negtrack2,1);
1440 
1441  AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1442  if(!vertexp1n2) {
1443  twoTrackArray2->Clear();
1444  negtrack2=0;
1445  continue;
1446  }
1447 
1448  if(f3Prong) {
1449  AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1450  io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,okForLcTopKpi,okForDsToKKpi,ok3Prong);
1451  if(ok3Prong) {
1452  AliAODVertex *v3Prong = 0x0;
1453  if(!fMakeReducedRHF) v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1454  if(!isLikeSign3Prong) {
1455  rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1459  if(fMakeReducedRHF){
1460  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1461  ((AliAODRecoDecayHF3Prong*)rd)->DeleteRecoD();
1462  }else{
1463  rd->SetSecondaryVtx(v3Prong);
1464  v3Prong->SetParent(rd);
1465  AddRefs(v3Prong,rd,event,threeTrackArray);
1466  }
1467  } else { // isLikeSign3Prong
1468  if(fLikeSign3prong){
1469  rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1473  if(fMakeReducedRHF){
1474  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1475  ((AliAODRecoDecayHF3Prong*)rd)->DeleteRecoD();
1476  }else{
1477  rd->SetSecondaryVtx(v3Prong);
1478  v3Prong->SetParent(rd);
1479  AddRefs(v3Prong,rd,event,threeTrackArray);
1480  }
1481  }
1482 
1483  }
1484  }
1485  if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1486  if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1487  }
1488  threeTrackArray->Clear();
1489  negtrack2 = 0;
1490  delete vertexp1n2;
1491 
1492  } // end 2nd loop on negative tracks
1493 
1494  twoTrackArray2->Clear();
1495 
1496  negtrack1 = 0;
1497  delete vertexp1n1;
1498  } // end 1st loop on negative tracks
1499 
1500  postrack1 = 0;
1501  } // end 1st loop on positive tracks
1502 
1503 
1504  // AliDebug(1,Form(" Total HF vertices in event = %d;",
1505  // (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
1506  if(fD0toKpi) {
1507  AliDebug(1,Form(" D0->Kpi in event = %d;",
1508  (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
1509  }
1510  if(fJPSItoEle) {
1511  AliDebug(1,Form(" JPSI->ee in event = %d;",
1512  (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
1513  }
1514  if(f3Prong) {
1515  AliDebug(1,Form(" Charm->3Prong in event = %d;",
1516  (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
1517  }
1518  if(f4Prong) {
1519  AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
1520  (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
1521  }
1522  if(fDstar) {
1523  AliDebug(1,Form(" D*->D0pi in event = %d;\n",
1524  (Int_t)aodDstarTClArr->GetEntriesFast()));
1525  }
1526  if(fCascades){
1527  AliDebug(1,Form(" cascades -> v0 + track in event = %d;\n",
1528  (Int_t)aodCascadesTClArr->GetEntriesFast()));
1529  }
1530  if(fLikeSign) {
1531  AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
1532  (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
1533  }
1534  if(fLikeSign3prong && f3Prong) {
1535  AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
1536  (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
1537  }
1538 
1539 
1540  twoTrackArray1->Delete(); delete twoTrackArray1;
1541  twoTrackArray2->Delete(); delete twoTrackArray2;
1542  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1543  twoTrackArrayV0->Delete(); delete twoTrackArrayV0;
1544  threeTrackArray->Clear();
1545  threeTrackArray->Delete(); delete threeTrackArray;
1546  fourTrackArray->Delete(); delete fourTrackArray;
1547  delete [] seleFlags; seleFlags=NULL;
1548  if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
1549  tracksAtVertex.Delete();
1550 
1551  if(fInputAOD) {
1552  seleTrksArray.Delete();
1553  if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
1554  }
1555 
1556 
1557  //printf("Trks: total %d sele %d\n",fnTrksTotal,fnSeleTrksTotal);
1558 
1559  return;
1560 }
1561 //----------------------------------------------------------------------------
1563  const AliVEvent *event,
1564  const TObjArray *trkArray) const
1565 {
1568  //AliCodeTimerAuto("",0);
1569 
1570  if(fInputAOD) {
1571  AddDaughterRefs(v,event,trkArray);
1572  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1573  }
1574 
1575  /*
1576  rd->SetListOfCutsRef((TList*)fListOfCuts);
1577  //fListOfCuts->Print();
1578  cout<<fListOfCuts<<endl;
1579  TList *l=(TList*)rd->GetListOfCuts();
1580  cout<<l<<endl;
1581  if(l) {l->Print(); }else{printf("error\n");}
1582  */
1583 
1584  return;
1585 }
1586 //---------------------------------------------------------------------------
1588  const AliVEvent *event,
1589  const TObjArray *trkArray) const
1590 {
1592  //AliCodeTimerAuto("",0);
1593 
1594  Int_t nDg = v->GetNDaughters();
1595  TObject *dg = 0;
1596  if(nDg) dg = v->GetDaughter(0);
1597 
1598  if(dg) return; // daughters already added
1599 
1600  Int_t nTrks = trkArray->GetEntriesFast();
1601 
1602  AliExternalTrackParam *track = 0;
1603  AliAODTrack *aodTrack = 0;
1604  Int_t id;
1605 
1606  for(Int_t i=0; i<nTrks; i++) {
1607  track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1608  id = (Int_t)track->GetID();
1609  //printf("---> %d\n",id);
1610  if(id<0) continue; // this track is a AliAODRecoDecay
1611  aodTrack = dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[id]));
1612  if(!aodTrack) AliFatal("Not a standard AOD");
1613  v->AddDaughter(aodTrack);
1614  }
1615 
1616  return;
1617 }
1618 //----------------------------------------------------------------------------
1620 {
1623  //
1624  //AliCodeTimerAuto("",0);
1625 
1626 
1627  TClonesArray *inputArray=(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
1628  if(!inputArray) return;
1629 
1630  AliAODTrack *track = 0;
1631  AliAODVertex *vertex = 0;
1632 
1633  Bool_t needtofix=kFALSE;
1634  for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1635  vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1636  for(Int_t id=0; id<vertex->GetNDaughters(); id++) {
1637  track = (AliAODTrack*)vertex->GetDaughter(id);
1638  if(!track->GetStatus()) needtofix=kTRUE;
1639  }
1640  if(needtofix) break;
1641  }
1642 
1643  if(!needtofix) return;
1644 
1645 
1646  printf("Fixing references\n");
1647 
1648  fAODMapSize = 100000;
1649  fAODMap = new Int_t[fAODMapSize];
1650  memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
1651 
1652  for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
1653  track = dynamic_cast<AliAODTrack*>(aod->GetTrack(i));
1654  if(!track) AliFatal("Not a standard AOD");
1655 
1656  // skip pure ITS SA tracks
1657  if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
1658 
1659  // skip tracks without ITS
1660  if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
1661 
1662  // TEMPORARY: check that the cov matrix is there
1663  Double_t covtest[21];
1664  if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1665  //
1666 
1667  Int_t ind = (Int_t)track->GetID();
1668  if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
1669  }
1670 
1671 
1672  Int_t ids[4]={-1,-1,-1,-1};
1673  for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1674  Bool_t cascade=kFALSE;
1675  vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1676  Int_t id=0;
1677  Int_t nDgs = vertex->GetNDaughters();
1678  for(id=0; id<nDgs; id++) {
1679  track = (AliAODTrack*)vertex->GetDaughter(id);
1680  if(track->Charge()==0) {cascade=kTRUE; continue;} // cascade
1681  ids[id]=(Int_t)track->GetID();
1682  vertex->RemoveDaughter(track);
1683  }
1684  if(cascade) continue;
1685  for(id=0; id<nDgs; id++) {
1686  if (ids[id]>-1 && ids[id] < fAODMapSize) {
1687  track = dynamic_cast<AliAODTrack*>(aod->GetTrack(fAODMap[ids[id]]));
1688  if(!track) AliFatal("Not a standard AOD");
1689  vertex->AddDaughter(track);
1690  }
1691  }
1692 
1693  }
1694 
1695  return;
1696 }
1697 
1698 //----------------------------------------------------------------------------
1699 AliAODTrack* AliAnalysisVertexingHF::GetProng(AliVEvent *event,AliAODRecoDecayHF *rd,Int_t iprong){
1700  if(!fAODMap)MapAODtracks(event);
1701  return (AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(iprong)]);
1702 }
1703 
1704 
1705 //----------------------------------------------------------------------------
1707  // method to retrieve daughters from trackID and reconstruct secondary vertex
1708  // save the TRefs to the candidate AliAODRecoDecayHF3Prong rd
1709  // and fill on-the-fly the data member of rd
1710  if(rd->GetIsFilled()!=0)return kTRUE;//if 0: reduced dAOD. skip if rd is already filled (1: standard dAOD, 2 already refilled)
1711  if(!fAODMap)MapAODtracks(event);//fill the AOD index map if it is not yet done
1712  TObjArray *threeTrackArray = new TObjArray(3);
1713 
1714  AliAODTrack *track1 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(0)]);//retrieve daughter from the trackID through the AOD index map
1715  if(!track1)return kFALSE;
1716  AliAODTrack *track2 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(1)]);
1717  if(!track2)return kFALSE;
1718  AliESDtrack *postrack1 = 0;
1719  AliESDtrack *negtrack1 = 0;
1720  postrack1 = new AliESDtrack(track1);
1721  negtrack1 = new AliESDtrack(track2);
1722 
1723  // DCA between the two tracks
1724  Double_t xdummy, ydummy;
1725  fBzkG = (Double_t)event->GetMagneticField();
1726  Double_t dca12 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
1727 
1728  const AliVVertex *vprimary = event->GetPrimaryVertex();
1729  Double_t pos[3];
1730  Double_t cov[6];
1731  vprimary->GetXYZ(pos);
1732  vprimary->GetCovarianceMatrix(cov);
1733  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1734  fV1->GetCovMatrix(cov);
1735  if(!fVertexerTracks)fVertexerTracks=new AliVertexerTracks(fBzkG);
1736 
1737  AliAODTrack *track3 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(2)]);
1738  if(!track3)return kFALSE;
1739  AliESDtrack *esdt3 = new AliESDtrack(track3);
1740 
1741  Double_t dca2;
1742  Double_t dca3;
1743  threeTrackArray->AddAt(postrack1,0);
1744  threeTrackArray->AddAt(negtrack1,1);
1745  threeTrackArray->AddAt(esdt3,2);
1746  dca2 = esdt3->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
1747  dca3 = esdt3->GetDCA(postrack1,fBzkG,xdummy,ydummy);
1748  Double_t dispersion;
1749 
1750  AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray, dispersion);
1751  if (!secVert3PrAOD) {
1752  threeTrackArray->Clear();
1753  threeTrackArray->Delete(); delete threeTrackArray;
1754  delete fV1; fV1=0;
1755  delete postrack1; postrack1=NULL;
1756  delete negtrack1; negtrack1=NULL;
1757  delete esdt3; esdt3=NULL;
1758  return kFALSE;
1759  }
1760 
1761  rd->SetNProngs();
1762  Double_t vtxRec=rd->GetDist12toPrim();
1763  Double_t vertexp2n1=rd->GetDist23toPrim();
1764  rd= Make3Prong(threeTrackArray, event, secVert3PrAOD,dispersion, vtxRec, vertexp2n1, dca12, dca2, dca3, rd);
1765  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1766  rd->SetIsFilled(2);
1767  threeTrackArray->Clear();
1768  threeTrackArray->Delete(); delete threeTrackArray;
1769  delete fV1; fV1=0;
1770  delete postrack1; postrack1=NULL;
1771  delete negtrack1; negtrack1=NULL;
1772  delete esdt3; esdt3=NULL;
1773  return kTRUE;
1774 }
1775 //___________________________
1777  // method to retrieve daughters from trackID and reconstruct secondary vertex
1778  // save the TRefs to the candidate AliAODRecoDecayHF2Prong rd
1779  // and fill on-the-fly the data member of rd
1780  if(rd->GetIsFilled()!=0)return kTRUE;//if 0: reduced dAOD. skip if rd is already filled (1:standard dAOD, 2 already refilled)
1781  if(!fAODMap)MapAODtracks(event);//fill the AOD index map if it is not yet done
1782 
1783  Double_t dispersion;
1784  TObjArray *twoTrackArray1 = new TObjArray(2);
1785 
1786  AliAODTrack *track1 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(0)]);//retrieve daughter from the trackID through the AOD index map
1787  if(!track1)return kFALSE;
1788  AliAODTrack *track2 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(1)]);
1789  if(!track2)return kFALSE;
1790 
1791  AliESDtrack *esdt1 = 0;
1792  AliESDtrack *esdt2 = 0;
1793  esdt1 = new AliESDtrack(track1);
1794  esdt2 = new AliESDtrack(track2);
1795 
1796  twoTrackArray1->AddAt(esdt1,0);
1797  twoTrackArray1->AddAt(esdt2,1);
1798  // DCA between the two tracks
1799  Double_t xdummy, ydummy;
1800  fBzkG = (Double_t)event->GetMagneticField();
1801  Double_t dca12 = esdt1->GetDCA(esdt2,fBzkG,xdummy,ydummy);
1802  const AliVVertex *vprimary = event->GetPrimaryVertex();
1803  Double_t pos[3];
1804  Double_t cov[6];
1805  vprimary->GetXYZ(pos);
1806  vprimary->GetCovarianceMatrix(cov);
1807  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1808  fV1->GetCovMatrix(cov);
1809  if(!fVertexerTracks)fVertexerTracks=new AliVertexerTracks(fBzkG);
1810 
1811 
1812  AliAODVertex *vtxRec = ReconstructSecondaryVertex(twoTrackArray1, dispersion);
1813  if(!vtxRec) {
1814  twoTrackArray1->Clear();
1815  twoTrackArray1->Delete(); delete twoTrackArray1;
1816  delete fV1; fV1=0;
1817  delete esdt1; esdt1=NULL;
1818  delete esdt2; esdt2=NULL;
1819  return kFALSE; }
1820  Bool_t okD0=kFALSE;
1821  Bool_t okJPSI=kFALSE;
1822  Bool_t okD0FromDstar=kFALSE;
1823  Bool_t refill =kTRUE;
1824  rd->SetNProngs();
1825  rd= Make2Prong(twoTrackArray1, event, vtxRec, dca12, okD0, okJPSI, okD0FromDstar,refill,rd);
1826  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1827  rd->SetIsFilled(2);
1828  delete fV1; fV1=0;
1829  twoTrackArray1->Clear();
1830  twoTrackArray1->Delete(); delete twoTrackArray1;
1831  delete esdt1; esdt1=NULL;
1832  delete esdt2; esdt2=NULL;
1833  return kTRUE;
1834 }
1835 //----------------------------------------------------------------------------
1837  // method to retrieve daughters from trackID
1838  // and fill on-the-fly the data member of rCasc and their AliAODRecoDecayHF2Prong daughters
1839  if(rCasc->GetIsFilled()!=0) return kTRUE;//if 0: reduced dAOD. skip if rd is already filled (1: standard dAOD, 2: already refilled)
1840  if(!fAODMap)MapAODtracks(event);//fill the AOD index map if it is not yet done
1841  TObjArray *twoTrackArrayCasc = new TObjArray(2);
1842 
1843  AliAODTrack *trackB =(AliAODTrack*)event->GetTrack(fAODMap[rCasc->GetProngID(0)]);//retrieve bachelor from the trackID through the AOD index map
1844  if(!trackB)return kFALSE;
1845 
1846  AliNeutralTrackParam *trackV0=NULL;
1847  AliAODv0 *v0 =NULL;
1848  AliAODRecoDecayHF2Prong *trackD0=NULL;
1849  rCasc->SetNProngs();
1850 
1851  if(DStar){
1852  TClonesArray *inputArrayD0=(TClonesArray*)event->GetList()->FindObject("D0toKpi");
1853  if(!inputArrayD0) return kFALSE;
1854  trackD0=(AliAODRecoDecayHF2Prong*)inputArrayD0->At(rCasc->GetProngID(1));
1855  if(!trackD0)return kFALSE;
1856  if(!FillRecoCand(event,trackD0)) return kFALSE; //fill missing info of the corresponding D0 daughter
1857 
1858  trackV0 = new AliNeutralTrackParam(trackD0);
1859 
1860  }else{//is a V0 candidate
1861  v0 = ((AliAODEvent*)event)->GetV0(rCasc->GetProngID(1));
1862  if(!v0) return kFALSE;
1863  // Define the V0 (neutral) track
1864  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
1865  if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
1866  }
1867 
1868  AliESDtrack *esdB = new AliESDtrack(trackB);
1869 
1870  twoTrackArrayCasc->AddAt(esdB,0);
1871  twoTrackArrayCasc->AddAt(trackV0,1);
1872 
1873  fBzkG = (Double_t)event->GetMagneticField();
1874  const AliVVertex *vprimary = event->GetPrimaryVertex();
1875 
1876  Double_t pos[3];
1877  Double_t cov[6];
1878  vprimary->GetXYZ(pos);
1879  vprimary->GetCovarianceMatrix(cov);
1880  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1881  fV1->GetCovMatrix(cov);
1882 
1883  Double_t dca = 0.;
1884  AliAODVertex *vtxCasc = 0x0;
1885  Double_t chi2perNDF = fV1->GetChi2toNDF();
1886  if (recoSecVtx) {
1887  Double_t dispersion, xdummy, ydummy;
1888  dca = esdB->GetDCA(trackV0,fBzkG,xdummy,ydummy);
1889  if (!fVertexerTracks) fVertexerTracks = new AliVertexerTracks(fBzkG);
1890  vtxCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
1891  } else {
1892  vtxCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
1893  }
1894  if(!vtxCasc) {
1895  twoTrackArrayCasc->Clear();
1896  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1897  delete fV1; fV1=0;
1898  delete esdB; esdB=NULL;
1899  delete vtxCasc;vtxCasc=NULL;
1900  trackB=NULL;
1901  delete trackV0; trackV0=NULL;
1902  if(!DStar){
1903  v0=NULL;
1904  }
1905  return kFALSE;
1906  }
1907  vtxCasc->SetParent(rCasc);
1908  rCasc->SetSecondaryVtx(vtxCasc);
1909  AddDaughterRefs(vtxCasc,(AliAODEvent*)event,twoTrackArrayCasc);
1910  if(DStar)vtxCasc->AddDaughter(trackD0);
1911  else vtxCasc->AddDaughter(v0);
1912  rCasc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1913 
1914  Bool_t refill =kTRUE;
1915 
1916  Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1917  // propagate tracks to secondary vertex, to compute inv. mass
1918  esdB->PropagateToDCA(vtxCasc,fBzkG,kVeryBig);
1919  trackV0->PropagateToDCA(vtxCasc,fBzkG,kVeryBig);
1920  Double_t momentum[3];
1921  esdB->GetPxPyPz(momentum);
1922  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1923  trackV0->GetPxPyPz(momentum);
1924  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1925 
1926  AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArrayCasc,event);
1927  if(!primVertexAOD){
1928  delete fV1; fV1=0;
1929  delete vtxCasc; vtxCasc=NULL;
1930  twoTrackArrayCasc->Clear();
1931  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1932  delete esdB; esdB=NULL;
1933  delete trackV0; trackV0=NULL;
1934  if(!DStar)v0=NULL;
1935  return kFALSE;
1936  }
1937  Double_t d0z0[2],covd0z0[3];
1938  esdB->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1939  d0[0] = d0z0[0];
1940  d0err[0] = TMath::Sqrt(covd0z0[0]);
1941  trackV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1942  d0[1] = d0z0[0];
1943  d0err[1] = TMath::Sqrt(covd0z0[0]);
1944  rCasc->SetPxPyPzProngs(2,px,py,pz);
1945  rCasc->SetDCA(dca);
1946  rCasc->Setd0Prongs(2,d0);
1947  rCasc->Setd0errProngs(2,d0err);
1948  rCasc->SetOwnPrimaryVtx(primVertexAOD);
1949  rCasc->SetCharge(esdB->Charge());
1950  // get PID info from ESD
1951  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1952  if(esdB->GetStatus()&AliESDtrack::kESDpid) esdB->GetESDpid(esdpid0);
1953  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1954  Double_t esdpid[10];
1955  for(Int_t i=0;i<5;i++) {
1956  esdpid[i] = esdpid0[i];
1957  esdpid[5+i] = esdpid1[i];
1958  }
1959  rCasc->SetPID(2,esdpid);
1960  rCasc->SetIsFilled(2);
1961 
1962 
1963  delete fV1; fV1=0;
1964  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1965  twoTrackArrayCasc->Clear();
1966  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1967  delete esdB; esdB=NULL;
1968  delete trackV0; trackV0=NULL;
1969  if(!DStar)v0=NULL;
1970  return kTRUE;
1971 }
1972 //---------------------------------------------------------------------------
1974 {
1979 
1980 
1981  if (rc->GetIsFilled()==0) return kFALSE; // Reduced dAOD with candidates not yet refilled
1982  if (!fAODMap) MapAODtracks(event); // Fill the AOD index map if it is not yet done
1983 
1984 
1985  // - Get bachelor and V0 from the cascade
1986  AliAODTrack *trackB = dynamic_cast<AliAODTrack*> (rc->GetBachelor());
1987  if (!trackB) return kFALSE;
1988 
1989  AliAODv0 *v0 = dynamic_cast<AliAODv0*> (rc->Getv0());
1990  if (!v0) return kFALSE;
1991 
1992 
1993  // - Cast bachelor (AOD->ESD) and V0 (AliAODv0->AliNeutralTrackParam)
1994  AliESDtrack *esdB = new AliESDtrack(trackB);
1995  if (!esdB) return kFALSE;
1996 
1997  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
1998  if (!trackVV0) return kFALSE;
1999 
2000  AliNeutralTrackParam *trackV0 = new AliNeutralTrackParam(trackVV0);
2001  if (!trackV0) return kFALSE;
2002 
2003  TObjArray *twoTrackArrayCasc = new TObjArray(2);
2004  twoTrackArrayCasc->AddAt(esdB, 0);
2005  twoTrackArrayCasc->AddAt(trackV0, 1);
2006 
2007 
2008  Double_t dispersion, xdummy, ydummy, pos[3], cov[6];
2009 
2010 
2011  // - Some ingredients will be needed:
2012  // magnetic field
2013  fBzkG = (Double_t) event->GetMagneticField();
2014 
2015  // primary vertex
2016  const AliVVertex *vprimary = event->GetPrimaryVertex();
2017  vprimary->GetXYZ(pos);
2018  vprimary->GetCovarianceMatrix(cov);
2019  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2020 
2021  // vertexer
2022  if (!fVertexerTracks) fVertexerTracks = new AliVertexerTracks(fBzkG);
2023 
2024 
2025  // - Compute the DCA and the decay vertex
2026  Double_t dca = esdB->GetDCA(trackV0, fBzkG, xdummy, ydummy);
2027  AliAODVertex *vtxCasc = ReconstructSecondaryVertex(twoTrackArrayCasc, dispersion, kFALSE);
2028 
2029  if (!vtxCasc) {
2030  twoTrackArrayCasc->Clear();
2031  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
2032  delete fV1; fV1=0;
2033  delete esdB; esdB=0;
2034  delete vtxCasc; vtxCasc=0;
2035  delete trackV0; trackV0=0;
2036  trackB=0; v0=0;
2037  return kFALSE;
2038  }
2039 
2040 
2041  // - Linking the cascade with the new secondary vertex
2042  vtxCasc->SetParent(rc);
2043  rc->SetSecondaryVtx(vtxCasc);
2044  AddDaughterRefs(vtxCasc, (AliAODEvent*)event, twoTrackArrayCasc);
2045  vtxCasc->AddDaughter(v0);
2046  rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
2047 
2048 
2049  // - Propagate tracks to secondary vertex, to get momenta
2050  Double_t momentum[3], px[2], py[2], pz[2], d0[2], d0err[2], d0z0[2], covd0z0[3];
2051  esdB->PropagateToDCA(vtxCasc, fBzkG, kVeryBig);
2052  trackV0->PropagateToDCA(vtxCasc, fBzkG, kVeryBig);
2053  esdB->GetPxPyPz(momentum);
2054  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2055  trackV0->GetPxPyPz(momentum);
2056  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2057 
2058 
2059  // - Propagate tracks to primary vertex, to get impact parameters
2060  AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArrayCasc, event);
2061  if (!primVertexAOD) {
2062  twoTrackArrayCasc->Clear();
2063  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
2064  delete fV1; fV1=0;
2065  delete esdB; esdB=0;
2066  delete vtxCasc; vtxCasc=0;
2067  delete trackV0; trackV0=0;
2068  trackB=0; v0=0;
2069  return kFALSE;
2070  }
2071 
2072  esdB->PropagateToDCA(primVertexAOD, fBzkG, kVeryBig, d0z0, covd0z0);
2073  d0[0] = d0z0[0];
2074  d0err[0] = TMath::Sqrt(covd0z0[0]);
2075  trackV0->PropagateToDCA(primVertexAOD, fBzkG, kVeryBig, d0z0, covd0z0);
2076  d0[1] = d0z0[0];
2077  d0err[1] = TMath::Sqrt(covd0z0[0]);
2078 
2079 
2080  // - Get PID info from ESD
2081  Double_t esdpid[10];
2082  Double_t esdpid0[5] = {0., 0., 0., 0., 0.};
2083  Double_t esdpid1[5] = {0., 0., 0., 0., 0.};
2084  if (esdB->GetStatus()&AliESDtrack::kESDpid) esdB->GetESDpid(esdpid0);
2085  for (Int_t ipid=0; ipid<5; ipid++) {
2086  esdpid[ipid] = esdpid0[ipid];
2087  esdpid[5+ipid] = esdpid1[ipid];
2088  }
2089 
2090 
2091  // - Set data members
2092  rc->SetOwnPrimaryVtx(primVertexAOD);
2093  rc->SetDCA(dca);
2094  rc->SetPxPyPzProngs(2, px, py, pz);
2095  rc->Setd0Prongs(2, d0);
2096  rc->Setd0errProngs(2, d0err);
2097  rc->SetCharge(esdB->Charge());
2098  rc->SetPID(2, esdpid);
2099 
2100  rc->SetIsFilled(2); // To clean the newly reconstructed secondary vertex with the CleanUp task
2101 
2102 
2103  twoTrackArrayCasc->Clear();
2104  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
2105  delete fV1; fV1=0;
2106  delete primVertexAOD; primVertexAOD=0;
2107  delete esdB; esdB=0;
2108  delete trackV0; trackV0=0;
2109  trackB=0; v0=0;
2110 
2111  return kTRUE;
2112 }
2113 
2114 //---------------------------------------------------------------------------
2116  TObjArray *twoTrackArray,AliVEvent *event,
2117  AliAODVertex *secVert,
2118  AliAODRecoDecayHF2Prong *rd2Prong,
2119  Double_t dca,
2120  Bool_t &okDstar)
2121 {
2124  //AliCodeTimerAuto("",0);
2125  UInt_t ntref=TProcessID::GetObjectCount();
2126  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2127  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2128  }
2129 
2130  okDstar = kFALSE;
2131 
2132  Bool_t dummy1,dummy2,dummy3;
2133 
2134  // We use Make2Prong to construct the AliAODRecoCascadeHF
2135  // (which inherits from AliAODRecoDecayHF2Prong)
2136  AliAODRecoCascadeHF *theCascade =
2137  (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
2138  dummy1,dummy2,dummy3);
2139  if(!theCascade) return 0x0;
2140 
2141  // charge
2142  AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2143  theCascade->SetCharge(trackPi->Charge());
2144 
2145  //--- selection cuts
2146  //
2147  AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
2148  if(fInputAOD){
2149  Int_t idSoftPi=(Int_t)trackPi->GetID();
2150  if (idSoftPi > -1 && idSoftPi < fAODMapSize) {
2151  AliAODTrack* trackPiAOD=dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[idSoftPi]));
2152  if(!trackPiAOD) AliFatal("Not a standard AOD");
2153  tmpCascade->GetSecondaryVtx()->AddDaughter(trackPiAOD);
2154  }
2155  }else{
2156  tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
2157  }
2158  tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
2159 
2160  AliAODVertex *primVertexAOD=0;
2162  // take event primary vertex
2163  primVertexAOD = PrimaryVertex();
2164  tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
2165  rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
2166  }
2167  // select D*->D0pi
2168  if(fDstar) {
2170  if(okDstar) theCascade->SetSelectionBit(AliRDHFCuts::kDstarCuts);
2171  }
2172  tmpCascade->GetSecondaryVtx()->RemoveDaughters();
2173  tmpCascade->UnsetOwnPrimaryVtx();
2174  delete tmpCascade; tmpCascade=NULL;
2176  rd2Prong->UnsetOwnPrimaryVtx();
2177  }
2178  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2179  //---
2180 
2181 
2182  return theCascade;
2183 }
2184 
2185 
2186 //----------------------------------------------------------------------------
2188  TObjArray *twoTrackArray,AliVEvent *event,
2189  AliAODVertex *secVert,
2190  AliAODv0 *v0,
2191  Double_t dca,
2192  Bool_t &okCascades)
2193 {
2196  //AliCodeTimerAuto("",0);
2197  UInt_t ntref=TProcessID::GetObjectCount();
2198  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2199  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2200  }
2201 
2202  // preselection to reduce the number of TRefs
2203  Double_t px[2],py[2],pz[2];
2204  AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2205  AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
2206  // propagate tracks to secondary vertex, to compute inv. mass
2207  postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2208  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2209  Double_t momentum[3];
2210  postrack->GetPxPyPz(momentum);
2211  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2212  negtrack->GetPxPyPz(momentum);
2213  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2214  if(!SelectInvMassAndPtCascade(px,py,pz)) return 0x0;
2215  Double_t dummyd0[2]={0,0};
2216  Double_t dummyd0err[2]={0,0};
2217  AliAODRecoCascadeHF tmpCasc(0x0,postrack->Charge(),px, py, pz, dummyd0, dummyd0err,0.);
2218  Bool_t presel=fCutsLctoV0->PreSelect(&tmpCasc,v0,postrack);
2219  if(!presel) return 0x0;
2220 
2221  // AliDebug(2,Form(" building the cascade"));
2222  okCascades = kFALSE;
2223  Bool_t dummy1,dummy2,dummy3;
2224 
2225  // We use Make2Prong to construct the AliAODRecoCascadeHF
2226  // (which inherits from AliAODRecoDecayHF2Prong)
2227  AliAODRecoCascadeHF *theCascade =
2228  (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
2229  dummy1,dummy2,dummy3);
2230  if(!theCascade) return 0x0;
2231 
2232  // bachelor track and charge
2233  AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2234  theCascade->SetCharge(trackBachelor->Charge());
2235 
2236  //--- selection cuts
2237  //
2238 
2239  AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
2240  if(fInputAOD){
2241  Int_t idBachelor=(Int_t)trackBachelor->GetID();
2242  if (idBachelor > -1 && idBachelor < fAODMapSize) {
2243  AliAODTrack* trackBachelorAOD=dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[idBachelor]));
2244  if(!trackBachelorAOD) AliFatal("Not a standard AOD");
2245  tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelorAOD);
2246  }
2247  }else{
2248  tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
2249  }
2250  tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
2251 
2252  AliAODVertex *primVertexAOD=0;
2254  // take event primary vertex
2255  primVertexAOD = PrimaryVertex();
2256  if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
2257  tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
2258  }
2259 
2260  // select Cascades
2261  if (fCascades && fInputAOD) {
2262  if (fCutsLctoV0->IsSelected(tmpCascade, AliRDHFCuts::kCandidate)>0) {
2263  okCascades = kTRUE;
2265  }
2266  if (fCutsDplustoK0spi) {
2268  okCascades = kTRUE;
2270  }
2271  }
2272  if (fCutsDstoK0sK){
2273  if (fCutsDstoK0sK->IsSelected(tmpCascade, AliRDHFCuts::kCandidate)>0) {
2274  okCascades = kTRUE;
2276  }
2277  }
2278  }
2279  else {
2280  //AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied"));
2281  okCascades=kTRUE;
2282  } // no cuts implemented from ESDs
2283  tmpCascade->GetSecondaryVtx()->RemoveDaughters();
2284  tmpCascade->UnsetOwnPrimaryVtx();
2285  delete tmpCascade; tmpCascade=NULL;
2286  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2287  //---
2288 
2289  return theCascade;
2290 }
2291 
2292 //-----------------------------------------------------------------------------
2294  TObjArray *twoTrackArray,AliVEvent *event,
2295  AliAODVertex *secVert,Double_t dca,
2296  Bool_t &okD0,Bool_t &okJPSI,
2297  Bool_t &okD0fromDstar, Bool_t refill, AliAODRecoDecayHF2Prong *rd)
2298 {
2301  // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
2302  // AliCodeTimerAuto("",0);
2303  okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
2304 
2305  UInt_t ntref=TProcessID::GetObjectCount();
2306  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2307  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2308  }
2309 
2310  Double_t px[2],py[2],pz[2],d0[2],d0err[2];
2311  AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2312  AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
2313 
2314  // propagate tracks to secondary vertex, to compute inv. mass
2315  postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2316  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2317 
2318  Double_t momentum[3];
2319  postrack->GetPxPyPz(momentum);
2320  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2321  negtrack->GetPxPyPz(momentum);
2322  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2323 
2324  if(!refill){//skip if it is called in refill step because already checked
2325  // invariant mass cut (try to improve coding here..)
2326  Bool_t okMassCut=kFALSE;
2327  if(!okMassCut && fD0toKpi) if(SelectInvMassAndPtD0Kpi(px,py,pz)) okMassCut=kTRUE;
2328  if(!okMassCut && fJPSItoEle) if(SelectInvMassAndPtJpsiee(px,py,pz)) okMassCut=kTRUE;
2329  if(!okMassCut && fDstar) if(SelectInvMassAndPtDstarD0pi(px,py,pz)) okMassCut=kTRUE;
2330  if(!okMassCut && fCascades) if(SelectInvMassAndPtCascade(px,py,pz)) okMassCut=kTRUE;
2331  if(!okMassCut) {
2332  //AliDebug(2," candidate didn't pass mass cut");
2333  return 0x0;
2334  }
2335  }
2336  // primary vertex to be used by this candidate
2337  AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
2338  if(!primVertexAOD) return 0x0;
2339 
2340  Double_t d0z0[2],covd0z0[3];
2341  postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2342  d0[0] = d0z0[0];
2343  d0err[0] = TMath::Sqrt(covd0z0[0]);
2344  negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2345  d0[1] = d0z0[0];
2346  d0err[1] = TMath::Sqrt(covd0z0[0]);
2347  AliAODRecoDecayHF2Prong *the2Prong;
2348  // create the object AliAODRecoDecayHF2Prong
2349  if(!refill){
2350  the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
2351  the2Prong->SetOwnPrimaryVtx(primVertexAOD);
2352  UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
2353  the2Prong->SetProngIDs(2,id);
2354  if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
2355  // Add daughter references already here
2356  if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
2357 
2358  // select D0->Kpi
2359  if(fD0toKpi) {
2361  if(okD0) the2Prong->SetSelectionBit(AliRDHFCuts::kD0toKpiCuts);
2362  }
2363  //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
2364  // select J/psi from B
2365  if(fJPSItoEle) {
2367  }
2368  //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
2369  // select D0->Kpi from Dstar
2370  if(fDstar) {
2371  okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
2372  if(okD0fromDstar) the2Prong->SetSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
2373  }
2374  //if(fDebug && fDstar) printf(" %d\n",(Int_t)okD0fromDstar);
2375  }
2376  }else{
2377  the2Prong =rd;
2378  the2Prong->SetSecondaryVtx(secVert);
2379  secVert->SetParent(the2Prong);
2380  AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
2381  the2Prong->SetOwnPrimaryVtx(primVertexAOD);
2382  the2Prong->SetPxPyPzProngs(2,px,py,pz);
2383  the2Prong->SetDCA(dca);
2384  the2Prong->Setd0Prongs(2,d0);
2385  the2Prong->Setd0errProngs(2,d0err);
2386  the2Prong->SetCharge(0);
2387  }
2388  delete primVertexAOD; primVertexAOD=NULL;
2389 
2390  // remove the primary vertex (was used only for selection)
2392  the2Prong->UnsetOwnPrimaryVtx();
2393  }
2394 
2395  // get PID info from ESD
2396  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2397  if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
2398  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2399  if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2400  Double_t esdpid[10];
2401  for(Int_t i=0;i<5;i++) {
2402  esdpid[i] = esdpid0[i];
2403  esdpid[5+i] = esdpid1[i];
2404  }
2405  the2Prong->SetPID(2,esdpid);
2406  return the2Prong;
2407 }
2408 //----------------------------------------------------------------------------
2410  TObjArray *threeTrackArray,AliVEvent *event,
2411  AliAODVertex *secVert,Double_t dispersion,
2412  const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
2413  Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
2414  Bool_t useForLc, Bool_t useForDs, Bool_t &ok3Prong)
2415 {
2418  // E.Bruna, F.Prino
2419  // AliCodeTimerAuto("",0);
2420 
2421 
2422  UInt_t ntref=TProcessID::GetObjectCount();
2423  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2424  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2425  }
2426 
2427  ok3Prong=kFALSE;
2428  if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
2429 
2430  Double_t px[3],py[3],pz[3],d0[3],d0err[3];
2431  Double_t momentum[3];
2432 
2433 
2434  AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
2435  AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
2436  AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
2437 
2438  postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2439  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2440  postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2441  postrack1->GetPxPyPz(momentum);
2442  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2443  negtrack->GetPxPyPz(momentum);
2444  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2445  postrack2->GetPxPyPz(momentum);
2446  px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
2447 
2448  // invariant mass cut for D+, Ds, Lc
2449  Bool_t okMassCut=kFALSE;
2450  if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
2451  if(!okMassCut && f3Prong) if(SelectInvMassAndPt3prong(px,py,pz)) okMassCut=kTRUE;
2452  if(!okMassCut) {
2453  //AliDebug(2," candidate didn't pass mass cut");
2454  return 0x0;
2455  }
2456 
2457  // primary vertex to be used by this candidate
2458  AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
2459  if(!primVertexAOD) return 0x0;
2460 
2461  Double_t d0z0[2],covd0z0[3];
2462  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2463  d0[0]=d0z0[0];
2464  d0err[0] = TMath::Sqrt(covd0z0[0]);
2465  negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2466  d0[1]=d0z0[0];
2467  d0err[1] = TMath::Sqrt(covd0z0[0]);
2468  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2469  d0[2]=d0z0[0];
2470  d0err[2] = TMath::Sqrt(covd0z0[0]);
2471 
2472 
2473  // create the object AliAODRecoDecayHF3Prong
2474  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2475  Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
2476  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]));
2477  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]));
2478  Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
2479 
2480 
2481  // construct the candidate passing a NULL pointer for the secondary vertex to avoid creation of TRef
2482  AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(0x0,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
2483  the3Prong->SetOwnPrimaryVtx(primVertexAOD);
2484  // add a pointer to the secondary vertex via SetOwnSecondaryVtx (no TRef created)
2485  AliAODVertex* ownsecv=secVert->CloneWithoutRefs();
2486  the3Prong->SetOwnSecondaryVtx(ownsecv);
2487  UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
2488  the3Prong->SetProngIDs(3,id);
2489 
2490  delete primVertexAOD; primVertexAOD=NULL;
2491 
2492  // disable PID, which requires the TRefs to the daughter tracks
2493  fCutsDplustoKpipi->SetUsePID(kFALSE);
2494  fCutsDstoKKpi->SetUsePID(kFALSE);
2495  fCutsLctopKpi->SetUsePID(kFALSE);
2496 
2497  // select D+->Kpipi, Ds->KKpi, Lc->pKpi
2498  if(f3Prong) {
2499  ok3Prong = kFALSE;
2500 
2502  ok3Prong = kTRUE;
2504  }
2505  if(useForDs && fOKInvMassDs){
2507  ok3Prong = kTRUE;
2509  }
2510  }
2511  if(useForLc && fOKInvMassLc){
2513  ok3Prong = kTRUE;
2515  }
2516  }
2517  }
2518  //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
2519  the3Prong->UnsetOwnSecondaryVtx();
2521  the3Prong->UnsetOwnPrimaryVtx();
2522  }
2523 
2524  // Add TRefs to secondary vertex and daughter tracks only for candidates passing the filtering cuts
2525  if(ok3Prong && fInputAOD){
2526  the3Prong->SetSecondaryVtx(secVert);
2527  AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
2528  }
2529 
2530  // get PID info from ESD
2531  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2532  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2533  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2534  if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2535  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2536  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2537 
2538  Double_t esdpid[15];
2539  for(Int_t i=0;i<5;i++) {
2540  esdpid[i] = esdpid0[i];
2541  esdpid[5+i] = esdpid1[i];
2542  esdpid[10+i] = esdpid2[i];
2543  }
2544  the3Prong->SetPID(3,esdpid);
2545 
2546  return the3Prong;
2547 }
2548 //----------------------------------------------------------------------------
2550  TObjArray *threeTrackArray,AliVEvent *event,
2551  AliAODVertex *secVert,Double_t dispersion,
2552  Double32_t dist12, Double32_t dist23,
2553  Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
2555 {
2556  // Fill on-the-fly the 3prong data member missing info
2557  // set TRef of vertex and daughters
2558  // do not recalculate the two-track secondary vertex dist12 and dist23
2559  // because they are stored in the AliAODRecoDecayHF3Prong candidate
2560  // reconstructed in the FindCandidate step
2561  // do not check if it is a Lambdac, Ds or Dplus because it is already check in the FindCandidate step
2562  // and the info is stored
2563  // AliCodeTimerAuto("",0);
2564 
2565  UInt_t ntref=TProcessID::GetObjectCount();
2566  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2567  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2568  }
2569 
2570  Double_t px[3],py[3],pz[3],d0[3],d0err[3];
2571  Double_t momentum[3];
2572 
2573 
2574  AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
2575  AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
2576  AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
2577 
2578  postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2579  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2580  postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2581  postrack1->GetPxPyPz(momentum);
2582  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2583  negtrack->GetPxPyPz(momentum);
2584  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2585  postrack2->GetPxPyPz(momentum);
2586  px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
2587  // primary vertex to be used by this candidate
2588  AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
2589  if(!primVertexAOD) return 0x0;
2590  Double_t d0z0[2],covd0z0[3];
2591  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2592  d0[0]=d0z0[0];
2593  d0err[0] = TMath::Sqrt(covd0z0[0]);
2594  negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2595  d0[1]=d0z0[0];
2596  d0err[1] = TMath::Sqrt(covd0z0[0]);
2597  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2598  d0[2]=d0z0[0];
2599  d0err[2] = TMath::Sqrt(covd0z0[0]);
2600 
2601  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2602  Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
2603  Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
2604 
2605  rd->SetSecondaryVtx(secVert);
2606  secVert->SetParent(rd);
2607  AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
2608 
2609  rd->SetPxPyPzProngs(3,px,py,pz);
2610  rd->SetDCAs(3,dca);
2611  rd->Setd0Prongs(3,d0);
2612  rd->Setd0errProngs(3,d0err);
2613  rd->SetCharge(charge);
2614  rd->SetOwnPrimaryVtx(primVertexAOD);
2615  rd->SetSigmaVert(dispersion);
2616  delete primVertexAOD; primVertexAOD=NULL;
2617 
2619  rd->UnsetOwnPrimaryVtx();
2620  }
2621 
2622  // get PID info from ESD
2623  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2624  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2625  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2626  if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2627  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2628  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2629 
2630  Double_t esdpid[15];
2631  for(Int_t i=0;i<5;i++) {
2632  esdpid[i] = esdpid0[i];
2633  esdpid[5+i] = esdpid1[i];
2634  esdpid[10+i] = esdpid2[i];
2635  }
2636  rd->SetPID(3,esdpid);
2637  return rd;
2638 }
2639 //----------------------------------------------------------------------------
2641  TObjArray *fourTrackArray,AliVEvent *event,
2642  AliAODVertex *secVert,
2643  const AliAODVertex *vertexp1n1,
2644  const AliAODVertex *vertexp1n1p2,
2645  Double_t dcap1n1,Double_t dcap1n2,
2646  Double_t dcap2n1,Double_t dcap2n2,
2647  Bool_t &ok4Prong)
2648 {
2651  // G.E.Bruno, R.Romita
2652  // AliCodeTimerAuto("",0);
2653 
2654  UInt_t ntref=TProcessID::GetObjectCount();
2655  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2656  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2657  }
2658 
2659  ok4Prong=kFALSE;
2660  if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
2661 
2662  Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
2663 
2664  AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
2665  AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
2666  AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
2667  AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
2668 
2669  postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2670  negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2671  postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2672  negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2673 
2674  Double_t momentum[3];
2675  postrack1->GetPxPyPz(momentum);
2676  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2677  negtrack1->GetPxPyPz(momentum);
2678  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2679  postrack2->GetPxPyPz(momentum);
2680  px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
2681  negtrack2->GetPxPyPz(momentum);
2682  px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
2683 
2684  // invariant mass cut for rho or D0 (try to improve coding here..)
2685  Bool_t okMassCut=kFALSE;
2686  if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
2687  if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) { //no PID, to be implemented with PID
2688  if(SelectInvMassAndPt4prong(px,py,pz)) okMassCut=kTRUE;
2689  }
2690  if(!okMassCut) {
2691  //if(fDebug) printf(" candidate didn't pass mass cut\n");
2692  //printf(" candidate didn't pass mass cut\n");
2693  return 0x0;
2694  }
2695 
2696  // primary vertex to be used by this candidate
2697  AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
2698  if(!primVertexAOD) return 0x0;
2699 
2700  Double_t d0z0[2],covd0z0[3];
2701  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2702  d0[0]=d0z0[0];
2703  d0err[0] = TMath::Sqrt(covd0z0[0]);
2704  negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2705  d0[1]=d0z0[0];
2706  d0err[1] = TMath::Sqrt(covd0z0[0]);
2707  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2708  d0[2]=d0z0[0];
2709  d0err[2] = TMath::Sqrt(covd0z0[0]);
2710  negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2711  d0[3]=d0z0[0];
2712  d0err[3] = TMath::Sqrt(covd0z0[0]);
2713 
2714 
2715  // create the object AliAODRecoDecayHF4Prong
2716  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2717  Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
2718  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]));
2719  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]));
2720  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]));
2721  Short_t charge=0;
2722  AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
2723  the4Prong->SetOwnPrimaryVtx(primVertexAOD);
2724  UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
2725  the4Prong->SetProngIDs(4,id);
2726 
2727  delete primVertexAOD; primVertexAOD=NULL;
2728 
2730 
2731 
2733  the4Prong->UnsetOwnPrimaryVtx();
2734  }
2735 
2736 
2737  // get PID info from ESD
2738  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2739  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2740  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2741  if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
2742  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2743  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2744  Double_t esdpid3[5]={0.,0.,0.,0.,0.};
2745  if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
2746 
2747  Double_t esdpid[20];
2748  for(Int_t i=0;i<5;i++) {
2749  esdpid[i] = esdpid0[i];
2750  esdpid[5+i] = esdpid1[i];
2751  esdpid[10+i] = esdpid2[i];
2752  esdpid[15+i] = esdpid3[i];
2753  }
2754  the4Prong->SetPID(4,esdpid);
2755 
2756  return the4Prong;
2757 }
2758 //----------------------------------------------------------------------------------
2760  //assign and save in fAODMap the index of the AliAODTrack track
2761  //ordering them on the basis of selected criteria
2762 
2763  fAODMapSize = 100000;
2764  fAODMap = new Int_t[fAODMapSize];
2765  AliAODTrack *track=0;
2766  memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
2767  for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
2768  track = dynamic_cast<AliAODTrack*>(aod->GetTrack(i));
2769  if(!track) AliFatal("Not a standard AOD");
2770  // skip pure ITS SA tracks
2771  if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2772 
2773  // skip tracks without ITS
2774  if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
2775 
2776  // TEMPORARY: check that the cov matrix is there
2777  Double_t covtest[21];
2778  if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2779  //
2780 
2781  Int_t ind = (Int_t)track->GetID();
2782  if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
2783  }
2784  return;
2785 }
2786 //-----------------------------------------------------------------------------
2787 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
2788  AliVEvent *event) const
2789 {
2791  //AliCodeTimerAuto("",0);
2792 
2793  AliESDVertex *vertexESD = 0;
2794  AliAODVertex *vertexAOD = 0;
2795 
2796 
2798  // primary vertex from the input event
2799 
2800  vertexESD = new AliESDVertex(*fV1);
2801 
2802  } else {
2803  // primary vertex specific to this candidate
2804 
2805  Int_t nTrks = trkArray->GetEntriesFast();
2806  AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
2807 
2809  // recalculating the vertex
2810 
2811  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
2812  Float_t diamondcovxy[3];
2813  event->GetDiamondCovXY(diamondcovxy);
2814  Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
2815  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
2816  AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
2817  vertexer->SetVtxStart(diamond);
2818  delete diamond; diamond=NULL;
2819  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
2820  vertexer->SetOnlyFitter();
2821  }
2822  Int_t skipped[1000];
2823  Int_t nTrksToSkip=0,id;
2824  AliExternalTrackParam *t = 0;
2825  for(Int_t i=0; i<nTrks; i++) {
2826  t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
2827  id = (Int_t)t->GetID();
2828  if(id<0) continue;
2829  skipped[nTrksToSkip++] = id;
2830  }
2831  // TEMPORARY FIX
2832  // For AOD, skip also tracks without covariance matrix
2833  if(fInputAOD) {
2834  Double_t covtest[21];
2835  for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
2836  AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
2837  if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
2838  id = (Int_t)vtrack->GetID();
2839  if(id<0) continue;
2840  skipped[nTrksToSkip++] = id;
2841  }
2842  }
2843  }
2844  for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
2845  //
2846  vertexer->SetSkipTracks(nTrksToSkip,skipped);
2847  vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
2848 
2849  } else if(fRmTrksFromPrimVtx && nTrks>0) {
2850  // removing the prongs tracks
2851 
2852  TObjArray rmArray(nTrks);
2853  UShort_t *rmId = new UShort_t[nTrks];
2854  AliESDtrack *esdTrack = 0;
2855  AliESDtrack *t = 0;
2856  for(Int_t i=0; i<nTrks; i++) {
2857  t = (AliESDtrack*)trkArray->UncheckedAt(i);
2858  esdTrack = new AliESDtrack(*t);
2859  rmArray.AddLast(esdTrack);
2860  if(esdTrack->GetID()>=0) {
2861  rmId[i]=(UShort_t)esdTrack->GetID();
2862  } else {
2863  rmId[i]=9999;
2864  }
2865  }
2866  Float_t diamondxy[2]={static_cast<Float_t>(event->GetDiamondX()),static_cast<Float_t>(event->GetDiamondY())};
2867  vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
2868  delete [] rmId; rmId=NULL;
2869  rmArray.Delete();
2870 
2871  }
2872 
2873  if(!vertexESD) return vertexAOD;
2874  if(vertexESD->GetNContributors()<=0) {
2875  //AliDebug(2,"vertexing failed");
2876  delete vertexESD; vertexESD=NULL;
2877  return vertexAOD;
2878  }
2879 
2880  delete vertexer; vertexer=NULL;
2881 
2882  }
2883 
2884  // convert to AliAODVertex
2885  Double_t pos[3],cov[6],chi2perNDF;
2886  vertexESD->GetXYZ(pos); // position
2887  vertexESD->GetCovMatrix(cov); //covariance matrix
2888  chi2perNDF = vertexESD->GetChi2toNDF();
2889  delete vertexESD; vertexESD=NULL;
2890 
2891  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
2892 
2893  return vertexAOD;
2894 }
2895 //-----------------------------------------------------------------------------
2898 
2899  //printf("Preselections:\n");
2900  // fTrackFilter->Dump();
2901  if(fSecVtxWithKF) {
2902  printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
2903  } else {
2904  printf("Secondary vertex with AliVertexerTracks\n");
2905  }
2906  if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
2907  if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
2908  if(fD0toKpi) {
2909  printf("Reconstruct D0->Kpi candidates with cuts:\n");
2911  }
2912  if(fDstar) {
2913  printf("Reconstruct D*->D0pi candidates with cuts:\n");
2914  if(fFindVertexForDstar) {
2915  printf(" Reconstruct a secondary vertex for the D*\n");
2916  } else {
2917  printf(" Assume the D* comes from the primary vertex\n");
2918  }
2920  }
2921  if(fJPSItoEle) {
2922  printf("Reconstruct J/psi from B candidates with cuts:\n");
2924  }
2925  if(f3Prong) {
2926  printf("Reconstruct 3 prong candidates.\n");
2927  printf(" D+->Kpipi cuts:\n");
2929  printf(" Ds->KKpi cuts:\n");
2931  printf(" Lc->pKpi cuts:\n");
2933  }
2934  if(f4Prong) {
2935  printf("Reconstruct 4 prong candidates.\n");
2936  printf(" D0->Kpipipi cuts:\n");
2938  }
2939  if(fCascades) {
2940  printf("Reconstruct cascade candidates formed with v0s.\n");
2941  printf(" Lc -> k0s P & Lc -> L Pi cuts:\n");
2943  printf(" D+ -> K0s pi cuts:\n");
2945  printf(" Ds -> K0s K cuts:\n");
2947  }
2948 
2949  return;
2950 }
2951 //-----------------------------------------------------------------------------
2953  Double_t &dispersion,Bool_t useTRefArray) const
2954 {
2956  //AliCodeTimerAuto("",0);
2957 
2958  AliESDVertex *vertexESD = 0;
2959  AliAODVertex *vertexAOD = 0;
2960 
2961  if(!fSecVtxWithKF) { // AliVertexerTracks
2962 
2963  fVertexerTracks->SetVtxStart(fV1);
2964  vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
2965 
2966  if(!vertexESD) return vertexAOD;
2967 
2968  if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
2969  //AliDebug(2,"vertexing failed");
2970  delete vertexESD; vertexESD=NULL;
2971  return vertexAOD;
2972  }
2973 
2974  Double_t vertRadius2=vertexESD->GetX()*vertexESD->GetX()+vertexESD->GetY()*vertexESD->GetY();
2975  if(vertRadius2>8.){
2976  // vertex outside beam pipe, reject candidate to avoid propagation through material
2977  delete vertexESD; vertexESD=NULL;
2978  return vertexAOD;
2979  }
2980 
2981  } else { // Kalman Filter vertexer (AliKFParticle)
2982 
2983  AliKFParticle::SetField(fBzkG);
2984 
2985  AliKFVertex vertexKF;
2986 
2987  Int_t nTrks = trkArray->GetEntriesFast();
2988  for(Int_t i=0; i<nTrks; i++) {
2989  AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
2990  AliKFParticle daughterKF(*esdTrack,211);
2991  vertexKF.AddDaughter(daughterKF);
2992  }
2993  vertexESD = new AliESDVertex(vertexKF.Parameters(),
2994  vertexKF.CovarianceMatrix(),
2995  vertexKF.GetChi2(),
2996  vertexKF.GetNContributors());
2997 
2998  }
2999 
3000  // convert to AliAODVertex
3001  Double_t pos[3],cov[6],chi2perNDF;
3002  vertexESD->GetXYZ(pos); // position
3003  vertexESD->GetCovMatrix(cov); //covariance matrix
3004  chi2perNDF = vertexESD->GetChi2toNDF();
3005  dispersion = vertexESD->GetDispersion();
3006  delete vertexESD; vertexESD=NULL;
3007 
3008  Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
3009  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
3010 
3011  return vertexAOD;
3012 }
3013 //-----------------------------------------------------------------------------
3016  //AliCodeTimerAuto("",0);
3017 
3018  Int_t retval=kFALSE;
3019  Double_t momentum[3];
3020  Double_t px[3],py[3],pz[3];
3021  for(Int_t iTrack=0; iTrack<3; iTrack++){
3022  AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
3023  track->GetPxPyPz(momentum);
3024  px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
3025  }
3026  retval = SelectInvMassAndPt3prong(px,py,pz);
3027 
3028  return retval;
3029 }
3030 
3031 //-----------------------------------------------------------------------------
3034  //AliCodeTimerAuto("",0);
3035 
3036  Int_t retval=kFALSE;
3037  Double_t momentum[3];
3038  Double_t px[4],py[4],pz[4];
3039 
3040  for(Int_t iTrack=0; iTrack<4; iTrack++){
3041  AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
3042  track->GetPxPyPz(momentum);
3043  px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
3044  }
3045 
3046  retval = SelectInvMassAndPt4prong(px,py,pz);
3047 
3048  return retval;
3049 }
3050 //-----------------------------------------------------------------------------
3053  //AliCodeTimerAuto("",0);
3054 
3055  Int_t retval=kFALSE;
3056  Double_t momentum[3];
3057  Double_t px[2],py[2],pz[2];
3058 
3059  for(Int_t iTrack=0; iTrack<2; iTrack++){
3060  AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
3061  track->GetPxPyPz(momentum);
3062  px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
3063  }
3064  retval = SelectInvMassAndPtDstarD0pi(px,py,pz);
3065 
3066  return retval;
3067 }
3068 //-----------------------------------------------------------------------------
3070  Double_t *py,
3071  Double_t *pz){
3073  //AliCodeTimerAuto("",0);
3074 
3075  UInt_t pdg2[2];
3076  Int_t nprongs=2;
3077  Double_t minv2,mrange;
3078  Double_t lolim,hilim;
3079  Double_t minPt=0;
3080  Bool_t retval=kFALSE;
3081 
3082  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3083  fOKInvMassD0=kFALSE;
3084  // pt cut
3086  if(minPt>0.1)
3087  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3088  // mass cut
3089  mrange=fCutsD0toKpi->GetMassCut();
3090  lolim=fMassDzero-mrange;
3091  hilim=fMassDzero+mrange;
3092  pdg2[0]=211; pdg2[1]=321;
3093  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3094  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3095  retval=kTRUE;
3096  fOKInvMassD0=kTRUE;
3097  }
3098  pdg2[0]=321; pdg2[1]=211;
3099  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3100  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3101  retval=kTRUE;
3102  fOKInvMassD0=kTRUE;
3103  }
3104  return retval;
3105 }
3106 
3107 //-----------------------------------------------------------------------------
3109  Double_t *py,
3110  Double_t *pz){
3112  //AliCodeTimerAuto("",0);
3113 
3114  UInt_t pdg2[2];
3115  Int_t nprongs=2;
3116  Double_t minv2,mrange;
3117  Double_t lolim,hilim;
3118  Double_t minPt=0;
3119  Bool_t retval=kFALSE;
3120 
3121  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3122  fOKInvMassJpsi=kFALSE;
3123  // pt cut
3125  if(minPt>0.1)
3126  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3127  // mass cut
3128  mrange=fCutsJpsitoee->GetMassCut();
3129  lolim=fMassJpsi-mrange;
3130  hilim=fMassJpsi+mrange;
3131 
3132  pdg2[0]=11; pdg2[1]=11;
3133  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3134  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3135  retval=kTRUE;
3136  fOKInvMassJpsi=kTRUE;
3137  }
3138 
3139  return retval;
3140 }
3141 //-----------------------------------------------------------------------------
3143  Double_t *py,
3144  Double_t *pz,
3145  Int_t pidLcStatus){
3147  //AliCodeTimerAuto("",0);
3148 
3149  UInt_t pdg3[3];
3150  Int_t nprongs=3;
3151  Double_t minv2,mrange;
3152  Double_t lolim,hilim;
3153  Double_t minPt=0;
3154  Bool_t retval=kFALSE;
3155 
3156 
3157  fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
3158  fOKInvMassDplus=kFALSE;
3159  fOKInvMassDs=kFALSE;
3160  fOKInvMassLc=kFALSE;
3161  // pt cut
3163  minPt=TMath::Min(minPt,fCutsLctopKpi->GetMinPtCandidate());
3164  if(minPt>0.1)
3165  if(fMassCalc3->Pt2() < minPt*minPt) return retval;
3166  // D+->Kpipi
3167  mrange=fCutsDplustoKpipi->GetMassCut();
3168  lolim=fMassDplus-mrange;
3169  hilim=fMassDplus+mrange;
3170  pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
3171  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3172  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3173  retval=kTRUE;
3174  fOKInvMassDplus=kTRUE;
3175  }
3176  // Ds+->KKpi
3177  mrange=fCutsDstoKKpi->GetMassCut();
3178  lolim=fMassDs-mrange;
3179  hilim=fMassDs+mrange;
3180  pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
3181  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3182  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3183  retval=kTRUE;
3184  fOKInvMassDs=kTRUE;
3185  }
3186  pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
3187  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3188  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3189  retval=kTRUE;
3190  fOKInvMassDs=kTRUE;
3191  }
3192  // Lc->pKpi
3193  mrange=fCutsLctopKpi->GetMassCut();
3194  lolim=fMassLambdaC-mrange;
3195  hilim=fMassLambdaC+mrange;
3196  if(pidLcStatus&1){
3197  pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
3198  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3199  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3200  retval=kTRUE;
3201  fOKInvMassLc=kTRUE;
3202  }
3203  }
3204  if(pidLcStatus&2){
3205  pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
3206  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3207  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3208  retval=kTRUE;
3209  fOKInvMassLc=kTRUE;
3210  }
3211  }
3212 
3213  return retval;
3214 }
3215 
3216 //-----------------------------------------------------------------------------
3218  Double_t *py,
3219  Double_t *pz){
3221  //AliCodeTimerAuto("",0);
3222 
3223  UInt_t pdg2[2];
3224  Int_t nprongs=2;
3225  Double_t minv2,mrange;
3226  Double_t lolim,hilim;
3227  Double_t minPt=0;
3228  Bool_t retval=kFALSE;
3229 
3230  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3231  fOKInvMassDstar=kFALSE;
3232  // pt cut
3234  if(minPt>0.1)
3235  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3236  // mass cut
3237  pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
3238  mrange=fCutsDStartoKpipi->GetMassCut();
3239  lolim=fMassDstar-mrange;
3240  hilim=fMassDstar+mrange;
3241  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3242  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3243  retval=kTRUE;
3244  fOKInvMassDstar=kTRUE;
3245  }
3246 
3247  return retval;
3248 }
3249 
3250 //-----------------------------------------------------------------------------
3252  Double_t *py,
3253  Double_t *pz){
3255  //AliCodeTimerAuto("",0);
3256 
3257  UInt_t pdg4[4];
3258  Int_t nprongs=4;
3259  Double_t minv2,mrange;
3260  Double_t lolim,hilim;
3261  Double_t minPt=0;
3262  Bool_t retval=kFALSE;
3263 
3264  // D0->Kpipipi without PID
3265  fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
3266  fOKInvMassD0to4p=kFALSE;
3267  // pt cut
3269  if(minPt>0.1)
3270  if(fMassCalc4->Pt2() < minPt*minPt) return retval;
3271  // mass cut
3272  mrange=fCutsD0toKpipipi->GetMassCut();
3273  lolim=fMassDzero-mrange;
3274  hilim=fMassDzero+mrange;
3275 
3276  pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
3277  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3278  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3279  retval=kTRUE;
3280  fOKInvMassD0to4p=kTRUE;
3281  }
3282 
3283  pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
3284  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3285  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3286  retval=kTRUE;
3287  fOKInvMassD0to4p=kTRUE;
3288  }
3289 
3290  pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
3291  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3292  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3293  retval=kTRUE;
3294  fOKInvMassD0to4p=kTRUE;
3295  }
3296 
3297  pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
3298  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3299  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3300  retval=kTRUE;
3301  fOKInvMassD0to4p=kTRUE;
3302  }
3303 
3304  return retval;
3305 }
3306 //-----------------------------------------------------------------------------
3308  Double_t *py,
3309  Double_t *pz){
3311  //AliCodeTimerAuto("",0);
3312 
3313  UInt_t pdg2[2];
3314  Int_t nprongs=2;
3315  Double_t minv2,mrange;
3316  Double_t lolim,hilim;
3317  Double_t minPt=0;
3318  Bool_t retval=kFALSE;
3319 
3320  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3321  minPt=fCutsLctoV0->GetMinPtCandidate();
3322  if(minPt>0.1)
3323  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3324 
3325  fOKInvMassLctoV0 = kFALSE;
3326 
3327  // LambdaC candidate
3328  if (fCutsLctoV0) {
3329  mrange = fCutsLctoV0->GetMassCut();
3330  lolim = fMassLambdaC - mrange;
3331  hilim = fMassLambdaC + mrange;
3332  pdg2[0] = 2212; pdg2[1] = 310;
3333  minv2 = fMassCalc2->InvMass2(2, pdg2);
3334  if ((minv2 > lolim*lolim) && (minv2 < hilim*hilim)) {
3335  retval = kTRUE;
3336  fOKInvMassLctoV0 = kTRUE;
3337  }
3338  pdg2[0] = 211; pdg2[1] = 3122;
3339  minv2 = fMassCalc2->InvMass2(2, pdg2);
3340  if ((minv2>lolim*lolim) && (minv2<hilim*hilim)) {
3341  retval = kTRUE;
3342  fOKInvMassLctoV0 = kTRUE;
3343  }
3344  }
3345 
3346  // D+ cascade candidate
3347  if (fCutsDplustoK0spi) {
3348  mrange = fCutsDplustoK0spi->GetMassCut();
3349  lolim = fMassDplus - mrange;
3350  hilim = fMassDplus + mrange;
3351  pdg2[0] = 211; pdg2[1] = 310;
3352  minv2 = fMassCalc2->InvMass2(2, pdg2);
3353  if ((minv2 > lolim*lolim) && (minv2 < hilim*hilim)) {
3354  retval = kTRUE;
3355  }
3356  }
3357 
3358  // Ds cascade candidate
3359  if (fCutsDstoK0sK) {
3360  mrange = fCutsDstoK0sK->GetMassCut();
3361  lolim = fMassDs - mrange;
3362  hilim = fMassDs + mrange;
3363  pdg2[0] = 321; pdg2[1] = 310;
3364  minv2 = fMassCalc2->InvMass2(2, pdg2);
3365  if ((minv2 > lolim*lolim) && (minv2 < hilim*hilim)) {
3366  retval = kTRUE;
3367  }
3368  }
3369 
3370  return retval;
3371 }
3372 //-----------------------------------------------------------------------------
3374  Int_t trkEntries,
3375  TObjArray &seleTrksArray,
3376  TObjArray &tracksAtVertex,
3377  Int_t &nSeleTrks,
3378  UChar_t *seleFlags,Int_t *evtNumber)
3379 {
3385  //AliCodeTimerAuto("",0);
3386 
3387  const AliVVertex *vprimary = event->GetPrimaryVertex();
3388 
3389  if(fV1) { delete fV1; fV1=NULL; }
3390  if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
3391 
3392  Int_t nindices=0;
3393  UShort_t *indices = 0;
3394  Double_t pos[3],cov[6];
3395  const Int_t entries = event->GetNumberOfTracks();
3396  AliCentrality* cent;
3397 
3398  if(!fInputAOD) { // ESD
3399  fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
3400  cent=((AliESDEvent*)event)->GetCentrality();
3401  } else { // AOD
3402  vprimary->GetXYZ(pos);
3403  vprimary->GetCovarianceMatrix(cov);
3404  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
3405  if(entries<=0) return;
3406  indices = new UShort_t[entries];
3407  memset(indices,0,sizeof(UShort_t)*entries);
3408  fAODMapSize = 100000;
3409  fAODMap = new Int_t[fAODMapSize];
3410  memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
3411  cent=((AliAODEvent*)event)->GetCentrality();
3412  }
3413  Float_t centperc=0.1;
3414  if(event->GetRunNumber()<244824){
3415  centperc=cent->GetCentralityPercentile("V0M");
3416  }else{
3417  AliMultSelection *multSelection = (AliMultSelection * ) event->FindListObject("MultSelection");
3418  if(multSelection){
3419  centperc=multSelection->GetMultiplicityPercentile("V0M");
3420  Int_t qual = multSelection->GetEvSelCode();
3421  if(qual == 199 ) centperc=0.1; // use central cuts as default
3422  }
3423  }
3424  Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE,okFor3Prong=kFALSE,okBachelor=kFALSE;
3425  nSeleTrks=0;
3426 
3427  // transfer ITS tracks from event to arrays
3428  for(Int_t i=0; i<entries; i++) {
3429  AliVTrack *track;
3430  track = (AliVTrack*)event->GetTrack(i);
3431 
3432  // skip pure ITS SA tracks
3433  if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
3434 
3435  // skip tracks without ITS
3436  if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
3437 
3438  // skip tracks with negative ID
3439  // (these are duplicated TPC-only AOD tracks, for jet analysis...)
3440  if(track->GetID()<0) continue;
3441 
3442  // TEMPORARY: check that the cov matrix is there
3443  Double_t covtest[21];
3444  if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
3445  //
3446 
3447  if(fInputAOD) {
3448  AliAODTrack *aodt = (AliAODTrack*)track;
3449  if(aodt->GetUsedForPrimVtxFit()) {
3450  indices[nindices]=aodt->GetID(); nindices++;
3451  }
3452  Int_t ind = (Int_t)aodt->GetID();
3453  if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
3454  }
3455 
3456  AliESDtrack *esdt = 0;
3457 
3458  if(!fInputAOD) {
3459  esdt = (AliESDtrack*)track;
3460  } else {
3461  esdt = new AliESDtrack(track);
3462  }
3463 
3464  // single track selection
3465  okDisplaced=kFALSE; okSoftPi=kFALSE; okFor3Prong=kFALSE; okBachelor=kFALSE;
3466  if(fMixEvent && i<trkEntries){
3467  evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
3468  const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
3469  Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
3470  eventVtx->GetXYZ(vtxPos);
3471  vprimary->GetXYZ(primPos);
3472  eventVtx->GetCovarianceMatrix(primCov);
3473  for(Int_t ind=0;ind<3;ind++){
3474  trasl[ind]=vtxPos[ind]-primPos[ind];
3475  }
3476 
3477  Bool_t isTransl=esdt->Translate(trasl,primCov);
3478  if(!isTransl) {
3479  delete esdt;
3480  esdt = NULL;
3481  continue;
3482  }
3483  }
3484 
3485  if(SingleTrkCuts(esdt,centperc,okDisplaced,okSoftPi,okFor3Prong,okBachelor) && nSeleTrks<trkEntries) {
3486  esdt->PropagateToDCA(fV1,fBzkG,kVeryBig);
3487  seleTrksArray.AddLast(esdt);
3488  tracksAtVertex.AddLast(new AliExternalTrackParam(*esdt));
3489  seleFlags[nSeleTrks]=0;
3490  if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
3491  if(okFor3Prong) SETBIT(seleFlags[nSeleTrks],kBit3Prong);
3492  if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
3493  if(okBachelor) SETBIT(seleFlags[nSeleTrks],kBitBachelor);
3494  // Check the PID
3495  SETBIT(seleFlags[nSeleTrks],kBitPionCompat);
3496  SETBIT(seleFlags[nSeleTrks],kBitKaonCompat);
3497  SETBIT(seleFlags[nSeleTrks],kBitProtonCompat);
3498  Bool_t useTPC=kTRUE;
3499  if(fUseTOFPID){
3500  Double_t nsigmatofPi= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kPion);
3501  if(nsigmatofPi>-990. && (nsigmatofPi<-fnSigmaTOFPionLow || nsigmatofPi>fnSigmaTOFPionHi)){
3502  CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
3503  }
3504  Double_t nsigmatofK= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kKaon);
3505  if(nsigmatofK>-990. && (nsigmatofK<-fnSigmaTOFKaonLow || nsigmatofK>fnSigmaTOFKaonHi)){
3506  CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
3507  }
3508  Double_t nsigmatofP= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kProton);
3509  if(nsigmatofP>-990. && (nsigmatofP<-fnSigmaTOFProtonLow || nsigmatofP>fnSigmaTOFProtonHi)){
3510  CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
3511  }
3512  if(fUseTPCPIDOnlyIfNoTOF && nsigmatofPi>-990.) useTPC=kFALSE;
3513  }
3514  if(useTPC && fUseTPCPID && esdt->P()<fMaxMomForTPCPid){
3515  Double_t nsigmatpcPi= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kPion);
3516  if(nsigmatpcPi>-990. && (nsigmatpcPi<-fnSigmaTPCPionLow || nsigmatpcPi>fnSigmaTPCPionHi)){
3517  CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
3518  }
3519  Double_t nsigmatpcK= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kKaon);
3520  if(nsigmatpcK>-990. && (nsigmatpcK<-fnSigmaTPCKaonLow || nsigmatpcK>fnSigmaTPCKaonHi)){
3521  CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
3522  }
3523  Double_t nsigmatpcP= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kProton);
3524  if(nsigmatpcP>-990. && (nsigmatpcP<-fnSigmaTPCProtonLow || nsigmatpcP>fnSigmaTPCProtonHi)){
3525  CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
3526  }
3527  }
3528  nSeleTrks++;
3529  } else {
3530  if(fInputAOD) delete esdt;
3531  esdt = NULL;
3532  continue;
3533  }
3534 
3535  } // end loop on tracks
3536 
3537  // primary vertex from AOD
3538  if(fInputAOD) {
3539  delete fV1; fV1=NULL;
3540  vprimary->GetXYZ(pos);
3541  vprimary->GetCovarianceMatrix(cov);
3542  Double_t chi2toNDF = vprimary->GetChi2perNDF();
3543  Int_t ncontr=nindices;
3544  if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
3545  Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
3546  fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
3547  fV1->SetTitle(vprimary->GetTitle());
3548  fV1->SetIndices(nindices,indices);
3549  }
3550  if(indices) { delete [] indices; indices=NULL; }
3551 
3552 
3553  return;
3554 }
3555 //-----------------------------------------------------------------------------
3557  //
3559  //
3560  if(fUsePidTag && cuts->GetPidHF()) {
3561  Bool_t usepid=cuts->GetIsUsePID();
3562  cuts->SetUsePID(kTRUE);
3563  if(cuts->IsSelectedPID(rd))
3564  rd->SetSelectionBit(bit);
3565  cuts->SetUsePID(usepid);
3566  }
3567  return;
3568 }
3569 //-----------------------------------------------------------------------------
3571  Float_t centralityperc,
3572  Bool_t &okDisplaced,
3573  Bool_t &okSoftPi,
3574  Bool_t &okFor3Prong,
3575  Bool_t &okBachelor) const
3576 {
3578 
3579  // this is needed to store the impact parameters
3580  //AliCodeTimerAuto("",0);
3581 
3582  if (!trk->PropagateToDCA(fV1,fBzkG,kVeryBig)) return kFALSE;
3583 
3584  trk->RelateToVertex(fV1,fBzkG,kVeryBig);
3585 
3586  UInt_t selectInfo;
3587  //
3588  // Track selection, displaced tracks -- 2 prong
3589  selectInfo = 0;
3590  if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
3592  // central PbPb events, tighter cuts
3593  selectInfo = fTrackFilter2prongCentral->IsSelected(trk);
3594  }else{
3595  // standard cuts
3596  if(fTrackFilter) {
3597  selectInfo = fTrackFilter->IsSelected(trk);
3598  }
3599  }
3600  if(selectInfo) okDisplaced=kTRUE;
3601 
3602  // Track selection, displaced tracks -- 3 prong
3603  selectInfo = 0;
3604  if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
3606  // central PbPb events, tighter cuts
3607  selectInfo = fTrackFilter3prongCentral->IsSelected(trk);
3608  }else{
3609  // standard cuts
3610  if(fTrackFilter) {
3611  selectInfo = fTrackFilter->IsSelected(trk);
3612  }
3613  }
3614  if(selectInfo) okFor3Prong=kTRUE;
3615 
3616  // Track selection, soft pions
3617  selectInfo = 0;
3618  if(fDstar && fTrackFilterSoftPi) {
3619  selectInfo = fTrackFilterSoftPi->IsSelected(trk);
3620  }
3621  if(selectInfo) okSoftPi=kTRUE;
3622 
3623  // Track selection, bachelor
3624  selectInfo = 0;
3625  if(fCascades){
3627  selectInfo = fTrackFilterBachelor->IsSelected(trk);
3628  }else{
3629  if(okDisplaced) selectInfo=1;
3630  }
3631  }
3632  if(selectInfo) okBachelor=kTRUE;
3633 
3634  if(okDisplaced || okSoftPi || okFor3Prong || okBachelor) return kTRUE;
3635 
3636  return kFALSE;
3637 }
3638 
3639 
3640 //-----------------------------------------------------------------------------
3641 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
3647  //
3648  //AliCodeTimerAuto("",0);
3649  Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
3650  AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
3651 
3652  // create the v0 neutral track to compute the DCA to the primary vertex
3653  Double_t xyz[3], pxpypz[3];
3654  esdV0->XvYvZv(xyz);
3655  esdV0->PxPyPz(pxpypz);
3656  Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
3657  AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
3658  if(!trackesdV0) {
3659  delete vertexV0;
3660  return 0;
3661  }
3662  Double_t d0z0[2],covd0z0[3];
3663  AliAODVertex *primVertexAOD = PrimaryVertex();
3664  trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
3665  Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
3666  // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
3667  Double_t dcaV0DaughterToPrimVertex[2];
3668  AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
3669  AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
3670  if( !posV0track || !negV0track) {
3671  if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
3672  delete vertexV0;
3673  delete primVertexAOD;
3674  return 0;
3675  }
3676  posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
3677  // if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
3678  // else
3679  dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
3680  negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
3681  // if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
3682  // else
3683  dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
3684  Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
3685  Double_t pmom[3],nmom[3];
3686  esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
3687  esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
3688 
3689  AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
3690  aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
3691 
3692  delete trackesdV0;
3693  delete primVertexAOD;
3694 
3695  return aodV0;
3696 }
3697 //-----------------------------------------------------------------------------
3698 void AliAnalysisVertexingHF::SetParametersAtVertex(AliESDtrack* esdt, const AliExternalTrackParam* extpar) const{
3700  //AliCodeTimerAuto("",0);
3701 
3702  const Double_t *par=extpar->GetParameter();
3703  const Double_t *cov=extpar->GetCovariance();
3704  Double_t alpha=extpar->GetAlpha();
3705  Double_t x=extpar->GetX();
3706  esdt->Set(x,alpha,par,cov);
3707  return;
3708 }
3709 //-----------------------------------------------------------------------------
3712 
3713  fMassDzero=TDatabasePDG::Instance()->GetParticle(421)->Mass();
3714  fMassDplus=TDatabasePDG::Instance()->GetParticle(411)->Mass();
3715  fMassDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
3716  fMassLambdaC=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
3717  fMassDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
3718  fMassJpsi=TDatabasePDG::Instance()->GetParticle(443)->Mass();
3719 }
3720 //-----------------------------------------------------------------------------
3722  //
3724  //
3725 
3726 
3727  //___ Check if the V0 type from AliRDHFCutsLctoV0 is the same as the one set in the ConfigVertexingHF.C for AliAnalysisVertexingHF
3728 
3729 
3731  printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHFCutsLctoV0 (%d)\n",fV0TypeForCascadeVertex,fCutsLctoV0->GetV0Type());
3732  return kFALSE;
3733  }
3734 
3736  printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHfCutsDplustoK0spi (%d)\n",fV0TypeForCascadeVertex,fCutsDplustoK0spi->GetV0Type());
3737  return kFALSE;
3738  }
3739 
3741  printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHFCutsDstoK0sK (%d)\n",fV0TypeForCascadeVertex,fCutsDstoK0sK->GetV0Type());
3742  return kFALSE;
3743  }
3744 
3745  return kTRUE;
3746 }
3747 //-----------------------------------------------------------------------------
3748 
Int_t charge
AliAODRecoDecay * fMassCalc4
for 3 prong
AliAODTrack * GetProng(AliVEvent *event, AliAODRecoDecayHF *rd, Int_t iprong)
Bool_t fUsePidTag
upper momentum limit to apply TPC PID
void SetSelectionBit(Int_t i)
selection map
double Double_t
Definition: External.C:58
void SetupPID(AliVEvent *event)
AliAnalysisFilter * fTrackFilterSoftPi
Track Filter for displaced vertices in PbPb central events (tighter cuts) for 3 prong (D+...
AliVertexerTracks * fVertexerTracks
Float_t GetMassCut(Int_t iPtBin=0) const
AliPIDResponse * fPidResponse
event mixing
AliRDHFCutsDstoK0sK * fCutsDstoK0sK
D+->Kpipi cuts.
AliAODVertex * ReconstructSecondaryVertex(TObjArray *trkArray, Double_t &dispersion, Bool_t useTRefArray=kTRUE) const
Double_t fnSigmaTOFPionHi
Low cut value on n. of sigmas for pi TOF PID.
Double_t fnSigmaTOFProtonHi
Low cut value on n. of sigmas for p TOF PID.
AliAODVertex * PrimaryVertex(const TObjArray *trkArray=0x0, AliVEvent *event=0x0) const
AliAnalysisFilter * fTrackFilter3prongCentral
Track Filter for displaced vertices in PbPb central events (tighter cuts) for 2 prong (D0->Kpi) ...
AliAnalysisVertexingHF & operator=(const AliAnalysisVertexingHF &source)
virtual void PrintAll() const
Float_t GetMassCut(Int_t iPtBin=0) const
AliAODv0 * Getv0() const
AliRDHFCutsDstoKKpi * fCutsDstoKKpi
Ds->K0s+K.
Int_t IsD0FromDStarSelected(Double_t pt, TObject *obj, Int_t selectionLevel, AliAODEvent *aod) const
Bool_t fSecVtxWithKF
z componenent of field in kG
AliRDHFCutsLctopKpi * fCutsLctopKpi
Ds->KKpi cuts.
Bool_t SelectInvMassAndPt3prong(Double_t *px, Double_t *py, Double_t *pz, Int_t pidLcStatus=3)
Int_t fnTrksTotal
triplet fullfilling Lc inv mass selection
Double_t fnSigmaTOFProtonLow
High cut value on n. of sigmas for p TPC PID.
Bool_t fOKInvMassDplus
pair fullfilling Jpsi inv mass selection
Double_t fnSigmaTPCKaonLow
High cut value on n. of sigmas for pi TOF PID.
Bool_t SingleTrkCuts(AliESDtrack *trk, Float_t centralityperc, Bool_t &okDisplaced, Bool_t &okSoftPi, Bool_t &ok3prong, Bool_t &okBachelor) const
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
Bool_t fUseTPCPIDOnlyIfNoTOF
switch use/not use TOF PID
AliRDHFCutsD0toKpi * fCutsD0toKpi
Track Filter for bachelor.
Bool_t fMixEvent
Like-sign triplets.
virtual void DeleteRecoD()
Float_t GetMassCut(Int_t iPtBin=0) const
Bool_t fUseTOFPID
switch use/not use TPC PID
Bool_t fOKInvMassD0to4p
combination fullfilling D* inv mass selection
Bool_t fOKInvMassLc
triplet fullfilling Ds inv mass selection
Float_t GetMassCut(Int_t iPtBin=0) const
Int_t * fAODMap
map between index and ID for AOD tracks
Double_t fnSigmaTPCPionLow
flag to control usage of PID tagging
Bool_t FillRecoCasc(AliVEvent *event, AliAODRecoCascadeHF *rc, Bool_t isDStar, Bool_t recoSecVtx=kFALSE)
virtual Int_t IsSelectedPID(AliAODRecoDecayHF *)
Definition: AliRDHFCuts.h:301
Bool_t fFindVertexForCascades
reconstruct a secondary vertex or assume it&#39;s from the primary vertex
Double_t fnSigmaTPCKaonHi
Low cut value on n. of sigmas for K TPC PID.
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:248
Bool_t fUseKaonPIDfor3Prong
PID response.
AliRDHFCutsDplustoK0spi * fCutsDplustoK0spi
J/psi->ee cuts.
AliRDHFCutsDplustoKpipi * fCutsDplustoKpipi
D+->K0s+pi.
AliRDHFCutsD0toKpipipi * fCutsD0toKpipipi
Lc –> v0 + bachelor cuts.
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Bool_t fLikeSign3prong
Like-sign pairs.
Bool_t SelectInvMassAndPtCascade(Double_t *px, Double_t *py, Double_t *pz)
AliAODRecoDecayHF4Prong * Make4Prong(TObjArray *fourTrackArray, AliVEvent *event, AliAODVertex *secVert, const AliAODVertex *vertexp1n1, const AliAODVertex *vertexp1n1p2, Double_t dcap1n1, Double_t dcap1n2, Double_t dcap2n1, Double_t dcap2n2, Bool_t &ok4Prong)
Class for cuts on AOD reconstructed D+->Kpipi.
Float_t GetDCACut(Int_t iPtBin=0) const
Double_t fnSigmaTPCProtonLow
High cut value on n. of sigmas for K TOF PID.
void SetParametersAtVertex(AliESDtrack *esdt, const AliExternalTrackParam *extpar) const
int Int_t
Definition: External.C:63
Float_t GetMassCut(Int_t iPtBin=0) const
Float_t GetDCACut(Int_t iPtBin=0) const
Double_t fMaxMomForTPCPid
use TPC PID only for tracks that without TOF
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
unsigned int UInt_t
Definition: External.C:33
Int_t GetIsFilled() const
float Float_t
Definition: External.C:68
Class for cuts on AOD reconstructed Ds->K0S+K.
Float_t GetMassCut(Int_t iPtBin=0) const
Float_t GetMassCut(Int_t iPtBin=0) const
AliAnalysisFilter * fTrackFilter2prongCentral
Track Filter for displaced vertices.
Double_t fnSigmaTOFKaonLow
High cut value on n. of sigmas for K TPC PID.
Float_t fMaxCentPercentileForTightCuts
High cut value on n. of sigmas for p TOF PID.
AliAODTrack * GetBachelor() const
Bool_t fOKInvMassDstar
triplet fullfilling Lc inv mass selection
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Double_t fBzkG
vertexer, to compute secondary vertices
Float_t GetDCACut(Int_t iPtBin=0) const
Bool_t fUsePIDforLc2V0
PID for Lambdac: 0=no, 1=proton, 2=p and pi.
void SetProngIDs(Int_t nIDs, UShort_t *id)
Class for cuts on AOD reconstructed D+->K0S+pi.
virtual void PrintAll() const
UShort_t GetProngID(Int_t ip) const
Float_t GetMassCut(Int_t iPtBin=0) const
Bool_t fUseKaonPIDforDs
PID for Lambdac 2 V0: 0=no, 1=proton,.
Bool_t fUseTPCPID
Kaon PID usage for Ds.
Bool_t fOKInvMassLctoV0
4tracks fullfilling D0 inv mass selection
void SetPrimaryVtxRef(TObject *vtx)
primary vertex
short Short_t
Definition: External.C:23
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
Bool_t fFindVertexForDstar
pointer to list of cuts for output file
Double_t fnSigmaTOFKaonHi
Low cut value on n. of sigmas for K TOF PID.
AliAODRecoDecay * fMassCalc2
to go faster in PbPb
Float_t GetDCACut(Int_t iPtBin=0) const
void FindCandidates(AliVEvent *event, TClonesArray *aodVerticesHFTClArr, TClonesArray *aodD0toKpiTClArr, TClonesArray *aodJPSItoEleTClArr, TClonesArray *aodCharm3ProngTClArr, TClonesArray *aodCharm4ProngTClArr, TClonesArray *aodDstarTClArr, TClonesArray *aodCascadesTClArr, TClonesArray *aodLikeSign2ProngTClArr, TClonesArray *aodLikeSign3ProngTClArr)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
void AddDaughterRefs(AliAODVertex *v, const AliVEvent *event, const TObjArray *trkArray) const
Bool_t PreSelect(TObject *obj, AliAODv0 *v0, AliVTrack *bachelorTrack)
Int_t fUsePIDforLc
Kaon PID usage for 3 prongs.
Bool_t IsEventSelected(AliVEvent *event)
Bool_t fLikeSign
cascades, Lc –> v0+track, D+ –> K0s+Pi, Ds –> K0s+K
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
void SetUsePID(Bool_t flag=kTRUE)
Definition: AliRDHFCuts.h:210
virtual void PrintAll() const
Float_t GetDCACut(Int_t iPtBin=0) const
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
AliRDHFCutsJpsitoee * fCutsJpsitoee
D0->Kpi cuts.
Bool_t SelectInvMassAndPtJpsiee(Double_t *px, Double_t *py, Double_t *pz)
Bool_t fOKInvMassD0
for 4 prong
void Setd0errProngs(Int_t nprongs, Double_t *d0)
Bool_t fOKInvMassJpsi
pair fullfilling D0 inv mass selection
Bool_t SelectInvMassAndPtDstarD0pi(Double_t *px, Double_t *py, Double_t *pz)
Int_t fV0TypeForCascadeVertex
reconstruct a secondary vertex or assume it&#39;s from the primary vertex
unsigned short UShort_t
Definition: External.C:28
AliAnalysisFilter * fTrackFilterBachelor
Track Filter for D* soft pion.
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Bool_t GetIsUsePID() const
Definition: AliRDHFCuts.h:271
Bool_t fMassCutBeforeVertexing
Select which V0 type we want to use for the cascas.
Int_t fAODMapSize
input from AOD (kTRUE) or ESD (kFALSE)
void FixReferences(AliAODEvent *aod)
void SetSelectionBitForPID(AliRDHFCuts *cuts, AliAODRecoDecayHF *rd, Int_t bit)
AliAODRecoCascadeHF * MakeCascade(TObjArray *twoTrackArray, AliVEvent *event, AliAODVertex *secVert, AliAODRecoDecayHF2Prong *rd2Prong, Double_t dca, Bool_t &okDstar)
void SelectTracksAndCopyVertex(const AliVEvent *event, Int_t trkEntries, TObjArray &seleTrksArray, TObjArray &tracksAtVertex, Int_t &nSeleTrks, UChar_t *seleFlags, Int_t *evtNumber)
bool Bool_t
Definition: External.C:53
void SetSigmaVert(Double_t sigmaVert)
Bool_t fRecoPrimVtxSkippingTrks
if kTRUE use KF vertexer, else AliVertexerTracks
AliAODRecoDecay * fMassCalc3
for 2 prong
TList * fListOfCuts
Dstar->D0pi cuts.
Bool_t RecoSecondaryVertexForCascades(AliVEvent *event, AliAODRecoCascadeHF *rc)
Float_t GetMassCut(Int_t iPtBin=0) const
Double_t fnSigmaTPCProtonHi
Low cut value on n. of sigmas for p TPC PID.
Double_t GetMinPtCandidate() const
Definition: AliRDHFCuts.h:285
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Float_t GetMassCut(Int_t iPtBin=0) const
AliRDHFCutsLctoV0 * fCutsLctoV0
Lc->pKpi cuts.
Float_t GetDCACut(Int_t iPtBin=0) const
void MapAODtracks(AliVEvent *aod)
AliAODRecoDecayHF3Prong * Make3Prong(TObjArray *threeTrackArray, AliVEvent *event, AliAODVertex *secVert, Double_t dispersion, const AliAODVertex *vertexp1n1, const AliAODVertex *vertexp2n1, Double_t dcap1n1, Double_t dcap2n1, Double_t dcap1p2, Bool_t useForLc, Bool_t useForDs, Bool_t &ok3Prong)
AliRDHFCutsDStartoKpipi * fCutsDStartoKpipi
D0->Kpipipi cuts.
AliAnalysisFilter * fTrackFilter
max. centrality percentile for using tight cuts
Float_t GetDCACut(Int_t iPtBin=0) const
Bool_t SelectInvMassAndPtD0Kpi(Double_t *px, Double_t *py, Double_t *pz)
Bool_t GetUsePID(Int_t iPtBin=0) const
AliAODv0 * TransformESDv0toAODv0(AliESDv0 *esdv0, TObjArray *twoTrackArrayV0)
Bool_t SelectInvMassAndPt4prong(Double_t *px, Double_t *py, Double_t *pz)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Double_t fnSigmaTPCPionHi
Low cut value on n. of sigmas for pi TPC PID.
Bool_t fD0toKpi
primary vertex
Bool_t fOKInvMassDs
triplet fullfilling D+ inv mass selection
void SetIsFilled(Int_t filled)
void AddRefs(AliAODVertex *v, AliAODRecoDecayHF *rd, const AliVEvent *event, const TObjArray *trkArray) const
AliAODRecoDecayHF2Prong * Make2Prong(TObjArray *twoTrackArray1, AliVEvent *event, AliAODVertex *secVert, Double_t dcap1n1, Bool_t &okD0, Bool_t &okJPSI, Bool_t &okD0fromDstar, Bool_t refill=kFALSE, AliAODRecoDecayHF2Prong *rd=0x0)
Double_t fnSigmaTOFPionLow
High cut value on n. of sigmas for pi TPC PID.