1 /**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //----------------------------------------------------------------------------
17 // Implementation of the heavy-flavour vertexing analysis class
18 // Candidates are stored in the AOD as objects deriving from AliAODRecoDecay.
19 // To be used as a task of AliAnalysisManager by means of the interface
20 // class AliAnalysisTaskSEVertexingHF.
21 // An example of usage in the macro AliAnalysisTaskSEVertexingHFTest.C.
23 // Contact: andrea.dainese@pd.infn.it
24 // Contributors: E.Bruna, G.E.Bruno, A.Dainese, C.Di Gliglio,
25 // F.Prino, R.Romita, X.M.Zhang
26 //----------------------------------------------------------------------------
27 #include <Riostream.h>
29 #include <TDatabasePDG.h>
33 #include "AliVEvent.h"
34 #include "AliVVertex.h"
35 #include "AliVTrack.h"
36 #include "AliVertexerTracks.h"
37 #include "AliKFVertex.h"
38 #include "AliESDEvent.h"
39 #include "AliESDVertex.h"
40 #include "AliExternalTrackParam.h"
41 #include "AliNeutralTrackParam.h"
42 #include "AliESDtrack.h"
43 #include "AliESDtrackCuts.h"
44 #include "AliAODEvent.h"
45 #include "AliAODRecoDecay.h"
46 #include "AliAODRecoDecayHF.h"
47 #include "AliAODRecoDecayHF2Prong.h"
48 #include "AliAODRecoDecayHF3Prong.h"
49 #include "AliAODRecoDecayHF4Prong.h"
50 #include "AliAODRecoCascadeHF.h"
51 #include "AliRDHFCutsD0toKpi.h"
52 #include "AliRDHFCutsJpsitoee.h"
53 #include "AliRDHFCutsDplustoKpipi.h"
54 #include "AliRDHFCutsDstoKKpi.h"
55 #include "AliRDHFCutsLctopKpi.h"
56 #include "AliRDHFCutsLctoV0.h"
57 #include "AliRDHFCutsD0toKpipipi.h"
58 #include "AliRDHFCutsDStartoKpipi.h"
59 #include "AliAnalysisFilter.h"
60 #include "AliAnalysisVertexingHF.h"
61 #include "AliMixedEvent.h"
65 ClassImp(AliAnalysisVertexingHF)
67 //----------------------------------------------------------------------------
68 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
73 fSecVtxWithKF(kFALSE),
74 fRecoPrimVtxSkippingTrks(kFALSE),
75 fRmTrksFromPrimVtx(kFALSE),
86 fTrackFilterSoftPi(0x0),
89 fCutsDplustoKpipi(0x0),
93 fCutsD0toKpipipi(0x0),
94 fCutsDStartoKpipi(0x0),
96 fFindVertexForDstar(kTRUE),
97 fFindVertexForCascades(kTRUE)
99 // Default constructor
107 SetD0to4ProngsCuts();
110 //--------------------------------------------------------------------------
111 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
113 fInputAOD(source.fInputAOD),
114 fAODMapSize(source.fAODMapSize),
115 fAODMap(source.fAODMap),
117 fSecVtxWithKF(source.fSecVtxWithKF),
118 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
119 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
121 fD0toKpi(source.fD0toKpi),
122 fJPSItoEle(source.fJPSItoEle),
123 f3Prong(source.f3Prong),
124 f4Prong(source.f4Prong),
125 fDstar(source.fDstar),
126 fCascades(source.fCascades),
127 fLikeSign(source.fLikeSign),
128 fMixEvent(source.fMixEvent),
129 fTrackFilter(source.fTrackFilter),
130 fTrackFilterSoftPi(source.fTrackFilterSoftPi),
131 fCutsD0toKpi(source.fCutsD0toKpi),
132 fCutsJpsitoee(source.fCutsJpsitoee),
133 fCutsDplustoKpipi(source.fCutsDplustoKpipi),
134 fCutsDstoKKpi(source.fCutsDstoKKpi),
135 fCutsLctopKpi(source.fCutsLctopKpi),
136 fCutsLctoV0(source.fCutsLctoV0),
137 fCutsD0toKpipipi(source.fCutsD0toKpipipi),
138 fCutsDStartoKpipi(source.fCutsDStartoKpipi),
139 fListOfCuts(source.fListOfCuts),
140 fFindVertexForDstar(source.fFindVertexForDstar),
141 fFindVertexForCascades(source.fFindVertexForCascades)
146 for(Int_t i=0; i<9; i++) fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
147 for(Int_t i=0; i<9; i++) fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
148 for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
149 for(Int_t i=0; i<14; i++) fDsCuts[i]=source.fDsCuts[i];
150 for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
151 for(Int_t i=0; i<8; i++) fLctoV0Cuts[i]=source.fLctoV0Cuts[i];
152 for(Int_t i=0; i<5; i++) fDstarCuts[i]=source.fDstarCuts[i];
153 for(Int_t i=0; i<9; i++) fD0to4ProngsCuts[i]=source.fD0to4ProngsCuts[i];
155 //--------------------------------------------------------------------------
156 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
159 // assignment operator
161 if(&source == this) return *this;
162 fInputAOD = source.fInputAOD;
163 fAODMapSize = source.fAODMapSize;
164 fBzkG = source.fBzkG;
165 fSecVtxWithKF = source.fSecVtxWithKF;
166 fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
167 fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
169 fD0toKpi = source.fD0toKpi;
170 fJPSItoEle = source.fJPSItoEle;
171 f3Prong = source.f3Prong;
172 f4Prong = source.f4Prong;
173 fDstar = source.fDstar;
174 fCascades = source.fCascades;
175 fLikeSign = source.fLikeSign;
176 fMixEvent = source.fMixEvent;
177 fTrackFilter = source.fTrackFilter;
178 fTrackFilterSoftPi = source.fTrackFilterSoftPi;
179 fCutsD0toKpi = source.fCutsD0toKpi;
180 fCutsJpsitoee = source.fCutsJpsitoee;
181 fCutsDplustoKpipi = source.fCutsDplustoKpipi;
182 fCutsDstoKKpi = source.fCutsDstoKKpi;
183 fCutsLctopKpi = source.fCutsLctopKpi;
184 fCutsLctoV0 = source.fCutsLctoV0;
185 fCutsD0toKpipipi = source.fCutsD0toKpipipi;
186 fCutsDStartoKpipi = source.fCutsDStartoKpipi;
187 fListOfCuts = source.fListOfCuts;
188 fFindVertexForDstar = source.fFindVertexForDstar;
189 fFindVertexForCascades = source.fFindVertexForCascades;
191 for(Int_t i=0; i<9; i++) fD0toKpiCuts[i]=source.fD0toKpiCuts[i];
192 for(Int_t i=0; i<9; i++) fBtoJPSICuts[i]=source.fBtoJPSICuts[i];
193 for(Int_t i=0; i<12; i++) fDplusCuts[i]=source.fDplusCuts[i];
194 for(Int_t i=0; i<14; i++) fDsCuts[i]=source.fDsCuts[i];
195 for(Int_t i=0; i<12; i++) fLcCuts[i]=source.fLcCuts[i];
196 for(Int_t i=0; i<8; i++) fLctoV0Cuts[i]=source.fLctoV0Cuts[i];
197 for(Int_t i=0; i<5; i++) fDstarCuts[i]=source.fDstarCuts[i];
198 for(Int_t i=0; i<9; i++) fD0to4ProngsCuts[i]=source.fD0to4ProngsCuts[i];
202 //----------------------------------------------------------------------------
203 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
205 if(fV1) { delete fV1; fV1=0; }
206 if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
207 if(fTrackFilterSoftPi) { delete fTrackFilterSoftPi; fTrackFilterSoftPi=0; }
208 if(fCutsD0toKpi) { delete fCutsD0toKpi; fCutsD0toKpi=0; }
209 if(fCutsJpsitoee) { delete fCutsJpsitoee; fCutsJpsitoee=0; }
210 if(fCutsDplustoKpipi) { delete fCutsDplustoKpipi; fCutsDplustoKpipi=0; }
211 if(fCutsDstoKKpi) { delete fCutsDstoKKpi; fCutsDstoKKpi=0; }
212 if(fCutsLctopKpi) { delete fCutsLctopKpi; fCutsLctopKpi=0; }
213 if(fCutsLctoV0) { delete fCutsLctoV0; fCutsLctoV0=0; }
214 if(fCutsD0toKpipipi) { delete fCutsD0toKpipipi; fCutsD0toKpipipi=0; }
215 if(fCutsDStartoKpipi) { delete fCutsDStartoKpipi; fCutsDStartoKpipi=0; }
216 if(fAODMap) { delete fAODMap; fAODMap=0; }
218 //----------------------------------------------------------------------------
219 TList *AliAnalysisVertexingHF::FillListOfCuts() {
220 // Fill list of analysis cuts
222 TList *list = new TList();
224 list->SetName("ListOfCuts");
227 AliRDHFCutsD0toKpi *cutsD0toKpi = new AliRDHFCutsD0toKpi(*fCutsD0toKpi);
228 list->Add(cutsD0toKpi);
231 AliRDHFCutsJpsitoee *cutsJpsitoee = new AliRDHFCutsJpsitoee(*fCutsJpsitoee);
232 list->Add(cutsJpsitoee);
234 if(fCutsDplustoKpipi) {
235 AliRDHFCutsDplustoKpipi *cutsDplustoKpipi = new AliRDHFCutsDplustoKpipi(*fCutsDplustoKpipi);
236 list->Add(cutsDplustoKpipi);
239 AliRDHFCutsDstoKKpi *cutsDstoKKpi = new AliRDHFCutsDstoKKpi(*fCutsDstoKKpi);
240 list->Add(cutsDstoKKpi);
243 AliRDHFCutsLctopKpi *cutsLctopKpi = new AliRDHFCutsLctopKpi(*fCutsLctopKpi);
244 list->Add(cutsLctopKpi);
247 AliRDHFCutsLctoV0 *cutsLctoV0 = new AliRDHFCutsLctoV0(*fCutsLctoV0);
248 list->Add(cutsLctoV0);
250 if(fCutsD0toKpipipi) {
251 AliRDHFCutsD0toKpipipi *cutsD0toKpipipi = new AliRDHFCutsD0toKpipipi(*fCutsD0toKpipipi);
252 list->Add(cutsD0toKpipipi);
254 if(fCutsDStartoKpipi) {
255 AliRDHFCutsDStartoKpipi *cutsDStartoKpipi = new AliRDHFCutsDStartoKpipi(*fCutsDStartoKpipi);
256 list->Add(cutsDStartoKpipi);
259 // keep a pointer to the list
264 //----------------------------------------------------------------------------
265 void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
266 TClonesArray *aodVerticesHFTClArr,
267 TClonesArray *aodD0toKpiTClArr,
268 TClonesArray *aodJPSItoEleTClArr,
269 TClonesArray *aodCharm3ProngTClArr,
270 TClonesArray *aodCharm4ProngTClArr,
271 TClonesArray *aodDstarTClArr,
272 TClonesArray *aodCascadesTClArr,
273 TClonesArray *aodLikeSign2ProngTClArr,
274 TClonesArray *aodLikeSign3ProngTClArr)
276 // Find heavy-flavour vertex candidates
278 // Output: AOD (additional branches added)
281 TString evtype = event->IsA()->GetName();
282 fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
283 } // if we do mixing AliVEvent is a AliMixedEvent
286 AliDebug(2,"Creating HF candidates from AOD");
288 AliDebug(2,"Creating HF candidates from ESD");
291 if(!aodVerticesHFTClArr) {
292 printf("ERROR: no aodVerticesHFTClArr");
295 if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
296 printf("ERROR: no aodD0toKpiTClArr");
299 if(fJPSItoEle && !aodJPSItoEleTClArr) {
300 printf("ERROR: no aodJPSItoEleTClArr");
303 if(f3Prong && !aodCharm3ProngTClArr) {
304 printf("ERROR: no aodCharm3ProngTClArr");
307 if(f4Prong && !aodCharm4ProngTClArr) {
308 printf("ERROR: no aodCharm4ProngTClArr");
311 if(fDstar && !aodDstarTClArr) {
312 printf("ERROR: no aodDstarTClArr");
315 if(fCascades && !aodCascadesTClArr){
316 printf("ERROR: no aodCascadesTClArr ");
319 if(fLikeSign && !aodLikeSign2ProngTClArr) {
320 printf("ERROR: no aodLikeSign2ProngTClArr");
323 if(fLikeSign && f3Prong && !aodLikeSign3ProngTClArr) {
324 printf("ERROR: no aodLikeSign2ProngTClArr");
328 // delete candidates from previous event and create references
329 Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0,iCascades=0,iLikeSign2Prong=0,iLikeSign3Prong=0;
330 aodVerticesHFTClArr->Delete();
331 iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
332 TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
333 if(fD0toKpi || fDstar) {
334 aodD0toKpiTClArr->Delete();
335 iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
338 aodJPSItoEleTClArr->Delete();
339 iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
342 aodCharm3ProngTClArr->Delete();
343 i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
346 aodCharm4ProngTClArr->Delete();
347 i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
350 aodDstarTClArr->Delete();
351 iDstar = aodDstarTClArr->GetEntriesFast();
354 aodCascadesTClArr->Delete();
355 iCascades = aodCascadesTClArr->GetEntriesFast();
358 aodLikeSign2ProngTClArr->Delete();
359 iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast();
361 if(fLikeSign && f3Prong) {
362 aodLikeSign3ProngTClArr->Delete();
363 iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast();
366 TClonesArray &aodD0toKpiRef = *aodD0toKpiTClArr;
367 TClonesArray &aodJPSItoEleRef = *aodJPSItoEleTClArr;
368 TClonesArray &aodCharm3ProngRef = *aodCharm3ProngTClArr;
369 TClonesArray &aodCharm4ProngRef = *aodCharm4ProngTClArr;
370 TClonesArray &aodDstarRef = *aodDstarTClArr;
371 TClonesArray &aodCascadesRef = *aodCascadesTClArr;
372 TClonesArray &aodLikeSign2ProngRef = *aodLikeSign2ProngTClArr;
373 TClonesArray &aodLikeSign3ProngRef = *aodLikeSign3ProngTClArr;
376 AliAODRecoDecayHF2Prong *io2Prong = 0;
377 AliAODRecoDecayHF3Prong *io3Prong = 0;
378 AliAODRecoDecayHF4Prong *io4Prong = 0;
379 AliAODRecoCascadeHF *ioCascade = 0;
381 Int_t iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries,iv0,nv0;
382 Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaV0,dcaCasc;
383 Bool_t okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
384 Bool_t okDstar=kFALSE,okD0fromDstar=kFALSE;
385 Bool_t okCascades=kFALSE;
386 AliESDtrack *postrack1 = 0;
387 AliESDtrack *postrack2 = 0;
388 AliESDtrack *negtrack1 = 0;
389 AliESDtrack *negtrack2 = 0;
390 AliESDtrack *trackPi = 0;
391 // AliESDtrack *posV0track = 0;
392 // AliESDtrack *negV0track = 0;
394 Double_t dcaMax = fD0toKpiCuts[1];
395 if(dcaMax < fBtoJPSICuts[1]) dcaMax=fBtoJPSICuts[1];
396 if(dcaMax < fDplusCuts[11]) dcaMax=fDplusCuts[11];
397 if(dcaMax < fD0to4ProngsCuts[1]) dcaMax=fD0to4ProngsCuts[1];
399 Float_t dcaMax = fCutsD0toKpi->GetDCACut();
400 if(fCutsJpsitoee) dcaMax=TMath::Max(dcaMax,fCutsJpsitoee->GetDCACut());
401 if(fCutsDplustoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDplustoKpipi->GetDCACut());
402 if(fCutsDstoKKpi) dcaMax=TMath::Max(dcaMax,fCutsDstoKKpi->GetDCACut());
403 if(fCutsLctopKpi) dcaMax=TMath::Max(dcaMax,fCutsLctopKpi->GetDCACut());
404 if(fCutsD0toKpipipi) dcaMax=TMath::Max(dcaMax,fCutsD0toKpipipi->GetDCACut());
405 if(fCutsDStartoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDStartoKpipi->GetDCACut());
407 AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
411 fBzkG = (Double_t)event->GetMagneticField();
413 trkEntries = (Int_t)event->GetNumberOfTracks();
414 AliDebug(1,Form(" Number of tracks: %d",trkEntries));
416 nv0 = (Int_t)event->GetNumberOfV0s();
417 AliDebug(1,Form(" Number of V0s: %d",nv0));
419 if( trkEntries<2 && (trkEntries<1 || nv0<1) ) {
420 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
425 if(!fCutsD0toKpi->IsEventSelected(event)) return;
427 // call function that applies sigle-track selection,
428 // for displaced tracks and soft pions (both charges) for D*,
429 // and retrieves primary vertex
430 TObjArray seleTrksArray(trkEntries);
431 UChar_t *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi
433 Int_t *evtNumber = new Int_t[trkEntries];
434 SelectTracksAndCopyVertex(event,seleTrksArray,nSeleTrks,seleFlags,evtNumber);
436 AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
438 TObjArray *twoTrackArray1 = new TObjArray(2);
439 TObjArray *twoTrackArray2 = new TObjArray(2);
440 TObjArray *twoTrackArrayV0 = new TObjArray(2);
441 TObjArray *twoTrackArrayCasc = new TObjArray(2);
442 TObjArray *threeTrackArray = new TObjArray(3);
443 TObjArray *fourTrackArray = new TObjArray(4);
446 Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
448 AliAODRecoDecayHF *rd = 0;
449 AliAODRecoCascadeHF *rc = 0;
453 // LOOP ON POSITIVE TRACKS
454 for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
456 if(iTrkP1%1==0) AliDebug(1,Form(" 1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));
458 // get track from tracks array
459 postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
461 if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
463 // Make cascades with V0+track
467 for(iv0=0; iv0<nv0; iv0++){
469 AliDebug(1,Form(" loop on v0s for track number %d and v0 number %d",iTrkP1,iv0));
473 v0 = ((AliAODEvent*)event)->GetV0(iv0);
475 esdV0 = ((AliESDEvent*)event)->GetV0(iv0);
477 if ( (!v0 || !v0->IsA()->InheritsFrom("AliAODv0") ) &&
478 (!esdV0 || !esdV0->IsA()->InheritsFrom("AliESDv0") ) ) continue;
481 // Get the tracks that form the V0
482 // ( parameters at primary vertex )
483 // and define an AliExternalTrackParam out of them
484 AliExternalTrackParam * posV0track;
485 AliExternalTrackParam * negV0track;
488 AliAODTrack *posVV0track = (AliAODTrack*)(v0->GetDaughter(0));
489 AliAODTrack *negVV0track = (AliAODTrack*)(v0->GetDaughter(1));
490 if( !posVV0track || !negVV0track ) continue;
492 // Apply some basic V0 daughter criteria
494 // bachelor must not be a v0-track
495 if (posVV0track->GetID() == postrack1->GetID() ||
496 negVV0track->GetID() == postrack1->GetID()) continue;
497 // reject like-sign v0
498 if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
499 // avoid ghost TPC tracks
500 if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
501 !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
502 // Get AliExternalTrackParam out of the AliAODTracks
503 Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign;
504 posVV0track->PxPyPz(pxpypz); posVV0track->XvYvZv(xyz);
505 posVV0track->GetCovarianceXYZPxPyPz(cv); sign=posVV0track->Charge();
506 posV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
507 negVV0track->PxPyPz(pxpypz); negVV0track->XvYvZv(xyz);
508 negVV0track->GetCovarianceXYZPxPyPz(cv); sign=negVV0track->Charge();
509 negV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
511 AliESDtrack *posVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetPindex() ));
512 AliESDtrack *negVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetNindex() ));
513 if( !posVV0track || !negVV0track ) continue;
515 // Apply some basic V0 daughter criteria
517 // bachelor must not be a v0-track
518 if (posVV0track->GetID() == postrack1->GetID() ||
519 negVV0track->GetID() == postrack1->GetID()) continue;
520 // reject like-sign v0
521 if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
522 // avoid ghost TPC tracks
523 if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
524 !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
525 // reject kinks (only necessary on AliESDtracks)
526 if (posVV0track->GetKinkIndex(0)>0 || negVV0track->GetKinkIndex(0)>0) continue;
527 // Get AliExternalTrackParam out of the AliESDtracks
528 posV0track = new AliExternalTrackParam(*posVV0track);
529 negV0track = new AliExternalTrackParam(*negVV0track);
531 // Define the AODv0 from ESDv0 if reading ESDs
532 v0 = TransformESDv0toAODv0(esdV0,twoTrackArrayV0);
534 if( !posV0track || !negV0track ){
535 AliDebug(1,Form(" Couldn't get the V0 daughters"));
539 // fill in the v0 two-external-track-param array
540 twoTrackArrayV0->AddAt(posV0track,0);
541 twoTrackArrayV0->AddAt(negV0track,1);
544 dcaV0 = v0->DcaV0Daughters();
546 // Define the V0 (neutral) track
547 AliNeutralTrackParam *trackV0;
549 const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
550 if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
552 Double_t xyz[3], pxpypz[3];
554 esdV0->PxPyPz(pxpypz);
555 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
556 trackV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
558 // Fill in the object array to create the cascade
559 twoTrackArrayCasc->AddAt(postrack1,0);
560 twoTrackArrayCasc->AddAt(trackV0,1);
562 // Compute the cascade vertex
563 AliAODVertex *vertexCasc = 0;
564 if(fFindVertexForCascades) {
565 // DCA between the two tracks
566 dcaCasc = postrack1->GetDCA(trackV0,fBzkG,xdummy,ydummy);
568 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
570 // assume Cascade decays at the primary vertex
571 Double_t pos[3],cov[6],chi2perNDF;
573 fV1->GetCovMatrix(cov);
574 chi2perNDF = fV1->GetChi2toNDF();
575 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
579 delete posV0track; posV0track=NULL;
580 delete negV0track; negV0track=NULL;
581 delete trackV0; trackV0=NULL;
582 if(!fInputAOD) {delete v0; v0=NULL;}
583 twoTrackArrayCasc->Clear();
587 // Create and store the Cascade if passed the cuts
588 ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,v0,dcaCasc,okCascades);
589 if(okCascades && ioCascade) {
590 AliDebug(1,Form("Storing a cascade object... "));
591 // add the vertex and the cascade to the AOD
592 AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
593 rc = new(aodCascadesRef[iCascades++])AliAODRecoCascadeHF(*ioCascade);
594 rc->SetSecondaryVtx(vCasc);
595 vCasc->SetParent(rc);
596 rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
597 if(!fInputAOD) vCasc->AddDaughter(v0); // just to fill ref #0 ??
598 AddRefs(vCasc,rc,event,twoTrackArrayCasc); // add the track (proton)
599 vCasc->AddDaughter(v0); // fill the 2prong V0
603 delete posV0track; posV0track=NULL;
604 delete negV0track; negV0track=NULL;
605 delete trackV0; trackV0=NULL;
606 twoTrackArrayV0->Clear();
607 twoTrackArrayCasc->Clear();
608 if(ioCascade) { delete ioCascade; ioCascade=NULL; }
609 if(vertexCasc) { delete vertexCasc; vertexCasc=NULL; }
610 if(!fInputAOD) {delete v0; v0=NULL;}
612 } // end loop on V0's
615 // If there is less than 2 particles exit
617 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
621 if(postrack1->Charge()<0 && !fLikeSign) continue;
623 // LOOP ON NEGATIVE TRACKS
624 for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
626 if(iTrkN1%1==0) AliDebug(1,Form(" 1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));
628 if(iTrkN1==iTrkP1) continue;
630 // get track from tracks array
631 negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
633 if(negtrack1->Charge()>0 && !fLikeSign) continue;
635 if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
638 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
641 if(postrack1->Charge()==negtrack1->Charge()) { // like-sign
642 isLikeSign2Prong=kTRUE;
643 if(!fLikeSign) continue;
644 if(iTrkN1<iTrkP1) continue; // this is needed to avoid double-counting of like-sign
645 } else { // unlike-sign
646 isLikeSign2Prong=kFALSE;
647 if(postrack1->Charge()<0 || negtrack1->Charge()>0) continue; // this is needed to avoid double-counting of unlike-sign
649 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
654 // back to primary vertex
655 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
656 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
658 // DCA between the two tracks
659 dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
660 if(dcap1n1>dcaMax) { negtrack1=0; continue; }
663 twoTrackArray1->AddAt(postrack1,0);
664 twoTrackArray1->AddAt(negtrack1,1);
665 AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
667 twoTrackArray1->Clear();
673 if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) {
675 io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
677 if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI) || (isLikeSign2Prong && (okD0 || okJPSI))) {
678 // add the vertex and the decay to the AOD
679 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
680 if(!isLikeSign2Prong) {
682 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
683 rd->SetSecondaryVtx(v2Prong);
684 v2Prong->SetParent(rd);
685 AddRefs(v2Prong,rd,event,twoTrackArray1);
688 rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
689 rd->SetSecondaryVtx(v2Prong);
690 if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
691 AddRefs(v2Prong,rd,event,twoTrackArray1);
693 } else { // isLikeSign2Prong
694 rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
695 rd->SetSecondaryVtx(v2Prong);
696 v2Prong->SetParent(rd);
697 AddRefs(v2Prong,rd,event,twoTrackArray1);
701 if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
702 // write references in io2Prong
704 AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
706 vertexp1n1->AddDaughter(postrack1);
707 vertexp1n1->AddDaughter(negtrack1);
709 io2Prong->SetSecondaryVtx(vertexp1n1);
710 //printf("---> %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
711 // create a track from the D0
712 AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
714 // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
715 for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
717 if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
719 if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
722 if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] ||
723 evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] ||
724 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
727 if(iTrkSoftPi%1==0) AliDebug(1,Form(" 1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));
729 trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
730 if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
732 // get track from tracks array
733 trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
734 trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
736 twoTrackArrayCasc->AddAt(trackPi,0);
737 twoTrackArrayCasc->AddAt(trackD0,1);
739 AliAODVertex *vertexCasc = 0;
741 if(fFindVertexForDstar) {
742 // DCA between the two tracks
743 dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
745 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
747 // assume Dstar decays at the primary vertex
748 Double_t pos[3],cov[6],chi2perNDF;
750 fV1->GetCovMatrix(cov);
751 chi2perNDF = fV1->GetChi2toNDF();
752 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
756 twoTrackArrayCasc->Clear();
761 ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
763 // add the D0 to the AOD (if not already done)
765 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
766 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
767 rd->SetSecondaryVtx(v2Prong);
768 v2Prong->SetParent(rd);
769 AddRefs(v2Prong,rd,event,twoTrackArray1);
770 okD0=kTRUE; // this is done to add it only once
772 // add the vertex and the cascade to the AOD
773 AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
774 rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
775 rc->SetSecondaryVtx(vCasc);
776 vCasc->SetParent(rc);
777 if(!fInputAOD) vCasc->AddDaughter(rd); // just to fill ref #0
778 AddRefs(vCasc,rc,event,twoTrackArrayCasc);
779 vCasc->AddDaughter(rd); // add the D0 (in ref #1)
781 twoTrackArrayCasc->Clear();
783 if(ioCascade) {delete ioCascade; ioCascade=NULL;}
784 delete vertexCasc; vertexCasc=NULL;
785 } // end loop on soft pi tracks
787 if(trackD0) {delete trackD0; trackD0=NULL;}
790 if(io2Prong) {delete io2Prong; io2Prong=NULL;}
793 twoTrackArray1->Clear();
794 if( (!f3Prong && !f4Prong) ||
795 (isLikeSign2Prong && !f3Prong) ) {
802 // 2nd LOOP ON POSITIVE TRACKS
803 for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
805 if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
807 if(iTrkP2%1==0) AliDebug(1,Form(" 2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));
809 // get track from tracks array
810 postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
812 if(postrack2->Charge()<0) continue;
814 if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
817 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
818 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
819 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
822 if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
823 if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
824 isLikeSign3Prong=kTRUE;
828 } else { // normal triplet (+-+)
829 isLikeSign3Prong=kFALSE;
831 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
832 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
833 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
837 // back to primary vertex
838 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
839 postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
840 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
841 //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
843 dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
844 if(dcap2n1>dcaMax) { postrack2=0; continue; }
845 dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
846 if(dcap1p2>dcaMax) { postrack2=0; continue; }
849 twoTrackArray2->AddAt(postrack2,0);
850 twoTrackArray2->AddAt(negtrack1,1);
851 AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
853 twoTrackArray2->Clear();
858 // 3 prong candidates
860 if(postrack2->Charge()>0) {
861 threeTrackArray->AddAt(postrack1,0);
862 threeTrackArray->AddAt(negtrack1,1);
863 threeTrackArray->AddAt(postrack2,2);
865 threeTrackArray->AddAt(negtrack1,0);
866 threeTrackArray->AddAt(postrack1,1);
867 threeTrackArray->AddAt(postrack2,2);
870 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
871 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
873 AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
874 if(!isLikeSign3Prong) {
875 rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
876 rd->SetSecondaryVtx(v3Prong);
877 v3Prong->SetParent(rd);
878 AddRefs(v3Prong,rd,event,threeTrackArray);
879 } else { // isLikeSign3Prong
880 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
881 rd->SetSecondaryVtx(v3Prong);
882 v3Prong->SetParent(rd);
883 AddRefs(v3Prong,rd,event,threeTrackArray);
886 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
887 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
890 // 4 prong candidates
892 // don't make 4 prong with like-sign pairs and triplets
893 && !isLikeSign2Prong && !isLikeSign3Prong
894 // track-to-track dca cuts already now
895 && dcap1n1 < fD0to4ProngsCuts[1]
896 && dcap2n1 < fD0to4ProngsCuts[1]) {
898 // back to primary vertex
899 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
900 postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
901 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
902 // Vertexing for these 3 (can be taken from above?)
903 threeTrackArray->AddAt(postrack1,0);
904 threeTrackArray->AddAt(negtrack1,1);
905 threeTrackArray->AddAt(postrack2,2);
906 AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
908 // 3rd LOOP ON NEGATIVE TRACKS (for 4 prong)
909 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
911 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
913 if(iTrkN2%1==0) AliDebug(1,Form(" 3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
915 // get track from tracks array
916 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
918 if(negtrack2->Charge()>0) continue;
920 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
922 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
923 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
924 evtNumber[iTrkP2]==evtNumber[iTrkN2] ||
925 evtNumber[iTrkP1]==evtNumber[iTrkN1] ||
926 evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
927 evtNumber[iTrkN1]==evtNumber[iTrkP2]) continue;
930 // back to primary vertex
931 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
932 postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
933 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
934 negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
935 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
936 if(dcap1n2 > fD0to4ProngsCuts[1]) { negtrack2=0; continue; }
937 dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
938 if(dcap2n2 > fD0to4ProngsCuts[1]) { negtrack2=0; continue; }
941 fourTrackArray->AddAt(postrack1,0);
942 fourTrackArray->AddAt(negtrack1,1);
943 fourTrackArray->AddAt(postrack2,2);
944 fourTrackArray->AddAt(negtrack2,3);
946 AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
947 io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
949 AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
950 rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
951 rd->SetSecondaryVtx(v4Prong);
952 v4Prong->SetParent(rd);
953 AddRefs(v4Prong,rd,event,fourTrackArray);
956 if(io4Prong) {delete io4Prong; io4Prong=NULL;}
957 if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;}
958 fourTrackArray->Clear();
961 } // end loop on negative tracks
963 threeTrackArray->Clear();
971 } // end 2nd loop on positive tracks
973 twoTrackArray2->Clear();
975 // 2nd LOOP ON NEGATIVE TRACKS (for 3 prong -+-)
976 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
978 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
980 if(iTrkN2%1==0) AliDebug(1,Form(" 2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
982 // get track from tracks array
983 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
985 if(negtrack2->Charge()>0) continue;
987 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
990 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
991 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
992 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
995 if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
996 if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
997 isLikeSign3Prong=kTRUE;
1001 } else { // normal triplet (-+-)
1002 isLikeSign3Prong=kFALSE;
1005 // back to primary vertex
1006 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1007 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1008 negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1009 //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
1011 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1012 if(dcap1n2>dcaMax) { negtrack2=0; continue; }
1013 dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1014 if(dcan1n2>dcaMax) { negtrack2=0; continue; }
1017 twoTrackArray2->AddAt(postrack1,0);
1018 twoTrackArray2->AddAt(negtrack2,1);
1020 AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1022 twoTrackArray2->Clear();
1028 threeTrackArray->AddAt(negtrack1,0);
1029 threeTrackArray->AddAt(postrack1,1);
1030 threeTrackArray->AddAt(negtrack2,2);
1031 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1032 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
1034 AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1035 if(!isLikeSign3Prong) {
1036 rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1037 rd->SetSecondaryVtx(v3Prong);
1038 v3Prong->SetParent(rd);
1039 AddRefs(v3Prong,rd,event,threeTrackArray);
1040 } else { // isLikeSign3Prong
1041 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1042 rd->SetSecondaryVtx(v3Prong);
1043 v3Prong->SetParent(rd);
1044 AddRefs(v3Prong,rd,event,threeTrackArray);
1047 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1048 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1054 } // end 2nd loop on negative tracks
1056 twoTrackArray2->Clear();
1060 } // end 1st loop on negative tracks
1063 } // end 1st loop on positive tracks
1066 AliDebug(1,Form(" Total HF vertices in event = %d;",
1067 (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
1069 AliDebug(1,Form(" D0->Kpi in event = %d;",
1070 (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
1073 AliDebug(1,Form(" JPSI->ee in event = %d;",
1074 (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
1077 AliDebug(1,Form(" Charm->3Prong in event = %d;",
1078 (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
1081 AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
1082 (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
1085 AliDebug(1,Form(" D*->D0pi in event = %d;\n",
1086 (Int_t)aodDstarTClArr->GetEntriesFast()));
1089 AliDebug(1,Form(" cascades -> v0 + track in event = %d;\n",
1090 (Int_t)aodCascadesTClArr->GetEntriesFast()));
1093 AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
1094 (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
1096 if(fLikeSign && f3Prong) {
1097 AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
1098 (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
1102 twoTrackArray1->Delete(); delete twoTrackArray1;
1103 twoTrackArray2->Delete(); delete twoTrackArray2;
1104 twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1105 twoTrackArrayV0->Delete(); delete twoTrackArrayV0;
1106 threeTrackArray->Clear();
1107 threeTrackArray->Delete(); delete threeTrackArray;
1108 fourTrackArray->Delete(); delete fourTrackArray;
1109 delete [] seleFlags; seleFlags=NULL;
1110 if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
1113 seleTrksArray.Delete();
1114 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
1119 //----------------------------------------------------------------------------
1120 void AliAnalysisVertexingHF::AddRefs(AliAODVertex *v,AliAODRecoDecayHF *rd,
1121 const AliVEvent *event,
1122 const TObjArray *trkArray) const
1124 // Add the AOD tracks as daughters of the vertex (TRef)
1125 // Also add the references to the primary vertex and to the cuts
1128 AddDaughterRefs(v,event,trkArray);
1129 rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1133 rd->SetListOfCutsRef((TList*)fListOfCuts);
1134 //fListOfCuts->Print();
1135 cout<<fListOfCuts<<endl;
1136 TList *l=(TList*)rd->GetListOfCuts();
1138 if(l) {l->Print(); }else{printf("error\n");}
1143 //----------------------------------------------------------------------------
1144 void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,
1145 const AliVEvent *event,
1146 const TObjArray *trkArray) const
1148 // Add the AOD tracks as daughters of the vertex (TRef)
1150 Int_t nDg = v->GetNDaughters();
1152 if(nDg) dg = v->GetDaughter(0);
1154 if(dg) return; // daughters already added
1156 Int_t nTrks = trkArray->GetEntriesFast();
1158 AliExternalTrackParam *track = 0;
1159 AliAODTrack *aodTrack = 0;
1162 for(Int_t i=0; i<nTrks; i++) {
1163 track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1164 id = (Int_t)track->GetID();
1165 //printf("---> %d\n",id);
1166 if(id<0) continue; // this track is a AliAODRecoDecay
1167 aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
1168 v->AddDaughter(aodTrack);
1173 //----------------------------------------------------------------------------
1174 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1175 TObjArray *twoTrackArray,AliVEvent *event,
1176 AliAODVertex *secVert,
1177 AliAODRecoDecayHF2Prong *rd2Prong,
1179 Bool_t &okDstar) const
1181 // Make the cascade as a 2Prong decay and check if it passes Dstar
1182 // reconstruction cuts
1186 Bool_t dummy1,dummy2,dummy3;
1188 // We use Make2Prong to construct the AliAODRecoCascadeHF
1189 // (which inherits from AliAODRecoDecayHF2Prong)
1190 AliAODRecoCascadeHF *theCascade =
1191 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1192 dummy1,dummy2,dummy3);
1193 if(!theCascade) return 0x0;
1196 AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1197 theCascade->SetCharge(trackPi->Charge());
1199 //--- selection cuts
1201 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1202 tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
1203 tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
1204 AliAODVertex *primVertexAOD=0;
1205 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1206 // take event primary vertex
1207 primVertexAOD = PrimaryVertex();
1208 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1209 rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
1213 //Bool_t testD0=kTRUE;
1214 //okDstar = tmpCascade->SelectDstar(fDstarCuts,fD0fromDstarCuts,testD0);
1215 okDstar = (Bool_t)fCutsDStartoKpipi->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1217 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1218 tmpCascade->UnsetOwnPrimaryVtx();
1219 delete tmpCascade; tmpCascade=NULL;
1220 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1221 rd2Prong->UnsetOwnPrimaryVtx();
1223 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1230 //----------------------------------------------------------------------------
1231 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1232 TObjArray *twoTrackArray,AliVEvent *event,
1233 AliAODVertex *secVert,
1236 Bool_t &okCascades) const
1239 // Make the cascade as a 2Prong decay and check if it passes
1240 // cascades reconstruction cuts
1242 // AliDebug(2,Form(" building the cascade"));
1244 Bool_t dummy1,dummy2,dummy3;
1246 // We use Make2Prong to construct the AliAODRecoCascadeHF
1247 // (which inherits from AliAODRecoDecayHF2Prong)
1248 AliAODRecoCascadeHF *theCascade =
1249 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1250 dummy1,dummy2,dummy3);
1251 if(!theCascade) return 0x0;
1253 // bachelor track and charge
1254 AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1255 theCascade->SetCharge(trackBachelor->Charge());
1257 //--- selection cuts
1259 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1260 tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
1261 tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
1262 AliAODVertex *primVertexAOD=0;
1263 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1264 // take event primary vertex
1265 primVertexAOD = PrimaryVertex();
1266 if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
1267 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1271 Bool_t okLcksp=0, okLcLpi=0;
1272 if(fCascades && fInputAOD){
1274 okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1275 if(okCascades==1) okLcksp=1;
1276 if(okCascades==2) okLcLpi=1;
1277 if(okCascades==3) { okLcksp=1; okLcLpi=1;}
1279 else okCascades = tmpCascade->SelectLctoV0(fLctoV0Cuts,okLcksp,okLcLpi);
1281 else { AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied")); okCascades=kTRUE; }// no cuts implemented from ESDs
1282 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1283 tmpCascade->UnsetOwnPrimaryVtx();
1284 delete tmpCascade; tmpCascade=NULL;
1285 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1291 //-----------------------------------------------------------------------------
1292 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
1293 TObjArray *twoTrackArray,AliVEvent *event,
1294 AliAODVertex *secVert,Double_t dca,
1295 Bool_t &okD0,Bool_t &okJPSI,
1296 Bool_t &okD0fromDstar) const
1298 // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
1299 // reconstruction cuts
1300 // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
1302 okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
1304 Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1305 AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1306 AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
1308 // propagate tracks to secondary vertex, to compute inv. mass
1309 postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1310 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1312 Double_t momentum[3];
1313 postrack->GetPxPyPz(momentum);
1314 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1315 negtrack->GetPxPyPz(momentum);
1316 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1319 // invariant mass cut (try to improve coding here..)
1320 Bool_t okMassCut=kFALSE;
1321 if(!okMassCut && fD0toKpi) if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
1322 if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
1323 if(!okMassCut && fDstar) if(SelectInvMass(3,2,px,py,pz)) okMassCut=kTRUE;
1325 AliDebug(2," candidate didn't pass mass cut");
1328 // primary vertex to be used by this candidate
1329 AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
1330 if(!primVertexAOD) return 0x0;
1333 Double_t d0z0[2],covd0z0[3];
1334 postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1336 d0err[0] = TMath::Sqrt(covd0z0[0]);
1337 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1339 d0err[1] = TMath::Sqrt(covd0z0[0]);
1341 // create the object AliAODRecoDecayHF2Prong
1342 AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
1343 the2Prong->SetOwnPrimaryVtx(primVertexAOD);
1344 UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1345 the2Prong->SetProngIDs(2,id);
1346 delete primVertexAOD; primVertexAOD=NULL;
1349 if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
1351 //Int_t checkD0,checkD0bar;
1352 //if(fD0toKpi) okD0 = the2Prong->SelectD0(fD0toKpiCuts,checkD0,checkD0bar);
1353 if(fD0toKpi) okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1354 //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
1355 // select J/psi from B
1357 //if(fJPSItoEle) okJPSI = the2Prong->SelectBtoJPSI(fBtoJPSICuts,checkJPSI);
1358 if(fJPSItoEle) okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1359 //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
1360 // select D0->Kpi from Dstar
1361 //if(fDstar) okD0fromDstar = the2Prong->SelectD0(fD0fromDstarCuts,checkD0,checkD0bar);
1362 if(fDstar) okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
1363 //if(fDebug && fDstar) printf(" %d\n",(Int_t)okD0fromDstar);
1367 // remove the primary vertex (was used only for selection)
1368 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1369 the2Prong->UnsetOwnPrimaryVtx();
1372 // get PID info from ESD
1373 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1374 if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1375 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1376 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1377 Double_t esdpid[10];
1378 for(Int_t i=0;i<5;i++) {
1379 esdpid[i] = esdpid0[i];
1380 esdpid[5+i] = esdpid1[i];
1382 the2Prong->SetPID(2,esdpid);
1386 //----------------------------------------------------------------------------
1387 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
1388 TObjArray *threeTrackArray,AliVEvent *event,
1389 AliAODVertex *secVert,Double_t dispersion,
1390 const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
1391 Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1392 Bool_t &ok3Prong) const
1394 // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1395 // reconstruction cuts
1400 if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
1402 Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1403 Double_t momentum[3];
1406 AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1407 AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1408 AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1410 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1411 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1412 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1413 postrack1->GetPxPyPz(momentum);
1414 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1415 negtrack->GetPxPyPz(momentum);
1416 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1417 postrack2->GetPxPyPz(momentum);
1418 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1420 // invariant mass cut for D+, Ds, Lc
1421 Bool_t okMassCut=kFALSE;
1422 if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
1424 AliDebug(2," candidate didn't pass mass cut");
1428 // primary vertex to be used by this candidate
1429 AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
1430 if(!primVertexAOD) return 0x0;
1432 Double_t d0z0[2],covd0z0[3];
1433 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1435 d0err[0] = TMath::Sqrt(covd0z0[0]);
1436 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1438 d0err[1] = TMath::Sqrt(covd0z0[0]);
1439 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1441 d0err[2] = TMath::Sqrt(covd0z0[0]);
1444 // create the object AliAODRecoDecayHF3Prong
1445 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1446 Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1447 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]));
1448 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]));
1449 Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
1450 AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1451 the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1452 UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1453 the3Prong->SetProngIDs(3,id);
1455 delete primVertexAOD; primVertexAOD=NULL;
1457 // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1462 //if(the3Prong->SelectDplus(fDplusCuts)) ok3Prong = kTRUE;
1463 //if(the3Prong->SelectDs(fDsCuts,ok1,ok2,dum1,dum2)) ok3Prong = kTRUE;
1464 //if(the3Prong->SelectLc(fLcCuts,ok1,ok2)) ok3Prong = kTRUE;
1466 if(fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1467 if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1468 if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1471 //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1473 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1474 the3Prong->UnsetOwnPrimaryVtx();
1477 // get PID info from ESD
1478 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1479 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1480 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1481 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1482 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1483 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1485 Double_t esdpid[15];
1486 for(Int_t i=0;i<5;i++) {
1487 esdpid[i] = esdpid0[i];
1488 esdpid[5+i] = esdpid1[i];
1489 esdpid[10+i] = esdpid2[i];
1491 the3Prong->SetPID(3,esdpid);
1495 //----------------------------------------------------------------------------
1496 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1497 TObjArray *fourTrackArray,AliVEvent *event,
1498 AliAODVertex *secVert,
1499 const AliAODVertex *vertexp1n1,
1500 const AliAODVertex *vertexp1n1p2,
1501 Double_t dcap1n1,Double_t dcap1n2,
1502 Double_t dcap2n1,Double_t dcap2n2,
1503 Bool_t &ok4Prong) const
1505 // Make 4Prong candidates and check if they pass D0toKpipipi
1506 // reconstruction cuts
1507 // G.E.Bruno, R.Romita
1510 if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
1512 Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1514 AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1515 AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1516 AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1517 AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1519 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1520 negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1521 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1522 negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1524 Double_t momentum[3];
1525 postrack1->GetPxPyPz(momentum);
1526 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1527 negtrack1->GetPxPyPz(momentum);
1528 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1529 postrack2->GetPxPyPz(momentum);
1530 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1531 negtrack2->GetPxPyPz(momentum);
1532 px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1534 // invariant mass cut for rho or D0 (try to improve coding here..)
1535 Bool_t okMassCut=kFALSE;
1536 if(!okMassCut && TMath::Abs(fD0to4ProngsCuts[8])<1.e-13){ //no PID, to be implemented with PID
1537 if(SelectInvMass(4,4,px,py,pz)) okMassCut=kTRUE;
1540 //if(fDebug) printf(" candidate didn't pass mass cut\n");
1541 //printf(" candidate didn't pass mass cut\n");
1545 // primary vertex to be used by this candidate
1546 AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
1547 if(!primVertexAOD) return 0x0;
1549 Double_t d0z0[2],covd0z0[3];
1550 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1552 d0err[0] = TMath::Sqrt(covd0z0[0]);
1553 negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1555 d0err[1] = TMath::Sqrt(covd0z0[0]);
1556 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1558 d0err[2] = TMath::Sqrt(covd0z0[0]);
1559 negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1561 d0err[3] = TMath::Sqrt(covd0z0[0]);
1564 // create the object AliAODRecoDecayHF4Prong
1565 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1566 Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
1567 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]));
1568 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]));
1569 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]));
1571 AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
1572 the4Prong->SetOwnPrimaryVtx(primVertexAOD);
1573 UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
1574 the4Prong->SetProngIDs(4,id);
1576 delete primVertexAOD; primVertexAOD=NULL;
1578 //Int_t checkD0,checkD0bar;
1579 //ok4Prong=the4Prong->SelectD0(fD0to4ProngsCuts,checkD0,checkD0bar);
1580 ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
1583 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1584 the4Prong->UnsetOwnPrimaryVtx();
1588 // get PID info from ESD
1589 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1590 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1591 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1592 if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
1593 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1594 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1595 Double_t esdpid3[5]={0.,0.,0.,0.,0.};
1596 if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
1598 Double_t esdpid[20];
1599 for(Int_t i=0;i<5;i++) {
1600 esdpid[i] = esdpid0[i];
1601 esdpid[5+i] = esdpid1[i];
1602 esdpid[10+i] = esdpid2[i];
1603 esdpid[15+i] = esdpid3[i];
1605 the4Prong->SetPID(4,esdpid);
1609 //-----------------------------------------------------------------------------
1610 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
1611 AliVEvent *event) const
1613 // Returns primary vertex to be used for this candidate
1615 AliESDVertex *vertexESD = 0;
1616 AliAODVertex *vertexAOD = 0;
1619 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1620 // primary vertex from the input event
1622 vertexESD = new AliESDVertex(*fV1);
1625 // primary vertex specific to this candidate
1627 Int_t nTrks = trkArray->GetEntriesFast();
1628 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1630 if(fRecoPrimVtxSkippingTrks) {
1631 // recalculating the vertex
1633 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1634 Float_t diamondcovxy[3];
1635 event->GetDiamondCovXY(diamondcovxy);
1636 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1637 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1638 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1639 vertexer->SetVtxStart(diamond);
1640 delete diamond; diamond=NULL;
1641 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1642 vertexer->SetOnlyFitter();
1644 Int_t skipped[1000];
1645 Int_t nTrksToSkip=0,id;
1646 AliExternalTrackParam *t = 0;
1647 for(Int_t i=0; i<nTrks; i++) {
1648 t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1649 id = (Int_t)t->GetID();
1651 skipped[nTrksToSkip++] = id;
1654 // For AOD, skip also tracks without covariance matrix
1656 Double_t covtest[21];
1657 for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1658 AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1659 if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1660 id = (Int_t)vtrack->GetID();
1662 skipped[nTrksToSkip++] = id;
1667 vertexer->SetSkipTracks(nTrksToSkip,skipped);
1668 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1670 } else if(fRmTrksFromPrimVtx) {
1671 // removing the prongs tracks
1673 TObjArray rmArray(nTrks);
1674 UShort_t *rmId = new UShort_t[nTrks];
1675 AliESDtrack *esdTrack = 0;
1677 for(Int_t i=0; i<nTrks; i++) {
1678 t = (AliESDtrack*)trkArray->UncheckedAt(i);
1679 esdTrack = new AliESDtrack(*t);
1680 rmArray.AddLast(esdTrack);
1681 if(esdTrack->GetID()>=0) {
1682 rmId[i]=(UShort_t)esdTrack->GetID();
1687 Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
1688 vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1689 delete [] rmId; rmId=NULL;
1694 if(!vertexESD) return vertexAOD;
1695 if(vertexESD->GetNContributors()<=0) {
1696 AliDebug(2,"vertexing failed");
1697 delete vertexESD; vertexESD=NULL;
1701 delete vertexer; vertexer=NULL;
1705 // convert to AliAODVertex
1706 Double_t pos[3],cov[6],chi2perNDF;
1707 vertexESD->GetXYZ(pos); // position
1708 vertexESD->GetCovMatrix(cov); //covariance matrix
1709 chi2perNDF = vertexESD->GetChi2toNDF();
1710 delete vertexESD; vertexESD=NULL;
1712 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1716 //-----------------------------------------------------------------------------
1717 void AliAnalysisVertexingHF::PrintStatus() const {
1718 // Print parameters being used
1720 //printf("Preselections:\n");
1721 // fTrackFilter->Dump();
1723 printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1725 printf("Secondary vertex with AliVertexerTracks\n");
1727 if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1728 if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1730 printf("Reconstruct D0->Kpi candidates with cuts:\n");
1731 if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
1734 printf("Reconstruct D*->D0pi candidates with cuts:\n");
1735 if(fFindVertexForDstar) {
1736 printf(" Reconstruct a secondary vertex for the D*\n");
1738 printf(" Assume the D* comes from the primary vertex\n");
1740 if(fCutsDStartoKpipi) fCutsDStartoKpipi->PrintAll();
1743 printf("Reconstruct J/psi from B candidates with cuts:\n");
1744 if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
1747 printf("Reconstruct 3 prong candidates.\n");
1748 printf(" D+->Kpipi cuts:\n");
1749 if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
1750 printf(" Ds->KKpi cuts:\n");
1751 if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
1752 printf(" Lc->pKpi cuts:\n");
1753 if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
1756 printf("Reconstruct 4 prong candidates.\n");
1757 printf(" D0->Kpipipi cuts:\n");
1758 if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
1761 printf("Reconstruct cascades candidates formed with v0s.\n");
1762 printf(" Lc -> k0s P & Lc -> L Pi cuts:\n");
1763 if(fCutsLctoV0) fCutsLctoV0->PrintAll();
1768 //-----------------------------------------------------------------------------
1769 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
1770 Double_t &dispersion,Bool_t useTRefArray) const
1772 // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1774 AliESDVertex *vertexESD = 0;
1775 AliAODVertex *vertexAOD = 0;
1777 if(!fSecVtxWithKF) { // AliVertexerTracks
1779 AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1780 vertexer->SetVtxStart(fV1);
1781 vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1782 delete vertexer; vertexer=NULL;
1784 if(!vertexESD) return vertexAOD;
1786 if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
1787 AliDebug(2,"vertexing failed");
1788 delete vertexESD; vertexESD=NULL;
1792 } else { // Kalman Filter vertexer (AliKFParticle)
1794 AliKFParticle::SetField(fBzkG);
1796 AliKFVertex vertexKF;
1798 Int_t nTrks = trkArray->GetEntriesFast();
1799 for(Int_t i=0; i<nTrks; i++) {
1800 AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1801 AliKFParticle daughterKF(*esdTrack,211);
1802 vertexKF.AddDaughter(daughterKF);
1804 vertexESD = new AliESDVertex(vertexKF.Parameters(),
1805 vertexKF.CovarianceMatrix(),
1807 vertexKF.GetNContributors());
1811 // convert to AliAODVertex
1812 Double_t pos[3],cov[6],chi2perNDF;
1813 vertexESD->GetXYZ(pos); // position
1814 vertexESD->GetCovMatrix(cov); //covariance matrix
1815 chi2perNDF = vertexESD->GetChi2toNDF();
1816 dispersion = vertexESD->GetDispersion();
1817 delete vertexESD; vertexESD=NULL;
1819 Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
1820 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
1824 //-----------------------------------------------------------------------------
1825 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1829 Double_t *pz) const {
1830 // Check invariant mass cut
1832 Short_t dummycharge=0;
1833 Double_t *dummyd0 = new Double_t[nprongs];
1834 for(Int_t ip=0;ip<nprongs;ip++) dummyd0[ip]=0.;
1835 AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
1836 delete [] dummyd0; dummyd0=NULL;
1838 UInt_t pdg2[2],pdg3[3],pdg4[4];
1841 Bool_t retval=kFALSE;
1845 pdg2[0]=211; pdg2[1]=321;
1846 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1847 minv = rd->InvMass(nprongs,pdg2);
1848 //if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1849 if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
1850 pdg2[0]=321; pdg2[1]=211;
1851 minv = rd->InvMass(nprongs,pdg2);
1852 //if(TMath::Abs(minv-mPDG)<fD0toKpiCuts[0]) retval=kTRUE;
1853 if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
1856 pdg2[0]=11; pdg2[1]=11;
1857 mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
1858 minv = rd->InvMass(nprongs,pdg2);
1859 //if(TMath::Abs(minv-mPDG)<fBtoJPSICuts[0]) retval=kTRUE;
1860 if(TMath::Abs(minv-mPDG)<fCutsJpsitoee->GetMassCut()) retval=kTRUE;
1862 case 2: // D+->Kpipi
1863 pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
1864 mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
1865 minv = rd->InvMass(nprongs,pdg3);
1866 //if(TMath::Abs(minv-mPDG)<fDplusCuts[0]) retval=kTRUE;
1867 if(TMath::Abs(minv-mPDG)<fCutsDplustoKpipi->GetMassCut()) retval=kTRUE;
1869 pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
1870 mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1871 minv = rd->InvMass(nprongs,pdg3);
1872 //if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1873 if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
1874 pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
1875 minv = rd->InvMass(nprongs,pdg3);
1876 //if(TMath::Abs(minv-mPDG)<fDsCuts[0]) retval=kTRUE;
1877 if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
1879 pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
1880 mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1881 minv = rd->InvMass(nprongs,pdg3);
1882 //if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1883 if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
1884 pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
1885 minv = rd->InvMass(nprongs,pdg3);
1886 //if(TMath::Abs(minv-mPDG)<fLcCuts[0]) retval=kTRUE;
1887 if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
1890 pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
1891 mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
1892 minv = rd->InvMass(nprongs,pdg2);
1893 //if(TMath::Abs(minv-mPDG)<fDstarCuts[0]) retval=kTRUE;
1894 if(TMath::Abs(minv-mPDG)<fCutsDStartoKpipi->GetMassCut()) retval=kTRUE;
1896 case 4: // D0->Kpipipi without PID
1897 pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
1898 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1899 minv = rd->InvMass(nprongs,pdg4);
1900 //if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1901 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1902 pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
1903 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1904 minv = rd->InvMass(nprongs,pdg4);
1905 //if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1906 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1907 pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
1908 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1909 minv = rd->InvMass(nprongs,pdg4);
1910 //if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1911 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1912 pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
1913 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1914 minv = rd->InvMass(nprongs,pdg4);
1915 //if(TMath::Abs(minv-mPDG)<fD0to4ProngsCuts[0]) retval=kTRUE;
1916 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1919 printf("SelectInvMass(): wrong decay selection\n");
1927 //-----------------------------------------------------------------------------
1928 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
1929 TObjArray &seleTrksArray,Int_t &nSeleTrks,
1930 UChar_t *seleFlags,Int_t *evtNumber)
1932 // Apply single-track preselection.
1933 // Fill a TObjArray with selected tracks (for displaced vertices or
1934 // soft pion from D*). Selection flag stored in seleFlags.
1935 // Create the AliESDVertex object (convert from AliAODVertex if necessary)
1936 // In case of AOD input, also fill fAODMap for track index<->ID
1938 const AliVVertex *vprimary = event->GetPrimaryVertex();
1940 if(fV1) { delete fV1; fV1=NULL; }
1941 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
1944 UShort_t *indices = 0;
1945 Double_t pos[3],cov[6];
1947 if(!fInputAOD) { // ESD
1948 fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
1950 vprimary->GetXYZ(pos);
1951 vprimary->GetCovarianceMatrix(cov);
1952 fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
1953 indices = new UShort_t[event->GetNumberOfTracks()];
1954 fAODMapSize = 100000;
1955 fAODMap = new Int_t[fAODMapSize];
1959 Int_t entries = (Int_t)event->GetNumberOfTracks();
1960 Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
1963 // transfer ITS tracks from event to arrays
1964 for(Int_t i=0; i<entries; i++) {
1966 track = (AliVTrack*)event->GetTrack(i);
1968 // skip pure ITS SA tracks
1969 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
1971 // TEMPORARY: check that the cov matrix is there
1972 Double_t covtest[21];
1973 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1977 AliAODTrack *aodt = (AliAODTrack*)track;
1978 if(aodt->GetUsedForPrimVtxFit()) {
1979 indices[nindices]=aodt->GetID(); nindices++;
1981 fAODMap[(Int_t)aodt->GetID()] = i;
1984 AliESDtrack *esdt = 0;
1987 esdt = (AliESDtrack*)track;
1989 esdt = new AliESDtrack(track);
1992 // single track selection
1993 okDisplaced=kFALSE; okSoftPi=kFALSE;
1995 evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
1996 const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
1997 Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
1998 eventVtx->GetXYZ(vtxPos);
1999 vprimary->GetXYZ(primPos);
2000 eventVtx->GetCovarianceMatrix(primCov);
2001 for(Int_t ind=0;ind<3;ind++){
2002 trasl[ind]=vtxPos[ind]-primPos[ind];
2005 Bool_t isTransl=esdt->Translate(trasl,primCov);
2013 if(SingleTrkCuts(esdt,okDisplaced,okSoftPi)) {
2014 seleTrksArray.AddLast(esdt);
2015 seleFlags[nSeleTrks]=0;
2016 if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2017 if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2020 if(fInputAOD) delete esdt;
2025 } // end loop on tracks
2027 // primary vertex from AOD
2029 delete fV1; fV1=NULL;
2030 vprimary->GetXYZ(pos);
2031 vprimary->GetCovarianceMatrix(cov);
2032 Double_t chi2toNDF = vprimary->GetChi2perNDF();
2033 Int_t ncontr=nindices;
2034 if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
2035 Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
2036 fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2037 fV1->SetTitle(vprimary->GetTitle());
2038 fV1->SetIndices(nindices,indices);
2040 if(indices) { delete [] indices; indices=NULL; }
2045 //-----------------------------------------------------------------------------
2046 void AliAnalysisVertexingHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
2047 Double_t cut2,Double_t cut3,Double_t cut4,
2048 Double_t cut5,Double_t cut6,
2049 Double_t cut7,Double_t cut8)
2051 // Set the cuts for D0 selection
2052 fD0toKpiCuts[0] = cut0;
2053 fD0toKpiCuts[1] = cut1;
2054 fD0toKpiCuts[2] = cut2;
2055 fD0toKpiCuts[3] = cut3;
2056 fD0toKpiCuts[4] = cut4;
2057 fD0toKpiCuts[5] = cut5;
2058 fD0toKpiCuts[6] = cut6;
2059 fD0toKpiCuts[7] = cut7;
2060 fD0toKpiCuts[8] = cut8;
2064 //-----------------------------------------------------------------------------
2065 void AliAnalysisVertexingHF::SetD0toKpiCuts(const Double_t cuts[9])
2067 // Set the cuts for D0 selection
2069 for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
2073 //-----------------------------------------------------------------------------
2074 void AliAnalysisVertexingHF::SetD0fromDstarCuts(Double_t cut0,Double_t cut1,
2075 Double_t cut2,Double_t cut3,Double_t cut4,
2076 Double_t cut5,Double_t cut6,
2077 Double_t cut7,Double_t cut8)
2079 // Set the cuts for D0 from D* selection
2080 fD0fromDstarCuts[0] = cut0;
2081 fD0fromDstarCuts[1] = cut1;
2082 fD0fromDstarCuts[2] = cut2;
2083 fD0fromDstarCuts[3] = cut3;
2084 fD0fromDstarCuts[4] = cut4;
2085 fD0fromDstarCuts[5] = cut5;
2086 fD0fromDstarCuts[6] = cut6;
2087 fD0fromDstarCuts[7] = cut7;
2088 fD0fromDstarCuts[8] = cut8;
2092 //-----------------------------------------------------------------------------
2093 void AliAnalysisVertexingHF::SetD0fromDstarCuts(const Double_t cuts[9])
2095 // Set the cuts for D0 from D* selection
2097 for(Int_t i=0; i<9; i++) fD0fromDstarCuts[i] = cuts[i];
2101 //-----------------------------------------------------------------------------
2102 void AliAnalysisVertexingHF::SetDstarCuts(Double_t cut0,Double_t cut1,
2103 Double_t cut2,Double_t cut3,
2106 // Set the cuts for D* selection
2107 fDstarCuts[0] = cut0;
2108 fDstarCuts[1] = cut1;
2109 fDstarCuts[2] = cut2;
2110 fDstarCuts[3] = cut3;
2111 fDstarCuts[4] = cut4;
2115 //-----------------------------------------------------------------------------
2116 void AliAnalysisVertexingHF::SetDstarCuts(const Double_t cuts[5])
2118 // Set the cuts for D* selection
2120 for(Int_t i=0; i<5; i++) fDstarCuts[i] = cuts[i];
2124 //-----------------------------------------------------------------------------
2125 void AliAnalysisVertexingHF::SetBtoJPSICuts(Double_t cut0,Double_t cut1,
2126 Double_t cut2,Double_t cut3,Double_t cut4,
2127 Double_t cut5,Double_t cut6,
2128 Double_t cut7,Double_t cut8)
2130 // Set the cuts for J/psi from B selection
2131 fBtoJPSICuts[0] = cut0;
2132 fBtoJPSICuts[1] = cut1;
2133 fBtoJPSICuts[2] = cut2;
2134 fBtoJPSICuts[3] = cut3;
2135 fBtoJPSICuts[4] = cut4;
2136 fBtoJPSICuts[5] = cut5;
2137 fBtoJPSICuts[6] = cut6;
2138 fBtoJPSICuts[7] = cut7;
2139 fBtoJPSICuts[8] = cut8;
2143 //-----------------------------------------------------------------------------
2144 void AliAnalysisVertexingHF::SetBtoJPSICuts(const Double_t cuts[9])
2146 // Set the cuts for J/psi from B selection
2148 for(Int_t i=0; i<9; i++) fBtoJPSICuts[i] = cuts[i];
2152 //-----------------------------------------------------------------------------
2153 void AliAnalysisVertexingHF::SetDplusCuts(Double_t cut0,Double_t cut1,
2154 Double_t cut2,Double_t cut3,Double_t cut4,
2155 Double_t cut5,Double_t cut6,
2156 Double_t cut7,Double_t cut8,
2157 Double_t cut9,Double_t cut10,Double_t cut11)
2159 // Set the cuts for Dplus->Kpipi selection
2160 fDplusCuts[0] = cut0;
2161 fDplusCuts[1] = cut1;
2162 fDplusCuts[2] = cut2;
2163 fDplusCuts[3] = cut3;
2164 fDplusCuts[4] = cut4;
2165 fDplusCuts[5] = cut5;
2166 fDplusCuts[6] = cut6;
2167 fDplusCuts[7] = cut7;
2168 fDplusCuts[8] = cut8;
2169 fDplusCuts[9] = cut9;
2170 fDplusCuts[10] = cut10;
2171 fDplusCuts[11] = cut11;
2175 //-----------------------------------------------------------------------------
2176 void AliAnalysisVertexingHF::SetDplusCuts(const Double_t cuts[12])
2178 // Set the cuts for Dplus->Kpipi selection
2180 for(Int_t i=0; i<12; i++) fDplusCuts[i] = cuts[i];
2184 //-----------------------------------------------------------------------------
2185 void AliAnalysisVertexingHF::SetDsCuts(Double_t cut0,Double_t cut1,
2186 Double_t cut2,Double_t cut3,
2187 Double_t cut4,Double_t cut5,
2188 Double_t cut6,Double_t cut7,
2189 Double_t cut8,Double_t cut9,
2190 Double_t cut10,Double_t cut11,
2191 Double_t cut12,Double_t cut13)
2193 // Set the cuts for Ds->KKpi selection
2204 fDsCuts[10] = cut10;
2205 fDsCuts[11] = cut11;
2206 fDsCuts[12] = cut12;
2207 fDsCuts[13] = cut13;
2211 //-----------------------------------------------------------------------------
2212 void AliAnalysisVertexingHF::SetDsCuts(const Double_t cuts[13])
2214 // Set the cuts for Ds->KKpi selection
2216 for(Int_t i=0; i<14; i++) fDsCuts[i] = cuts[i];
2220 //-----------------------------------------------------------------------------
2221 void AliAnalysisVertexingHF::SetLcCuts(Double_t cut0,Double_t cut1,
2222 Double_t cut2,Double_t cut3,Double_t cut4,
2223 Double_t cut5,Double_t cut6,
2224 Double_t cut7,Double_t cut8,
2225 Double_t cut9,Double_t cut10,Double_t cut11)
2227 // Set the cuts for Lc->pKpi selection
2238 fLcCuts[10] = cut10;
2239 fLcCuts[11] = cut11;
2243 //-----------------------------------------------------------------------------
2244 void AliAnalysisVertexingHF::SetLcCuts(const Double_t cuts[12])
2246 // Set the cuts for Lc->pKpi selection
2248 for(Int_t i=0; i<12; i++) fLcCuts[i] = cuts[i];
2253 //-----------------------------------------------------------------------------
2254 void AliAnalysisVertexingHF::SetLctoV0Cuts(Double_t cut0,Double_t cut1,
2255 Double_t cut2,Double_t cut3,Double_t cut4,
2256 Double_t cut5,Double_t cut6,
2257 Double_t cut7,Double_t cut8)
2259 // Set the cuts for Lc->V0+bachelor selection
2260 fLctoV0Cuts[0] = cut0;
2261 fLctoV0Cuts[1] = cut1;
2262 fLctoV0Cuts[2] = cut2;
2263 fLctoV0Cuts[3] = cut3;
2264 fLctoV0Cuts[4] = cut4;
2265 fLctoV0Cuts[5] = cut5;
2266 fLctoV0Cuts[6] = cut6;
2267 fLctoV0Cuts[7] = cut7;
2268 fLctoV0Cuts[8] = cut8;
2273 //-----------------------------------------------------------------------------
2274 void AliAnalysisVertexingHF::SetLctoV0Cuts(const Double_t cuts[8])
2276 // Set the cuts for Lc-> V0 + bachelor selection
2278 for(Int_t i=0; i<8; i++) fLctoV0Cuts[i] = cuts[i];
2283 //-----------------------------------------------------------------------------
2284 void AliAnalysisVertexingHF::SetD0to4ProngsCuts(Double_t cut0,Double_t cut1,
2285 Double_t cut2,Double_t cut3,Double_t cut4,
2286 Double_t cut5,Double_t cut6,
2287 Double_t cut7,Double_t cut8)
2289 // Set the cuts for D0->Kpipipi selection
2291 fD0to4ProngsCuts[0] = cut0;
2292 fD0to4ProngsCuts[1] = cut1;
2293 fD0to4ProngsCuts[2] = cut2;
2294 fD0to4ProngsCuts[3] = cut3;
2295 fD0to4ProngsCuts[4] = cut4;
2296 fD0to4ProngsCuts[5] = cut5;
2297 fD0to4ProngsCuts[6] = cut6;
2298 fD0to4ProngsCuts[7] = cut7;
2299 fD0to4ProngsCuts[8] = cut8;
2303 //-----------------------------------------------------------------------------
2304 void AliAnalysisVertexingHF::SetD0to4ProngsCuts(const Double_t cuts[9])
2306 // Set the cuts for D0->Kpipipi selection
2308 for(Int_t i=0; i<9; i++) fD0to4ProngsCuts[i] = cuts[i];
2312 //-----------------------------------------------------------------------------
2313 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2314 Bool_t &okDisplaced,Bool_t &okSoftPi) const
2316 // Check if track passes some kinematical cuts
2318 // this is needed to store the impact parameters
2319 trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2323 // Track selection, displaced tracks
2326 selectInfo = fTrackFilter->IsSelected(trk);
2328 if(selectInfo) okDisplaced=kTRUE;
2329 // Track selection, soft pions
2331 if(fDstar && fTrackFilterSoftPi) {
2332 selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2334 if(selectInfo) okSoftPi=kTRUE;
2336 if(okDisplaced || okSoftPi) return kTRUE;
2342 //-----------------------------------------------------------------------------
2343 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2345 // Transform ESDv0 to AODv0
2347 // this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2348 // and creates an AODv0 out of them
2350 Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
2351 AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2353 // create the v0 neutral track to compute the DCA to the primary vertex
2354 Double_t xyz[3], pxpypz[3];
2356 esdV0->PxPyPz(pxpypz);
2357 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2358 AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2359 if(!trackesdV0) return 0;
2360 Double_t d0z0[2],covd0z0[3];
2361 AliAODVertex *primVertexAOD = PrimaryVertex();
2362 trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2363 Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2364 // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2365 Double_t dcaV0DaughterToPrimVertex[2];
2366 AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2367 AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
2368 if( !posV0track || !negV0track) {
2369 if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2372 posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2373 // if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2375 dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
2376 negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2377 // if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2379 dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
2380 Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
2381 Double_t pmom[3],nmom[3];
2382 esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2383 esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2385 AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2386 aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2388 if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2389 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2393 //-----------------------------------------------------------------------------