AliPhysics  608b256 (608b256)
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 //----------------------------------------------------------------------------
1706  // method to retrieve daughters from trackID and reconstruct secondary vertex
1707  // save the TRefs to the candidate AliAODRecoDecayHF3Prong rd
1708  // and fill on-the-fly the data member of rd
1709  if(rd->GetIsFilled()!=0)return kTRUE;//if 0: reduced dAOD. skip if rd is already filled (1: standard dAOD, 2 already refilled)
1710  if(!fAODMap)MapAODtracks(event);//fill the AOD index map if it is not yet done
1711  TObjArray *threeTrackArray = new TObjArray(3);
1712 
1713  AliAODTrack *track1 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(0)]);//retrieve daughter from the trackID through the AOD index map
1714  if(!track1)return kFALSE;
1715  AliAODTrack *track2 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(1)]);
1716  if(!track2)return kFALSE;
1717  AliESDtrack *postrack1 = 0;
1718  AliESDtrack *negtrack1 = 0;
1719  postrack1 = new AliESDtrack(track1);
1720  negtrack1 = new AliESDtrack(track2);
1721 
1722  // DCA between the two tracks
1723  Double_t xdummy, ydummy;
1724  fBzkG = (Double_t)event->GetMagneticField();
1725  Double_t dca12 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
1726 
1727  const AliVVertex *vprimary = event->GetPrimaryVertex();
1728  Double_t pos[3];
1729  Double_t cov[6];
1730  vprimary->GetXYZ(pos);
1731  vprimary->GetCovarianceMatrix(cov);
1732  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1733  fV1->GetCovMatrix(cov);
1734  if(!fVertexerTracks)fVertexerTracks=new AliVertexerTracks(fBzkG);
1735 
1736  AliAODTrack *track3 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(2)]);
1737  if(!track3)return kFALSE;
1738  AliESDtrack *esdt3 = new AliESDtrack(track3);
1739 
1740  Double_t dca2;
1741  Double_t dca3;
1742  threeTrackArray->AddAt(postrack1,0);
1743  threeTrackArray->AddAt(negtrack1,1);
1744  threeTrackArray->AddAt(esdt3,2);
1745  dca2 = esdt3->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
1746  dca3 = esdt3->GetDCA(postrack1,fBzkG,xdummy,ydummy);
1747  Double_t dispersion;
1748 
1749  AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray, dispersion);
1750  if (!secVert3PrAOD) {
1751  threeTrackArray->Clear();
1752  threeTrackArray->Delete(); delete threeTrackArray;
1753  delete fV1; fV1=0;
1754  delete postrack1; postrack1=NULL;
1755  delete negtrack1; negtrack1=NULL;
1756  delete esdt3; esdt3=NULL;
1757  return kFALSE;
1758  }
1759 
1760  rd->SetNProngs();
1761  Double_t vtxRec=rd->GetDist12toPrim();
1762  Double_t vertexp2n1=rd->GetDist23toPrim();
1763  rd= Make3Prong(threeTrackArray, event, secVert3PrAOD,dispersion, vtxRec, vertexp2n1, dca12, dca2, dca3, rd);
1764  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1765  rd->SetIsFilled(2);
1766  threeTrackArray->Clear();
1767  threeTrackArray->Delete(); delete threeTrackArray;
1768  delete fV1; fV1=0;
1769  delete postrack1; postrack1=NULL;
1770  delete negtrack1; negtrack1=NULL;
1771  delete esdt3; esdt3=NULL;
1772  return kTRUE;
1773 }
1774 //___________________________
1776  // method to retrieve daughters from trackID and reconstruct secondary vertex
1777  // save the TRefs to the candidate AliAODRecoDecayHF2Prong rd
1778  // and fill on-the-fly the data member of rd
1779  if(rd->GetIsFilled()!=0)return kTRUE;//if 0: reduced dAOD. skip if rd is already filled (1:standard dAOD, 2 already refilled)
1780  if(!fAODMap)MapAODtracks(event);//fill the AOD index map if it is not yet done
1781 
1782  Double_t dispersion;
1783  TObjArray *twoTrackArray1 = new TObjArray(2);
1784 
1785  AliAODTrack *track1 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(0)]);//retrieve daughter from the trackID through the AOD index map
1786  if(!track1)return kFALSE;
1787  AliAODTrack *track2 =(AliAODTrack*)event->GetTrack(fAODMap[rd->GetProngID(1)]);
1788  if(!track2)return kFALSE;
1789 
1790  AliESDtrack *esdt1 = 0;
1791  AliESDtrack *esdt2 = 0;
1792  esdt1 = new AliESDtrack(track1);
1793  esdt2 = new AliESDtrack(track2);
1794 
1795  twoTrackArray1->AddAt(esdt1,0);
1796  twoTrackArray1->AddAt(esdt2,1);
1797  // DCA between the two tracks
1798  Double_t xdummy, ydummy;
1799  fBzkG = (Double_t)event->GetMagneticField();
1800  Double_t dca12 = esdt1->GetDCA(esdt2,fBzkG,xdummy,ydummy);
1801  const AliVVertex *vprimary = event->GetPrimaryVertex();
1802  Double_t pos[3];
1803  Double_t cov[6];
1804  vprimary->GetXYZ(pos);
1805  vprimary->GetCovarianceMatrix(cov);
1806  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1807  fV1->GetCovMatrix(cov);
1808  if(!fVertexerTracks)fVertexerTracks=new AliVertexerTracks(fBzkG);
1809 
1810 
1811  AliAODVertex *vtxRec = ReconstructSecondaryVertex(twoTrackArray1, dispersion);
1812  if(!vtxRec) {
1813  twoTrackArray1->Clear();
1814  twoTrackArray1->Delete(); delete twoTrackArray1;
1815  delete fV1; fV1=0;
1816  delete esdt1; esdt1=NULL;
1817  delete esdt2; esdt2=NULL;
1818  return kFALSE; }
1819  Bool_t okD0=kFALSE;
1820  Bool_t okJPSI=kFALSE;
1821  Bool_t okD0FromDstar=kFALSE;
1822  Bool_t refill =kTRUE;
1823  rd->SetNProngs();
1824  rd= Make2Prong(twoTrackArray1, event, vtxRec, dca12, okD0, okJPSI, okD0FromDstar,refill,rd);
1825  rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1826  rd->SetIsFilled(2);
1827  delete fV1; fV1=0;
1828  twoTrackArray1->Clear();
1829  twoTrackArray1->Delete(); delete twoTrackArray1;
1830  delete esdt1; esdt1=NULL;
1831  delete esdt2; esdt2=NULL;
1832  return kTRUE;
1833 }
1834 //----------------------------------------------------------------------------
1836  // method to retrieve daughters from trackID
1837  // and fill on-the-fly the data member of rCasc and their AliAODRecoDecayHF2Prong daughters
1838  if(rCasc->GetIsFilled()!=0) return kTRUE;//if 0: reduced dAOD. skip if rd is already filled (1: standard dAOD, 2: already refilled)
1839  if(!fAODMap)MapAODtracks(event);//fill the AOD index map if it is not yet done
1840  TObjArray *twoTrackArrayCasc = new TObjArray(2);
1841 
1842  AliAODTrack *trackB =(AliAODTrack*)event->GetTrack(fAODMap[rCasc->GetProngID(0)]);//retrieve bachelor from the trackID through the AOD index map
1843  if(!trackB)return kFALSE;
1844 
1845  AliNeutralTrackParam *trackV0=NULL;
1846  AliAODv0 *v0 =NULL;
1847  AliAODRecoDecayHF2Prong *trackD0=NULL;
1848  rCasc->SetNProngs();
1849 
1850  if(DStar){
1851  TClonesArray *inputArrayD0=(TClonesArray*)event->GetList()->FindObject("D0toKpi");
1852  if(!inputArrayD0) return kFALSE;
1853  trackD0=(AliAODRecoDecayHF2Prong*)inputArrayD0->At(rCasc->GetProngID(1));
1854  if(!trackD0)return kFALSE;
1855  if(!FillRecoCand(event,trackD0)) return kFALSE; //fill missing info of the corresponding D0 daughter
1856 
1857  trackV0 = new AliNeutralTrackParam(trackD0);
1858 
1859  }else{//is a V0 candidate
1860  v0 = ((AliAODEvent*)event)->GetV0(rCasc->GetProngID(1));
1861  if(!v0) return kFALSE;
1862  // Define the V0 (neutral) track
1863  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
1864  if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
1865  }
1866 
1867  AliESDtrack *esdB = new AliESDtrack(trackB);
1868 
1869  twoTrackArrayCasc->AddAt(esdB,0);
1870  twoTrackArrayCasc->AddAt(trackV0,1);
1871 
1872  fBzkG = (Double_t)event->GetMagneticField();
1873  const AliVVertex *vprimary = event->GetPrimaryVertex();
1874 
1875  Double_t pos[3];
1876  Double_t cov[6];
1877  vprimary->GetXYZ(pos);
1878  vprimary->GetCovarianceMatrix(cov);
1879  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1880  fV1->GetCovMatrix(cov);
1881 
1882  Double_t dca = 0.;
1883  AliAODVertex *vtxCasc = 0x0;
1884  Double_t chi2perNDF = fV1->GetChi2toNDF();
1885  if (recoSecVtx) {
1886  Double_t dispersion, xdummy, ydummy;
1887  dca = esdB->GetDCA(trackV0,fBzkG,xdummy,ydummy);
1888  if (!fVertexerTracks) fVertexerTracks = new AliVertexerTracks(fBzkG);
1889  vtxCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
1890  } else {
1891  vtxCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
1892  }
1893  if(!vtxCasc) {
1894  twoTrackArrayCasc->Clear();
1895  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1896  delete fV1; fV1=0;
1897  delete esdB; esdB=NULL;
1898  delete vtxCasc;vtxCasc=NULL;
1899  trackB=NULL;
1900  delete trackV0; trackV0=NULL;
1901  if(!DStar){
1902  v0=NULL;
1903  }
1904  return kFALSE;
1905  }
1906  vtxCasc->SetParent(rCasc);
1907  rCasc->SetSecondaryVtx(vtxCasc);
1908  AddDaughterRefs(vtxCasc,(AliAODEvent*)event,twoTrackArrayCasc);
1909  if(DStar)vtxCasc->AddDaughter(trackD0);
1910  else vtxCasc->AddDaughter(v0);
1911  rCasc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1912 
1913  Bool_t refill =kTRUE;
1914 
1915  Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1916  // propagate tracks to secondary vertex, to compute inv. mass
1917  esdB->PropagateToDCA(vtxCasc,fBzkG,kVeryBig);
1918  trackV0->PropagateToDCA(vtxCasc,fBzkG,kVeryBig);
1919  Double_t momentum[3];
1920  esdB->GetPxPyPz(momentum);
1921  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1922  trackV0->GetPxPyPz(momentum);
1923  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1924 
1925  AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArrayCasc,event);
1926  if(!primVertexAOD){
1927  delete fV1; fV1=0;
1928  delete vtxCasc; vtxCasc=NULL;
1929  twoTrackArrayCasc->Clear();
1930  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1931  delete esdB; esdB=NULL;
1932  delete trackV0; trackV0=NULL;
1933  if(!DStar)v0=NULL;
1934  return kFALSE;
1935  }
1936  Double_t d0z0[2],covd0z0[3];
1937  esdB->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1938  d0[0] = d0z0[0];
1939  d0err[0] = TMath::Sqrt(covd0z0[0]);
1940  trackV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1941  d0[1] = d0z0[0];
1942  d0err[1] = TMath::Sqrt(covd0z0[0]);
1943  rCasc->SetPxPyPzProngs(2,px,py,pz);
1944  rCasc->SetDCA(dca);
1945  rCasc->Setd0Prongs(2,d0);
1946  rCasc->Setd0errProngs(2,d0err);
1947  rCasc->SetOwnPrimaryVtx(primVertexAOD);
1948  rCasc->SetCharge(esdB->Charge());
1949  // get PID info from ESD
1950  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1951  if(esdB->GetStatus()&AliESDtrack::kESDpid) esdB->GetESDpid(esdpid0);
1952  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1953  Double_t esdpid[10];
1954  for(Int_t i=0;i<5;i++) {
1955  esdpid[i] = esdpid0[i];
1956  esdpid[5+i] = esdpid1[i];
1957  }
1958  rCasc->SetPID(2,esdpid);
1959  rCasc->SetIsFilled(2);
1960 
1961 
1962  delete fV1; fV1=0;
1963  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1964  twoTrackArrayCasc->Clear();
1965  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1966  delete esdB; esdB=NULL;
1967  delete trackV0; trackV0=NULL;
1968  if(!DStar)v0=NULL;
1969  return kTRUE;
1970 }
1971 //---------------------------------------------------------------------------
1973 {
1978 
1979 
1980  if (rc->GetIsFilled()==0) return kFALSE; // Reduced dAOD with candidates not yet refilled
1981  if (!fAODMap) MapAODtracks(event); // Fill the AOD index map if it is not yet done
1982 
1983 
1984  // - Get bachelor and V0 from the cascade
1985  AliAODTrack *trackB = dynamic_cast<AliAODTrack*> (rc->GetBachelor());
1986  if (!trackB) return kFALSE;
1987 
1988  AliAODv0 *v0 = dynamic_cast<AliAODv0*> (rc->Getv0());
1989  if (!v0) return kFALSE;
1990 
1991 
1992  // - Cast bachelor (AOD->ESD) and V0 (AliAODv0->AliNeutralTrackParam)
1993  AliESDtrack *esdB = new AliESDtrack(trackB);
1994  if (!esdB) return kFALSE;
1995 
1996  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
1997  if (!trackVV0) return kFALSE;
1998 
1999  AliNeutralTrackParam *trackV0 = new AliNeutralTrackParam(trackVV0);
2000  if (!trackV0) return kFALSE;
2001 
2002  TObjArray *twoTrackArrayCasc = new TObjArray(2);
2003  twoTrackArrayCasc->AddAt(esdB, 0);
2004  twoTrackArrayCasc->AddAt(trackV0, 1);
2005 
2006 
2007  Double_t dispersion, xdummy, ydummy, pos[3], cov[6];
2008 
2009 
2010  // - Some ingredients will be needed:
2011  // magnetic field
2012  fBzkG = (Double_t) event->GetMagneticField();
2013 
2014  // primary vertex
2015  const AliVVertex *vprimary = event->GetPrimaryVertex();
2016  vprimary->GetXYZ(pos);
2017  vprimary->GetCovarianceMatrix(cov);
2018  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2019 
2020  // vertexer
2021  if (!fVertexerTracks) fVertexerTracks = new AliVertexerTracks(fBzkG);
2022 
2023 
2024  // - Compute the DCA and the decay vertex
2025  Double_t dca = esdB->GetDCA(trackV0, fBzkG, xdummy, ydummy);
2026  AliAODVertex *vtxCasc = ReconstructSecondaryVertex(twoTrackArrayCasc, dispersion, kFALSE);
2027 
2028  if (!vtxCasc) {
2029  twoTrackArrayCasc->Clear();
2030  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
2031  delete fV1; fV1=0;
2032  delete esdB; esdB=0;
2033  delete vtxCasc; vtxCasc=0;
2034  delete trackV0; trackV0=0;
2035  trackB=0; v0=0;
2036  return kFALSE;
2037  }
2038 
2039 
2040  // - Linking the cascade with the new secondary vertex
2041  vtxCasc->SetParent(rc);
2042  rc->SetSecondaryVtx(vtxCasc);
2043  AddDaughterRefs(vtxCasc, (AliAODEvent*)event, twoTrackArrayCasc);
2044  vtxCasc->AddDaughter(v0);
2045  rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
2046 
2047 
2048  // - Propagate tracks to secondary vertex, to get momenta
2049  Double_t momentum[3], px[2], py[2], pz[2], d0[2], d0err[2], d0z0[2], covd0z0[3];
2050  esdB->PropagateToDCA(vtxCasc, fBzkG, kVeryBig);
2051  trackV0->PropagateToDCA(vtxCasc, fBzkG, kVeryBig);
2052  esdB->GetPxPyPz(momentum);
2053  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2054  trackV0->GetPxPyPz(momentum);
2055  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2056 
2057 
2058  // - Propagate tracks to primary vertex, to get impact parameters
2059  AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArrayCasc, event);
2060  if (!primVertexAOD) {
2061  twoTrackArrayCasc->Clear();
2062  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
2063  delete fV1; fV1=0;
2064  delete esdB; esdB=0;
2065  delete vtxCasc; vtxCasc=0;
2066  delete trackV0; trackV0=0;
2067  trackB=0; v0=0;
2068  return kFALSE;
2069  }
2070 
2071  esdB->PropagateToDCA(primVertexAOD, fBzkG, kVeryBig, d0z0, covd0z0);
2072  d0[0] = d0z0[0];
2073  d0err[0] = TMath::Sqrt(covd0z0[0]);
2074  trackV0->PropagateToDCA(primVertexAOD, fBzkG, kVeryBig, d0z0, covd0z0);
2075  d0[1] = d0z0[0];
2076  d0err[1] = TMath::Sqrt(covd0z0[0]);
2077 
2078 
2079  // - Get PID info from ESD
2080  Double_t esdpid[10];
2081  Double_t esdpid0[5] = {0., 0., 0., 0., 0.};
2082  Double_t esdpid1[5] = {0., 0., 0., 0., 0.};
2083  if (esdB->GetStatus()&AliESDtrack::kESDpid) esdB->GetESDpid(esdpid0);
2084  for (Int_t ipid=0; ipid<5; ipid++) {
2085  esdpid[ipid] = esdpid0[ipid];
2086  esdpid[5+ipid] = esdpid1[ipid];
2087  }
2088 
2089 
2090  // - Set data members
2091  rc->SetOwnPrimaryVtx(primVertexAOD);
2092  rc->SetDCA(dca);
2093  rc->SetPxPyPzProngs(2, px, py, pz);
2094  rc->Setd0Prongs(2, d0);
2095  rc->Setd0errProngs(2, d0err);
2096  rc->SetCharge(esdB->Charge());
2097  rc->SetPID(2, esdpid);
2098 
2099  rc->SetIsFilled(2); // To clean the newly reconstructed secondary vertex with the CleanUp task
2100 
2101 
2102  twoTrackArrayCasc->Clear();
2103  twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
2104  delete fV1; fV1=0;
2105  delete primVertexAOD; primVertexAOD=0;
2106  delete esdB; esdB=0;
2107  delete trackV0; trackV0=0;
2108  trackB=0; v0=0;
2109 
2110  return kTRUE;
2111 }
2112 
2113 //---------------------------------------------------------------------------
2115  TObjArray *twoTrackArray,AliVEvent *event,
2116  AliAODVertex *secVert,
2117  AliAODRecoDecayHF2Prong *rd2Prong,
2118  Double_t dca,
2119  Bool_t &okDstar)
2120 {
2123  //AliCodeTimerAuto("",0);
2124  UInt_t ntref=TProcessID::GetObjectCount();
2125  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2126  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2127  }
2128 
2129  okDstar = kFALSE;
2130 
2131  Bool_t dummy1,dummy2,dummy3;
2132 
2133  // We use Make2Prong to construct the AliAODRecoCascadeHF
2134  // (which inherits from AliAODRecoDecayHF2Prong)
2135  AliAODRecoCascadeHF *theCascade =
2136  (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
2137  dummy1,dummy2,dummy3);
2138  if(!theCascade) return 0x0;
2139 
2140  // charge
2141  AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2142  theCascade->SetCharge(trackPi->Charge());
2143 
2144  //--- selection cuts
2145  //
2146  AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
2147  if(fInputAOD){
2148  Int_t idSoftPi=(Int_t)trackPi->GetID();
2149  if (idSoftPi > -1 && idSoftPi < fAODMapSize) {
2150  AliAODTrack* trackPiAOD=dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[idSoftPi]));
2151  if(!trackPiAOD) AliFatal("Not a standard AOD");
2152  tmpCascade->GetSecondaryVtx()->AddDaughter(trackPiAOD);
2153  }
2154  }else{
2155  tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
2156  }
2157  tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
2158 
2159  AliAODVertex *primVertexAOD=0;
2161  // take event primary vertex
2162  primVertexAOD = PrimaryVertex();
2163  tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
2164  rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
2165  }
2166  // select D*->D0pi
2167  if(fDstar) {
2169  if(okDstar) theCascade->SetSelectionBit(AliRDHFCuts::kDstarCuts);
2170  }
2171  tmpCascade->GetSecondaryVtx()->RemoveDaughters();
2172  tmpCascade->UnsetOwnPrimaryVtx();
2173  delete tmpCascade; tmpCascade=NULL;
2175  rd2Prong->UnsetOwnPrimaryVtx();
2176  }
2177  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2178  //---
2179 
2180 
2181  return theCascade;
2182 }
2183 
2184 
2185 //----------------------------------------------------------------------------
2187  TObjArray *twoTrackArray,AliVEvent *event,
2188  AliAODVertex *secVert,
2189  AliAODv0 *v0,
2190  Double_t dca,
2191  Bool_t &okCascades)
2192 {
2195  //AliCodeTimerAuto("",0);
2196  UInt_t ntref=TProcessID::GetObjectCount();
2197  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2198  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2199  }
2200 
2201  // preselection to reduce the number of TRefs
2202  Double_t px[2],py[2],pz[2];
2203  AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2204  AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
2205  // propagate tracks to secondary vertex, to compute inv. mass
2206  postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2207  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2208  Double_t momentum[3];
2209  postrack->GetPxPyPz(momentum);
2210  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2211  negtrack->GetPxPyPz(momentum);
2212  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2213  if(!SelectInvMassAndPtCascade(px,py,pz)) return 0x0;
2214  Double_t dummyd0[2]={0,0};
2215  Double_t dummyd0err[2]={0,0};
2216  AliAODRecoCascadeHF tmpCasc(0x0,postrack->Charge(),px, py, pz, dummyd0, dummyd0err,0.);
2217  Bool_t presel=fCutsLctoV0->PreSelect(&tmpCasc,v0,postrack);
2218  if(!presel) return 0x0;
2219 
2220  // AliDebug(2,Form(" building the cascade"));
2221  okCascades = kFALSE;
2222  Bool_t dummy1,dummy2,dummy3;
2223 
2224  // We use Make2Prong to construct the AliAODRecoCascadeHF
2225  // (which inherits from AliAODRecoDecayHF2Prong)
2226  AliAODRecoCascadeHF *theCascade =
2227  (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
2228  dummy1,dummy2,dummy3);
2229  if(!theCascade) return 0x0;
2230 
2231  // bachelor track and charge
2232  AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2233  theCascade->SetCharge(trackBachelor->Charge());
2234 
2235  //--- selection cuts
2236  //
2237 
2238  AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
2239  if(fInputAOD){
2240  Int_t idBachelor=(Int_t)trackBachelor->GetID();
2241  if (idBachelor > -1 && idBachelor < fAODMapSize) {
2242  AliAODTrack* trackBachelorAOD=dynamic_cast<AliAODTrack*>(event->GetTrack(fAODMap[idBachelor]));
2243  if(!trackBachelorAOD) AliFatal("Not a standard AOD");
2244  tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelorAOD);
2245  }
2246  }else{
2247  tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
2248  }
2249  tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
2250 
2251  AliAODVertex *primVertexAOD=0;
2253  // take event primary vertex
2254  primVertexAOD = PrimaryVertex();
2255  if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
2256  tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
2257  }
2258 
2259  // select Cascades
2260  if (fCascades && fInputAOD) {
2261  if (fCutsLctoV0->IsSelected(tmpCascade, AliRDHFCuts::kCandidate)>0) {
2262  okCascades = kTRUE;
2264  }
2265  if (fCutsDplustoK0spi) {
2267  okCascades = kTRUE;
2269  }
2270  }
2271  if (fCutsDstoK0sK){
2272  if (fCutsDstoK0sK->IsSelected(tmpCascade, AliRDHFCuts::kCandidate)>0) {
2273  okCascades = kTRUE;
2275  }
2276  }
2277  }
2278  else {
2279  //AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied"));
2280  okCascades=kTRUE;
2281  } // no cuts implemented from ESDs
2282  tmpCascade->GetSecondaryVtx()->RemoveDaughters();
2283  tmpCascade->UnsetOwnPrimaryVtx();
2284  delete tmpCascade; tmpCascade=NULL;
2285  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2286  //---
2287 
2288  return theCascade;
2289 }
2290 
2291 //-----------------------------------------------------------------------------
2293  TObjArray *twoTrackArray,AliVEvent *event,
2294  AliAODVertex *secVert,Double_t dca,
2295  Bool_t &okD0,Bool_t &okJPSI,
2296  Bool_t &okD0fromDstar, Bool_t refill, AliAODRecoDecayHF2Prong *rd)
2297 {
2300  // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
2301  // AliCodeTimerAuto("",0);
2302  okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
2303 
2304  UInt_t ntref=TProcessID::GetObjectCount();
2305  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2306  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2307  }
2308 
2309  Double_t px[2],py[2],pz[2],d0[2],d0err[2];
2310  AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
2311  AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
2312 
2313  // propagate tracks to secondary vertex, to compute inv. mass
2314  postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2315  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2316 
2317  Double_t momentum[3];
2318  postrack->GetPxPyPz(momentum);
2319  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2320  negtrack->GetPxPyPz(momentum);
2321  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2322 
2323  if(!refill){//skip if it is called in refill step because already checked
2324  // invariant mass cut (try to improve coding here..)
2325  Bool_t okMassCut=kFALSE;
2326  if(!okMassCut && fD0toKpi) if(SelectInvMassAndPtD0Kpi(px,py,pz)) okMassCut=kTRUE;
2327  if(!okMassCut && fJPSItoEle) if(SelectInvMassAndPtJpsiee(px,py,pz)) okMassCut=kTRUE;
2328  if(!okMassCut && fDstar) if(SelectInvMassAndPtDstarD0pi(px,py,pz)) okMassCut=kTRUE;
2329  if(!okMassCut && fCascades) if(SelectInvMassAndPtCascade(px,py,pz)) okMassCut=kTRUE;
2330  if(!okMassCut) {
2331  //AliDebug(2," candidate didn't pass mass cut");
2332  return 0x0;
2333  }
2334  }
2335  // primary vertex to be used by this candidate
2336  AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
2337  if(!primVertexAOD) return 0x0;
2338 
2339  Double_t d0z0[2],covd0z0[3];
2340  postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2341  d0[0] = d0z0[0];
2342  d0err[0] = TMath::Sqrt(covd0z0[0]);
2343  negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2344  d0[1] = d0z0[0];
2345  d0err[1] = TMath::Sqrt(covd0z0[0]);
2346  AliAODRecoDecayHF2Prong *the2Prong;
2347  // create the object AliAODRecoDecayHF2Prong
2348  if(!refill){
2349  the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
2350  the2Prong->SetOwnPrimaryVtx(primVertexAOD);
2351  UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
2352  the2Prong->SetProngIDs(2,id);
2353  if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
2354  // Add daughter references already here
2355  if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
2356 
2357  // select D0->Kpi
2358  if(fD0toKpi) {
2360  if(okD0) the2Prong->SetSelectionBit(AliRDHFCuts::kD0toKpiCuts);
2361  }
2362  //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
2363  // select J/psi from B
2364  if(fJPSItoEle) {
2366  }
2367  //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
2368  // select D0->Kpi from Dstar
2369  if(fDstar) {
2370  okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
2371  if(okD0fromDstar) the2Prong->SetSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
2372  }
2373  //if(fDebug && fDstar) printf(" %d\n",(Int_t)okD0fromDstar);
2374  }
2375  }else{
2376  the2Prong =rd;
2377  the2Prong->SetSecondaryVtx(secVert);
2378  secVert->SetParent(the2Prong);
2379  AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
2380  the2Prong->SetOwnPrimaryVtx(primVertexAOD);
2381  the2Prong->SetPxPyPzProngs(2,px,py,pz);
2382  the2Prong->SetDCA(dca);
2383  the2Prong->Setd0Prongs(2,d0);
2384  the2Prong->Setd0errProngs(2,d0err);
2385  the2Prong->SetCharge(0);
2386  }
2387  delete primVertexAOD; primVertexAOD=NULL;
2388 
2389  // remove the primary vertex (was used only for selection)
2391  the2Prong->UnsetOwnPrimaryVtx();
2392  }
2393 
2394  // get PID info from ESD
2395  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2396  if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
2397  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2398  if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2399  Double_t esdpid[10];
2400  for(Int_t i=0;i<5;i++) {
2401  esdpid[i] = esdpid0[i];
2402  esdpid[5+i] = esdpid1[i];
2403  }
2404  the2Prong->SetPID(2,esdpid);
2405  return the2Prong;
2406 }
2407 //----------------------------------------------------------------------------
2409  TObjArray *threeTrackArray,AliVEvent *event,
2410  AliAODVertex *secVert,Double_t dispersion,
2411  const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
2412  Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
2413  Bool_t useForLc, Bool_t useForDs, Bool_t &ok3Prong)
2414 {
2417  // E.Bruna, F.Prino
2418  // AliCodeTimerAuto("",0);
2419 
2420 
2421  UInt_t ntref=TProcessID::GetObjectCount();
2422  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2423  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2424  }
2425 
2426  ok3Prong=kFALSE;
2427  if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
2428 
2429  Double_t px[3],py[3],pz[3],d0[3],d0err[3];
2430  Double_t momentum[3];
2431 
2432 
2433  AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
2434  AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
2435  AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
2436 
2437  postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2438  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2439  postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2440  postrack1->GetPxPyPz(momentum);
2441  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2442  negtrack->GetPxPyPz(momentum);
2443  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2444  postrack2->GetPxPyPz(momentum);
2445  px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
2446 
2447  // invariant mass cut for D+, Ds, Lc
2448  Bool_t okMassCut=kFALSE;
2449  if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
2450  if(!okMassCut && f3Prong) if(SelectInvMassAndPt3prong(px,py,pz)) okMassCut=kTRUE;
2451  if(!okMassCut) {
2452  //AliDebug(2," candidate didn't pass mass cut");
2453  return 0x0;
2454  }
2455 
2456  // primary vertex to be used by this candidate
2457  AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
2458  if(!primVertexAOD) return 0x0;
2459 
2460  Double_t d0z0[2],covd0z0[3];
2461  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2462  d0[0]=d0z0[0];
2463  d0err[0] = TMath::Sqrt(covd0z0[0]);
2464  negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2465  d0[1]=d0z0[0];
2466  d0err[1] = TMath::Sqrt(covd0z0[0]);
2467  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2468  d0[2]=d0z0[0];
2469  d0err[2] = TMath::Sqrt(covd0z0[0]);
2470 
2471 
2472  // create the object AliAODRecoDecayHF3Prong
2473  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2474  Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
2475  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]));
2476  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]));
2477  Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
2478 
2479 
2480  // construct the candidate passing a NULL pointer for the secondary vertex to avoid creation of TRef
2481  AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(0x0,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
2482  the3Prong->SetOwnPrimaryVtx(primVertexAOD);
2483  // add a pointer to the secondary vertex via SetOwnSecondaryVtx (no TRef created)
2484  AliAODVertex* ownsecv=secVert->CloneWithoutRefs();
2485  the3Prong->SetOwnSecondaryVtx(ownsecv);
2486  UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
2487  the3Prong->SetProngIDs(3,id);
2488 
2489  delete primVertexAOD; primVertexAOD=NULL;
2490 
2491  // disable PID, which requires the TRefs to the daughter tracks
2492  fCutsDplustoKpipi->SetUsePID(kFALSE);
2493  fCutsDstoKKpi->SetUsePID(kFALSE);
2494  fCutsLctopKpi->SetUsePID(kFALSE);
2495 
2496  // select D+->Kpipi, Ds->KKpi, Lc->pKpi
2497  if(f3Prong) {
2498  ok3Prong = kFALSE;
2499 
2501  ok3Prong = kTRUE;
2503  }
2504  if(useForDs && fOKInvMassDs){
2506  ok3Prong = kTRUE;
2508  }
2509  }
2510  if(useForLc && fOKInvMassLc){
2512  ok3Prong = kTRUE;
2514  }
2515  }
2516  }
2517  //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
2518  the3Prong->UnsetOwnSecondaryVtx();
2520  the3Prong->UnsetOwnPrimaryVtx();
2521  }
2522 
2523  // Add TRefs to secondary vertex and daughter tracks only for candidates passing the filtering cuts
2524  if(ok3Prong && fInputAOD){
2525  the3Prong->SetSecondaryVtx(secVert);
2526  AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
2527  }
2528 
2529  // get PID info from ESD
2530  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2531  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2532  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2533  if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2534  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2535  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2536 
2537  Double_t esdpid[15];
2538  for(Int_t i=0;i<5;i++) {
2539  esdpid[i] = esdpid0[i];
2540  esdpid[5+i] = esdpid1[i];
2541  esdpid[10+i] = esdpid2[i];
2542  }
2543  the3Prong->SetPID(3,esdpid);
2544 
2545  return the3Prong;
2546 }
2547 //----------------------------------------------------------------------------
2549  TObjArray *threeTrackArray,AliVEvent *event,
2550  AliAODVertex *secVert,Double_t dispersion,
2551  Double32_t dist12, Double32_t dist23,
2552  Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
2554 {
2555  // Fill on-the-fly the 3prong data member missing info
2556  // set TRef of vertex and daughters
2557  // do not recalculate the two-track secondary vertex dist12 and dist23
2558  // because they are stored in the AliAODRecoDecayHF3Prong candidate
2559  // reconstructed in the FindCandidate step
2560  // do not check if it is a Lambdac, Ds or Dplus because it is already check in the FindCandidate step
2561  // and the info is stored
2562  // AliCodeTimerAuto("",0);
2563 
2564  UInt_t ntref=TProcessID::GetObjectCount();
2565  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2566  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2567  }
2568 
2569  Double_t px[3],py[3],pz[3],d0[3],d0err[3];
2570  Double_t momentum[3];
2571 
2572 
2573  AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
2574  AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
2575  AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
2576 
2577  postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2578  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
2579  postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2580  postrack1->GetPxPyPz(momentum);
2581  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2582  negtrack->GetPxPyPz(momentum);
2583  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2584  postrack2->GetPxPyPz(momentum);
2585  px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
2586  // primary vertex to be used by this candidate
2587  AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
2588  if(!primVertexAOD) return 0x0;
2589  Double_t d0z0[2],covd0z0[3];
2590  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2591  d0[0]=d0z0[0];
2592  d0err[0] = TMath::Sqrt(covd0z0[0]);
2593  negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2594  d0[1]=d0z0[0];
2595  d0err[1] = TMath::Sqrt(covd0z0[0]);
2596  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2597  d0[2]=d0z0[0];
2598  d0err[2] = TMath::Sqrt(covd0z0[0]);
2599 
2600  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2601  Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
2602  Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
2603 
2604  rd->SetSecondaryVtx(secVert);
2605  secVert->SetParent(rd);
2606  AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
2607 
2608  rd->SetPxPyPzProngs(3,px,py,pz);
2609  rd->SetDCAs(3,dca);
2610  rd->Setd0Prongs(3,d0);
2611  rd->Setd0errProngs(3,d0err);
2612  rd->SetCharge(charge);
2613  rd->SetOwnPrimaryVtx(primVertexAOD);
2614  rd->SetSigmaVert(dispersion);
2615  delete primVertexAOD; primVertexAOD=NULL;
2616 
2618  rd->UnsetOwnPrimaryVtx();
2619  }
2620 
2621  // get PID info from ESD
2622  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2623  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2624  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2625  if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
2626  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2627  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2628 
2629  Double_t esdpid[15];
2630  for(Int_t i=0;i<5;i++) {
2631  esdpid[i] = esdpid0[i];
2632  esdpid[5+i] = esdpid1[i];
2633  esdpid[10+i] = esdpid2[i];
2634  }
2635  rd->SetPID(3,esdpid);
2636  return rd;
2637 }
2638 //----------------------------------------------------------------------------
2640  TObjArray *fourTrackArray,AliVEvent *event,
2641  AliAODVertex *secVert,
2642  const AliAODVertex *vertexp1n1,
2643  const AliAODVertex *vertexp1n1p2,
2644  Double_t dcap1n1,Double_t dcap1n2,
2645  Double_t dcap2n1,Double_t dcap2n2,
2646  Bool_t &ok4Prong)
2647 {
2650  // G.E.Bruno, R.Romita
2651  // AliCodeTimerAuto("",0);
2652 
2653  UInt_t ntref=TProcessID::GetObjectCount();
2654  if(ntref>16776216){// This number is 2^24-1000. The maximum number of TRef for a given TProcesssID is 2^24=16777216.
2655  AliFatal(Form("Number of TRef created (=%d), close to limit (16777216)",ntref));
2656  }
2657 
2658  ok4Prong=kFALSE;
2659  if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
2660 
2661  Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
2662 
2663  AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
2664  AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
2665  AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
2666  AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
2667 
2668  postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2669  negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
2670  postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2671  negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
2672 
2673  Double_t momentum[3];
2674  postrack1->GetPxPyPz(momentum);
2675  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
2676  negtrack1->GetPxPyPz(momentum);
2677  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
2678  postrack2->GetPxPyPz(momentum);
2679  px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
2680  negtrack2->GetPxPyPz(momentum);
2681  px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
2682 
2683  // invariant mass cut for rho or D0 (try to improve coding here..)
2684  Bool_t okMassCut=kFALSE;
2685  if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
2686  if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) { //no PID, to be implemented with PID
2687  if(SelectInvMassAndPt4prong(px,py,pz)) okMassCut=kTRUE;
2688  }
2689  if(!okMassCut) {
2690  //if(fDebug) printf(" candidate didn't pass mass cut\n");
2691  //printf(" candidate didn't pass mass cut\n");
2692  return 0x0;
2693  }
2694 
2695  // primary vertex to be used by this candidate
2696  AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
2697  if(!primVertexAOD) return 0x0;
2698 
2699  Double_t d0z0[2],covd0z0[3];
2700  postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2701  d0[0]=d0z0[0];
2702  d0err[0] = TMath::Sqrt(covd0z0[0]);
2703  negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2704  d0[1]=d0z0[0];
2705  d0err[1] = TMath::Sqrt(covd0z0[0]);
2706  postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2707  d0[2]=d0z0[0];
2708  d0err[2] = TMath::Sqrt(covd0z0[0]);
2709  negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2710  d0[3]=d0z0[0];
2711  d0err[3] = TMath::Sqrt(covd0z0[0]);
2712 
2713 
2714  // create the object AliAODRecoDecayHF4Prong
2715  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
2716  Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
2717  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]));
2718  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]));
2719  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]));
2720  Short_t charge=0;
2721  AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
2722  the4Prong->SetOwnPrimaryVtx(primVertexAOD);
2723  UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
2724  the4Prong->SetProngIDs(4,id);
2725 
2726  delete primVertexAOD; primVertexAOD=NULL;
2727 
2729 
2730 
2732  the4Prong->UnsetOwnPrimaryVtx();
2733  }
2734 
2735 
2736  // get PID info from ESD
2737  Double_t esdpid0[5]={0.,0.,0.,0.,0.};
2738  if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
2739  Double_t esdpid1[5]={0.,0.,0.,0.,0.};
2740  if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
2741  Double_t esdpid2[5]={0.,0.,0.,0.,0.};
2742  if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
2743  Double_t esdpid3[5]={0.,0.,0.,0.,0.};
2744  if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
2745 
2746  Double_t esdpid[20];
2747  for(Int_t i=0;i<5;i++) {
2748  esdpid[i] = esdpid0[i];
2749  esdpid[5+i] = esdpid1[i];
2750  esdpid[10+i] = esdpid2[i];
2751  esdpid[15+i] = esdpid3[i];
2752  }
2753  the4Prong->SetPID(4,esdpid);
2754 
2755  return the4Prong;
2756 }
2757 //----------------------------------------------------------------------------------
2759  //assign and save in fAODMap the index of the AliAODTrack track
2760  //ordering them on the basis of selected criteria
2761 
2762  fAODMapSize = 100000;
2763  fAODMap = new Int_t[fAODMapSize];
2764  AliAODTrack *track=0;
2765  memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
2766  for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
2767  track = dynamic_cast<AliAODTrack*>(aod->GetTrack(i));
2768  if(!track) AliFatal("Not a standard AOD");
2769  // skip pure ITS SA tracks
2770  if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2771 
2772  // skip tracks without ITS
2773  if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
2774 
2775  // TEMPORARY: check that the cov matrix is there
2776  Double_t covtest[21];
2777  if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2778  //
2779 
2780  Int_t ind = (Int_t)track->GetID();
2781  if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
2782  }
2783  return;
2784 }
2785 //-----------------------------------------------------------------------------
2786 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
2787  AliVEvent *event) const
2788 {
2790  //AliCodeTimerAuto("",0);
2791 
2792  AliESDVertex *vertexESD = 0;
2793  AliAODVertex *vertexAOD = 0;
2794 
2795 
2797  // primary vertex from the input event
2798 
2799  vertexESD = new AliESDVertex(*fV1);
2800 
2801  } else {
2802  // primary vertex specific to this candidate
2803 
2804  Int_t nTrks = trkArray->GetEntriesFast();
2805  AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
2806 
2808  // recalculating the vertex
2809 
2810  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
2811  Float_t diamondcovxy[3];
2812  event->GetDiamondCovXY(diamondcovxy);
2813  Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
2814  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
2815  AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
2816  vertexer->SetVtxStart(diamond);
2817  delete diamond; diamond=NULL;
2818  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
2819  vertexer->SetOnlyFitter();
2820  }
2821  Int_t skipped[1000];
2822  Int_t nTrksToSkip=0,id;
2823  AliExternalTrackParam *t = 0;
2824  for(Int_t i=0; i<nTrks; i++) {
2825  t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
2826  id = (Int_t)t->GetID();
2827  if(id<0) continue;
2828  skipped[nTrksToSkip++] = id;
2829  }
2830  // TEMPORARY FIX
2831  // For AOD, skip also tracks without covariance matrix
2832  if(fInputAOD) {
2833  Double_t covtest[21];
2834  for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
2835  AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
2836  if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
2837  id = (Int_t)vtrack->GetID();
2838  if(id<0) continue;
2839  skipped[nTrksToSkip++] = id;
2840  }
2841  }
2842  }
2843  for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
2844  //
2845  vertexer->SetSkipTracks(nTrksToSkip,skipped);
2846  vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
2847 
2848  } else if(fRmTrksFromPrimVtx && nTrks>0) {
2849  // removing the prongs tracks
2850 
2851  TObjArray rmArray(nTrks);
2852  UShort_t *rmId = new UShort_t[nTrks];
2853  AliESDtrack *esdTrack = 0;
2854  AliESDtrack *t = 0;
2855  for(Int_t i=0; i<nTrks; i++) {
2856  t = (AliESDtrack*)trkArray->UncheckedAt(i);
2857  esdTrack = new AliESDtrack(*t);
2858  rmArray.AddLast(esdTrack);
2859  if(esdTrack->GetID()>=0) {
2860  rmId[i]=(UShort_t)esdTrack->GetID();
2861  } else {
2862  rmId[i]=9999;
2863  }
2864  }
2865  Float_t diamondxy[2]={static_cast<Float_t>(event->GetDiamondX()),static_cast<Float_t>(event->GetDiamondY())};
2866  vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
2867  delete [] rmId; rmId=NULL;
2868  rmArray.Delete();
2869 
2870  }
2871 
2872  if(!vertexESD) return vertexAOD;
2873  if(vertexESD->GetNContributors()<=0) {
2874  //AliDebug(2,"vertexing failed");
2875  delete vertexESD; vertexESD=NULL;
2876  return vertexAOD;
2877  }
2878 
2879  delete vertexer; vertexer=NULL;
2880 
2881  }
2882 
2883  // convert to AliAODVertex
2884  Double_t pos[3],cov[6],chi2perNDF;
2885  vertexESD->GetXYZ(pos); // position
2886  vertexESD->GetCovMatrix(cov); //covariance matrix
2887  chi2perNDF = vertexESD->GetChi2toNDF();
2888  delete vertexESD; vertexESD=NULL;
2889 
2890  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
2891 
2892  return vertexAOD;
2893 }
2894 //-----------------------------------------------------------------------------
2897 
2898  //printf("Preselections:\n");
2899  // fTrackFilter->Dump();
2900  if(fSecVtxWithKF) {
2901  printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
2902  } else {
2903  printf("Secondary vertex with AliVertexerTracks\n");
2904  }
2905  if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
2906  if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
2907  if(fD0toKpi) {
2908  printf("Reconstruct D0->Kpi candidates with cuts:\n");
2910  }
2911  if(fDstar) {
2912  printf("Reconstruct D*->D0pi candidates with cuts:\n");
2913  if(fFindVertexForDstar) {
2914  printf(" Reconstruct a secondary vertex for the D*\n");
2915  } else {
2916  printf(" Assume the D* comes from the primary vertex\n");
2917  }
2919  }
2920  if(fJPSItoEle) {
2921  printf("Reconstruct J/psi from B candidates with cuts:\n");
2923  }
2924  if(f3Prong) {
2925  printf("Reconstruct 3 prong candidates.\n");
2926  printf(" D+->Kpipi cuts:\n");
2928  printf(" Ds->KKpi cuts:\n");
2930  printf(" Lc->pKpi cuts:\n");
2932  }
2933  if(f4Prong) {
2934  printf("Reconstruct 4 prong candidates.\n");
2935  printf(" D0->Kpipipi cuts:\n");
2937  }
2938  if(fCascades) {
2939  printf("Reconstruct cascade candidates formed with v0s.\n");
2940  printf(" Lc -> k0s P & Lc -> L Pi cuts:\n");
2942  printf(" D+ -> K0s pi cuts:\n");
2944  printf(" Ds -> K0s K cuts:\n");
2946  }
2947 
2948  return;
2949 }
2950 //-----------------------------------------------------------------------------
2952  Double_t &dispersion,Bool_t useTRefArray) const
2953 {
2955  //AliCodeTimerAuto("",0);
2956 
2957  AliESDVertex *vertexESD = 0;
2958  AliAODVertex *vertexAOD = 0;
2959 
2960  if(!fSecVtxWithKF) { // AliVertexerTracks
2961 
2962  fVertexerTracks->SetVtxStart(fV1);
2963  vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
2964 
2965  if(!vertexESD) return vertexAOD;
2966 
2967  if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
2968  //AliDebug(2,"vertexing failed");
2969  delete vertexESD; vertexESD=NULL;
2970  return vertexAOD;
2971  }
2972 
2973  Double_t vertRadius2=vertexESD->GetX()*vertexESD->GetX()+vertexESD->GetY()*vertexESD->GetY();
2974  if(vertRadius2>8.){
2975  // vertex outside beam pipe, reject candidate to avoid propagation through material
2976  delete vertexESD; vertexESD=NULL;
2977  return vertexAOD;
2978  }
2979 
2980  } else { // Kalman Filter vertexer (AliKFParticle)
2981 
2982  AliKFParticle::SetField(fBzkG);
2983 
2984  AliKFVertex vertexKF;
2985 
2986  Int_t nTrks = trkArray->GetEntriesFast();
2987  for(Int_t i=0; i<nTrks; i++) {
2988  AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
2989  AliKFParticle daughterKF(*esdTrack,211);
2990  vertexKF.AddDaughter(daughterKF);
2991  }
2992  vertexESD = new AliESDVertex(vertexKF.Parameters(),
2993  vertexKF.CovarianceMatrix(),
2994  vertexKF.GetChi2(),
2995  vertexKF.GetNContributors());
2996 
2997  }
2998 
2999  // convert to AliAODVertex
3000  Double_t pos[3],cov[6],chi2perNDF;
3001  vertexESD->GetXYZ(pos); // position
3002  vertexESD->GetCovMatrix(cov); //covariance matrix
3003  chi2perNDF = vertexESD->GetChi2toNDF();
3004  dispersion = vertexESD->GetDispersion();
3005  delete vertexESD; vertexESD=NULL;
3006 
3007  Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
3008  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
3009 
3010  return vertexAOD;
3011 }
3012 //-----------------------------------------------------------------------------
3015  //AliCodeTimerAuto("",0);
3016 
3017  Int_t retval=kFALSE;
3018  Double_t momentum[3];
3019  Double_t px[3],py[3],pz[3];
3020  for(Int_t iTrack=0; iTrack<3; iTrack++){
3021  AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
3022  track->GetPxPyPz(momentum);
3023  px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
3024  }
3025  retval = SelectInvMassAndPt3prong(px,py,pz);
3026 
3027  return retval;
3028 }
3029 
3030 //-----------------------------------------------------------------------------
3033  //AliCodeTimerAuto("",0);
3034 
3035  Int_t retval=kFALSE;
3036  Double_t momentum[3];
3037  Double_t px[4],py[4],pz[4];
3038 
3039  for(Int_t iTrack=0; iTrack<4; iTrack++){
3040  AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
3041  track->GetPxPyPz(momentum);
3042  px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
3043  }
3044 
3045  retval = SelectInvMassAndPt4prong(px,py,pz);
3046 
3047  return retval;
3048 }
3049 //-----------------------------------------------------------------------------
3052  //AliCodeTimerAuto("",0);
3053 
3054  Int_t retval=kFALSE;
3055  Double_t momentum[3];
3056  Double_t px[2],py[2],pz[2];
3057 
3058  for(Int_t iTrack=0; iTrack<2; iTrack++){
3059  AliESDtrack *track = (AliESDtrack*)trkArray->UncheckedAt(iTrack);
3060  track->GetPxPyPz(momentum);
3061  px[iTrack] = momentum[0]; py[iTrack] = momentum[1]; pz[iTrack] = momentum[2];
3062  }
3063  retval = SelectInvMassAndPtDstarD0pi(px,py,pz);
3064 
3065  return retval;
3066 }
3067 //-----------------------------------------------------------------------------
3069  Double_t *py,
3070  Double_t *pz){
3072  //AliCodeTimerAuto("",0);
3073 
3074  UInt_t pdg2[2];
3075  Int_t nprongs=2;
3076  Double_t minv2,mrange;
3077  Double_t lolim,hilim;
3078  Double_t minPt=0;
3079  Bool_t retval=kFALSE;
3080 
3081  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3082  fOKInvMassD0=kFALSE;
3083  // pt cut
3085  if(minPt>0.1)
3086  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3087  // mass cut
3088  mrange=fCutsD0toKpi->GetMassCut();
3089  lolim=fMassDzero-mrange;
3090  hilim=fMassDzero+mrange;
3091  pdg2[0]=211; pdg2[1]=321;
3092  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3093  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3094  retval=kTRUE;
3095  fOKInvMassD0=kTRUE;
3096  }
3097  pdg2[0]=321; pdg2[1]=211;
3098  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3099  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3100  retval=kTRUE;
3101  fOKInvMassD0=kTRUE;
3102  }
3103  return retval;
3104 }
3105 
3106 //-----------------------------------------------------------------------------
3108  Double_t *py,
3109  Double_t *pz){
3111  //AliCodeTimerAuto("",0);
3112 
3113  UInt_t pdg2[2];
3114  Int_t nprongs=2;
3115  Double_t minv2,mrange;
3116  Double_t lolim,hilim;
3117  Double_t minPt=0;
3118  Bool_t retval=kFALSE;
3119 
3120  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3121  fOKInvMassJpsi=kFALSE;
3122  // pt cut
3124  if(minPt>0.1)
3125  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3126  // mass cut
3127  mrange=fCutsJpsitoee->GetMassCut();
3128  lolim=fMassJpsi-mrange;
3129  hilim=fMassJpsi+mrange;
3130 
3131  pdg2[0]=11; pdg2[1]=11;
3132  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3133  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3134  retval=kTRUE;
3135  fOKInvMassJpsi=kTRUE;
3136  }
3137 
3138  return retval;
3139 }
3140 //-----------------------------------------------------------------------------
3142  Double_t *py,
3143  Double_t *pz,
3144  Int_t pidLcStatus){
3146  //AliCodeTimerAuto("",0);
3147 
3148  UInt_t pdg3[3];
3149  Int_t nprongs=3;
3150  Double_t minv2,mrange;
3151  Double_t lolim,hilim;
3152  Double_t minPt=0;
3153  Bool_t retval=kFALSE;
3154 
3155 
3156  fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
3157  fOKInvMassDplus=kFALSE;
3158  fOKInvMassDs=kFALSE;
3159  fOKInvMassLc=kFALSE;
3160  // pt cut
3162  minPt=TMath::Min(minPt,fCutsLctopKpi->GetMinPtCandidate());
3163  if(minPt>0.1)
3164  if(fMassCalc3->Pt2() < minPt*minPt) return retval;
3165  // D+->Kpipi
3166  mrange=fCutsDplustoKpipi->GetMassCut();
3167  lolim=fMassDplus-mrange;
3168  hilim=fMassDplus+mrange;
3169  pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
3170  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3171  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3172  retval=kTRUE;
3173  fOKInvMassDplus=kTRUE;
3174  }
3175  // Ds+->KKpi
3176  mrange=fCutsDstoKKpi->GetMassCut();
3177  lolim=fMassDs-mrange;
3178  hilim=fMassDs+mrange;
3179  pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
3180  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3181  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3182  retval=kTRUE;
3183  fOKInvMassDs=kTRUE;
3184  }
3185  pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
3186  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3187  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3188  retval=kTRUE;
3189  fOKInvMassDs=kTRUE;
3190  }
3191  // Lc->pKpi
3192  mrange=fCutsLctopKpi->GetMassCut();
3193  lolim=fMassLambdaC-mrange;
3194  hilim=fMassLambdaC+mrange;
3195  if(pidLcStatus&1){
3196  pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
3197  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3198  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3199  retval=kTRUE;
3200  fOKInvMassLc=kTRUE;
3201  }
3202  }
3203  if(pidLcStatus&2){
3204  pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
3205  minv2 = fMassCalc3->InvMass2(nprongs,pdg3);
3206  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3207  retval=kTRUE;
3208  fOKInvMassLc=kTRUE;
3209  }
3210  }
3211 
3212  return retval;
3213 }
3214 
3215 //-----------------------------------------------------------------------------
3217  Double_t *py,
3218  Double_t *pz){
3220  //AliCodeTimerAuto("",0);
3221 
3222  UInt_t pdg2[2];
3223  Int_t nprongs=2;
3224  Double_t minv2,mrange;
3225  Double_t lolim,hilim;
3226  Double_t minPt=0;
3227  Bool_t retval=kFALSE;
3228 
3229  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3230  fOKInvMassDstar=kFALSE;
3231  // pt cut
3233  if(minPt>0.1)
3234  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3235  // mass cut
3236  pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
3237  mrange=fCutsDStartoKpipi->GetMassCut();
3238  lolim=fMassDstar-mrange;
3239  hilim=fMassDstar+mrange;
3240  minv2 = fMassCalc2->InvMass2(nprongs,pdg2);
3241  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3242  retval=kTRUE;
3243  fOKInvMassDstar=kTRUE;
3244  }
3245 
3246  return retval;
3247 }
3248 
3249 //-----------------------------------------------------------------------------
3251  Double_t *py,
3252  Double_t *pz){
3254  //AliCodeTimerAuto("",0);
3255 
3256  UInt_t pdg4[4];
3257  Int_t nprongs=4;
3258  Double_t minv2,mrange;
3259  Double_t lolim,hilim;
3260  Double_t minPt=0;
3261  Bool_t retval=kFALSE;
3262 
3263  // D0->Kpipipi without PID
3264  fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
3265  fOKInvMassD0to4p=kFALSE;
3266  // pt cut
3268  if(minPt>0.1)
3269  if(fMassCalc4->Pt2() < minPt*minPt) return retval;
3270  // mass cut
3271  mrange=fCutsD0toKpipipi->GetMassCut();
3272  lolim=fMassDzero-mrange;
3273  hilim=fMassDzero+mrange;
3274 
3275  pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
3276  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3277  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3278  retval=kTRUE;
3279  fOKInvMassD0to4p=kTRUE;
3280  }
3281 
3282  pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
3283  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3284  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3285  retval=kTRUE;
3286  fOKInvMassD0to4p=kTRUE;
3287  }
3288 
3289  pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
3290  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3291  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3292  retval=kTRUE;
3293  fOKInvMassD0to4p=kTRUE;
3294  }
3295 
3296  pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
3297  minv2 = fMassCalc4->InvMass2(nprongs,pdg4);
3298  if(minv2>lolim*lolim && minv2<hilim*hilim ){
3299  retval=kTRUE;
3300  fOKInvMassD0to4p=kTRUE;
3301  }
3302 
3303  return retval;
3304 }
3305 //-----------------------------------------------------------------------------
3307  Double_t *py,
3308  Double_t *pz){
3310  //AliCodeTimerAuto("",0);
3311 
3312  UInt_t pdg2[2];
3313  Int_t nprongs=2;
3314  Double_t minv2,mrange;
3315  Double_t lolim,hilim;
3316  Double_t minPt=0;
3317  Bool_t retval=kFALSE;
3318 
3319  fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
3320  minPt=fCutsLctoV0->GetMinPtCandidate();
3321  if(minPt>0.1)
3322  if(fMassCalc2->Pt2() < minPt*minPt) return retval;
3323 
3324  fOKInvMassLctoV0 = kFALSE;
3325 
3326  // LambdaC candidate
3327  if (fCutsLctoV0) {
3328  mrange = fCutsLctoV0->GetMassCut();
3329  lolim = fMassLambdaC - mrange;
3330  hilim = fMassLambdaC + mrange;
3331  pdg2[0] = 2212; pdg2[1] = 310;
3332  minv2 = fMassCalc2->InvMass2(2, pdg2);
3333  if ((minv2 > lolim*lolim) && (minv2 < hilim*hilim)) {
3334  retval = kTRUE;
3335  fOKInvMassLctoV0 = kTRUE;
3336  }
3337  pdg2[0] = 211; pdg2[1] = 3122;
3338  minv2 = fMassCalc2->InvMass2(2, pdg2);
3339  if ((minv2>lolim*lolim) && (minv2<hilim*hilim)) {
3340  retval = kTRUE;
3341  fOKInvMassLctoV0 = kTRUE;
3342  }
3343  }
3344 
3345  // D+ cascade candidate
3346  if (fCutsDplustoK0spi) {
3347  mrange = fCutsDplustoK0spi->GetMassCut();
3348  lolim = fMassDplus - mrange;
3349  hilim = fMassDplus + mrange;
3350  pdg2[0] = 211; pdg2[1] = 310;
3351  minv2 = fMassCalc2->InvMass2(2, pdg2);
3352  if ((minv2 > lolim*lolim) && (minv2 < hilim*hilim)) {
3353  retval = kTRUE;
3354  }
3355  }
3356 
3357  // Ds cascade candidate
3358  if (fCutsDstoK0sK) {
3359  mrange = fCutsDstoK0sK->GetMassCut();
3360  lolim = fMassDs - mrange;
3361  hilim = fMassDs + mrange;
3362  pdg2[0] = 321; pdg2[1] = 310;
3363  minv2 = fMassCalc2->InvMass2(2, pdg2);
3364  if ((minv2 > lolim*lolim) && (minv2 < hilim*hilim)) {
3365  retval = kTRUE;
3366  }
3367  }
3368 
3369  return retval;
3370 }
3371 //-----------------------------------------------------------------------------
3373  Int_t trkEntries,
3374  TObjArray &seleTrksArray,
3375  TObjArray &tracksAtVertex,
3376  Int_t &nSeleTrks,
3377  UChar_t *seleFlags,Int_t *evtNumber)
3378 {
3384  //AliCodeTimerAuto("",0);
3385 
3386  const AliVVertex *vprimary = event->GetPrimaryVertex();
3387 
3388  if(fV1) { delete fV1; fV1=NULL; }
3389  if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
3390 
3391  Int_t nindices=0;
3392  UShort_t *indices = 0;
3393  Double_t pos[3],cov[6];
3394  const Int_t entries = event->GetNumberOfTracks();
3395  AliCentrality* cent;
3396 
3397  if(!fInputAOD) { // ESD
3398  fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
3399  cent=((AliESDEvent*)event)->GetCentrality();
3400  } else { // AOD
3401  vprimary->GetXYZ(pos);
3402  vprimary->GetCovarianceMatrix(cov);
3403  fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
3404  if(entries<=0) return;
3405  indices = new UShort_t[entries];
3406  memset(indices,0,sizeof(UShort_t)*entries);
3407  fAODMapSize = 100000;
3408  fAODMap = new Int_t[fAODMapSize];
3409  memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
3410  cent=((AliAODEvent*)event)->GetCentrality();
3411  }
3412  Float_t centperc=0.1;
3413  if(event->GetRunNumber()<244824){
3414  centperc=cent->GetCentralityPercentile("V0M");
3415  }else{
3416  AliMultSelection *multSelection = (AliMultSelection * ) event->FindListObject("MultSelection");
3417  if(multSelection){
3418  centperc=multSelection->GetMultiplicityPercentile("V0M");
3419  Int_t qual = multSelection->GetEvSelCode();
3420  if(qual == 199 ) centperc=0.1; // use central cuts as default
3421  }
3422  }
3423  Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE,okFor3Prong=kFALSE,okBachelor=kFALSE;
3424  nSeleTrks=0;
3425 
3426  // transfer ITS tracks from event to arrays
3427  for(Int_t i=0; i<entries; i++) {
3428  AliVTrack *track;
3429  track = (AliVTrack*)event->GetTrack(i);
3430 
3431  // skip pure ITS SA tracks
3432  if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
3433 
3434  // skip tracks without ITS
3435  if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
3436 
3437  // skip tracks with negative ID
3438  // (these are duplicated TPC-only AOD tracks, for jet analysis...)
3439  if(track->GetID()<0) continue;
3440 
3441  // TEMPORARY: check that the cov matrix is there
3442  Double_t covtest[21];
3443  if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
3444  //
3445 
3446  if(fInputAOD) {
3447  AliAODTrack *aodt = (AliAODTrack*)track;
3448  if(aodt->GetUsedForPrimVtxFit()) {
3449  indices[nindices]=aodt->GetID(); nindices++;
3450  }
3451  Int_t ind = (Int_t)aodt->GetID();
3452  if (ind>-1 && ind < fAODMapSize) fAODMap[ind] = i;
3453  }
3454 
3455  AliESDtrack *esdt = 0;
3456 
3457  if(!fInputAOD) {
3458  esdt = (AliESDtrack*)track;
3459  } else {
3460  esdt = new AliESDtrack(track);
3461  }
3462 
3463  // single track selection
3464  okDisplaced=kFALSE; okSoftPi=kFALSE; okFor3Prong=kFALSE; okBachelor=kFALSE;
3465  if(fMixEvent && i<trkEntries){
3466  evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
3467  const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
3468  Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
3469  eventVtx->GetXYZ(vtxPos);
3470  vprimary->GetXYZ(primPos);
3471  eventVtx->GetCovarianceMatrix(primCov);
3472  for(Int_t ind=0;ind<3;ind++){
3473  trasl[ind]=vtxPos[ind]-primPos[ind];
3474  }
3475 
3476  Bool_t isTransl=esdt->Translate(trasl,primCov);
3477  if(!isTransl) {
3478  delete esdt;
3479  esdt = NULL;
3480  continue;
3481  }
3482  }
3483 
3484  if(SingleTrkCuts(esdt,centperc,okDisplaced,okSoftPi,okFor3Prong,okBachelor) && nSeleTrks<trkEntries) {
3485  esdt->PropagateToDCA(fV1,fBzkG,kVeryBig);
3486  seleTrksArray.AddLast(esdt);
3487  tracksAtVertex.AddLast(new AliExternalTrackParam(*esdt));
3488  seleFlags[nSeleTrks]=0;
3489  if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
3490  if(okFor3Prong) SETBIT(seleFlags[nSeleTrks],kBit3Prong);
3491  if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
3492  if(okBachelor) SETBIT(seleFlags[nSeleTrks],kBitBachelor);
3493  // Check the PID
3494  SETBIT(seleFlags[nSeleTrks],kBitPionCompat);
3495  SETBIT(seleFlags[nSeleTrks],kBitKaonCompat);
3496  SETBIT(seleFlags[nSeleTrks],kBitProtonCompat);
3497  Bool_t useTPC=kTRUE;
3498  if(fUseTOFPID){
3499  Double_t nsigmatofPi= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kPion);
3500  if(nsigmatofPi>-990. && (nsigmatofPi<-fnSigmaTOFPionLow || nsigmatofPi>fnSigmaTOFPionHi)){
3501  CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
3502  }
3503  Double_t nsigmatofK= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kKaon);
3504  if(nsigmatofK>-990. && (nsigmatofK<-fnSigmaTOFKaonLow || nsigmatofK>fnSigmaTOFKaonHi)){
3505  CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
3506  }
3507  Double_t nsigmatofP= fPidResponse->NumberOfSigmasTOF(esdt,AliPID::kProton);
3508  if(nsigmatofP>-990. && (nsigmatofP<-fnSigmaTOFProtonLow || nsigmatofP>fnSigmaTOFProtonHi)){
3509  CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
3510  }
3511  if(fUseTPCPIDOnlyIfNoTOF && nsigmatofPi>-990.) useTPC=kFALSE;
3512  }
3513  if(useTPC && fUseTPCPID && esdt->P()<fMaxMomForTPCPid){
3514  Double_t nsigmatpcPi= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kPion);
3515  if(nsigmatpcPi>-990. && (nsigmatpcPi<-fnSigmaTPCPionLow || nsigmatpcPi>fnSigmaTPCPionHi)){
3516  CLRBIT(seleFlags[nSeleTrks],kBitPionCompat);
3517  }
3518  Double_t nsigmatpcK= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kKaon);
3519  if(nsigmatpcK>-990. && (nsigmatpcK<-fnSigmaTPCKaonLow || nsigmatpcK>fnSigmaTPCKaonHi)){
3520  CLRBIT(seleFlags[nSeleTrks],kBitKaonCompat);
3521  }
3522  Double_t nsigmatpcP= fPidResponse->NumberOfSigmasTPC(esdt,AliPID::kProton);
3523  if(nsigmatpcP>-990. && (nsigmatpcP<-fnSigmaTPCProtonLow || nsigmatpcP>fnSigmaTPCProtonHi)){
3524  CLRBIT(seleFlags[nSeleTrks],kBitProtonCompat);
3525  }
3526  }
3527  nSeleTrks++;
3528  } else {
3529  if(fInputAOD) delete esdt;
3530  esdt = NULL;
3531  continue;
3532  }
3533 
3534  } // end loop on tracks
3535 
3536  // primary vertex from AOD
3537  if(fInputAOD) {
3538  delete fV1; fV1=NULL;
3539  vprimary->GetXYZ(pos);
3540  vprimary->GetCovarianceMatrix(cov);
3541  Double_t chi2toNDF = vprimary->GetChi2perNDF();
3542  Int_t ncontr=nindices;
3543  if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
3544  Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
3545  fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
3546  fV1->SetTitle(vprimary->GetTitle());
3547  fV1->SetIndices(nindices,indices);
3548  }
3549  if(indices) { delete [] indices; indices=NULL; }
3550 
3551 
3552  return;
3553 }
3554 //-----------------------------------------------------------------------------
3556  //
3558  //
3559  if(fUsePidTag && cuts->GetPidHF()) {
3560  Bool_t usepid=cuts->GetIsUsePID();
3561  cuts->SetUsePID(kTRUE);
3562  if(cuts->IsSelectedPID(rd))
3563  rd->SetSelectionBit(bit);
3564  cuts->SetUsePID(usepid);
3565  }
3566  return;
3567 }
3568 //-----------------------------------------------------------------------------
3570  Float_t centralityperc,
3571  Bool_t &okDisplaced,
3572  Bool_t &okSoftPi,
3573  Bool_t &okFor3Prong,
3574  Bool_t &okBachelor) const
3575 {
3577 
3578  // this is needed to store the impact parameters
3579  //AliCodeTimerAuto("",0);
3580 
3581  if (!trk->PropagateToDCA(fV1,fBzkG,kVeryBig)) return kFALSE;
3582 
3583  trk->RelateToVertex(fV1,fBzkG,kVeryBig);
3584 
3585  UInt_t selectInfo;
3586  //
3587  // Track selection, displaced tracks -- 2 prong
3588  selectInfo = 0;
3589  if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
3591  // central PbPb events, tighter cuts
3592  selectInfo = fTrackFilter2prongCentral->IsSelected(trk);
3593  }else{
3594  // standard cuts
3595  if(fTrackFilter) {
3596  selectInfo = fTrackFilter->IsSelected(trk);
3597  }
3598  }
3599  if(selectInfo) okDisplaced=kTRUE;
3600 
3601  // Track selection, displaced tracks -- 3 prong
3602  selectInfo = 0;
3603  if(centralityperc>=0 && fMaxCentPercentileForTightCuts>=0
3605  // central PbPb events, tighter cuts
3606  selectInfo = fTrackFilter3prongCentral->IsSelected(trk);
3607  }else{
3608  // standard cuts
3609  if(fTrackFilter) {
3610  selectInfo = fTrackFilter->IsSelected(trk);
3611  }
3612  }
3613  if(selectInfo) okFor3Prong=kTRUE;
3614 
3615  // Track selection, soft pions
3616  selectInfo = 0;
3617  if(fDstar && fTrackFilterSoftPi) {
3618  selectInfo = fTrackFilterSoftPi->IsSelected(trk);
3619  }
3620  if(selectInfo) okSoftPi=kTRUE;
3621 
3622  // Track selection, bachelor
3623  selectInfo = 0;
3624  if(fCascades){
3626  selectInfo = fTrackFilterBachelor->IsSelected(trk);
3627  }else{
3628  if(okDisplaced) selectInfo=1;
3629  }
3630  }
3631  if(selectInfo) okBachelor=kTRUE;
3632 
3633  if(okDisplaced || okSoftPi || okFor3Prong || okBachelor) return kTRUE;
3634 
3635  return kFALSE;
3636 }
3637 
3638 
3639 //-----------------------------------------------------------------------------
3640 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
3646  //
3647  //AliCodeTimerAuto("",0);
3648  Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
3649  AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
3650 
3651  // create the v0 neutral track to compute the DCA to the primary vertex
3652  Double_t xyz[3], pxpypz[3];
3653  esdV0->XvYvZv(xyz);
3654  esdV0->PxPyPz(pxpypz);
3655  Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
3656  AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
3657  if(!trackesdV0) {
3658  delete vertexV0;
3659  return 0;
3660  }
3661  Double_t d0z0[2],covd0z0[3];
3662  AliAODVertex *primVertexAOD = PrimaryVertex();
3663  trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
3664  Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
3665  // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
3666  Double_t dcaV0DaughterToPrimVertex[2];
3667  AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
3668  AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
3669  if( !posV0track || !negV0track) {
3670  if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
3671  delete vertexV0;
3672  delete primVertexAOD;
3673  return 0;
3674  }
3675  posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
3676  // if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
3677  // else
3678  dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
3679  negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
3680  // if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
3681  // else
3682  dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
3683  Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
3684  Double_t pmom[3],nmom[3];
3685  esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
3686  esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
3687 
3688  AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
3689  aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
3690 
3691  delete trackesdV0;
3692  delete primVertexAOD;
3693 
3694  return aodV0;
3695 }
3696 //-----------------------------------------------------------------------------
3697 void AliAnalysisVertexingHF::SetParametersAtVertex(AliESDtrack* esdt, const AliExternalTrackParam* extpar) const{
3699  //AliCodeTimerAuto("",0);
3700 
3701  const Double_t *par=extpar->GetParameter();
3702  const Double_t *cov=extpar->GetCovariance();
3703  Double_t alpha=extpar->GetAlpha();
3704  Double_t x=extpar->GetX();
3705  esdt->Set(x,alpha,par,cov);
3706  return;
3707 }
3708 //-----------------------------------------------------------------------------
3711 
3712  fMassDzero=TDatabasePDG::Instance()->GetParticle(421)->Mass();
3713  fMassDplus=TDatabasePDG::Instance()->GetParticle(411)->Mass();
3714  fMassDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
3715  fMassLambdaC=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
3716  fMassDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
3717  fMassJpsi=TDatabasePDG::Instance()->GetParticle(443)->Mass();
3718 }
3719 //-----------------------------------------------------------------------------
3721  //
3723  //
3724 
3725 
3726  //___ Check if the V0 type from AliRDHFCutsLctoV0 is the same as the one set in the ConfigVertexingHF.C for AliAnalysisVertexingHF
3727 
3728 
3730  printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHFCutsLctoV0 (%d)\n",fV0TypeForCascadeVertex,fCutsLctoV0->GetV0Type());
3731  return kFALSE;
3732  }
3733 
3735  printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHfCutsDplustoK0spi (%d)\n",fV0TypeForCascadeVertex,fCutsDplustoK0spi->GetV0Type());
3736  return kFALSE;
3737  }
3738 
3740  printf("ERROR: V0 type doesn not match in AliAnalysisVertexingHF (%d) required in AliRDHFCutsDstoK0sK (%d)\n",fV0TypeForCascadeVertex,fCutsDstoK0sK->GetV0Type());
3741  return kFALSE;
3742  }
3743 
3744  return kTRUE;
3745 }
3746 //-----------------------------------------------------------------------------
3747 
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:315
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:261
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:219
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)
virtual void PrintAll() const
Bool_t GetIsUsePID() const
Definition: AliRDHFCuts.h:284
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:298
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.