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 **************************************************************************/
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.
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>
31 #include <TDatabasePDG.h>
35 #include "AliVEvent.h"
36 #include "AliVVertex.h"
37 #include "AliVTrack.h"
38 #include "AliVertexerTracks.h"
39 #include "AliKFVertex.h"
40 #include "AliESDEvent.h"
41 #include "AliESDVertex.h"
42 #include "AliExternalTrackParam.h"
43 #include "AliNeutralTrackParam.h"
44 #include "AliESDtrack.h"
45 #include "AliESDtrackCuts.h"
46 #include "AliAODEvent.h"
47 #include "AliAODRecoDecay.h"
48 #include "AliAODRecoDecayHF.h"
49 #include "AliAODRecoDecayHF2Prong.h"
50 #include "AliAODRecoDecayHF3Prong.h"
51 #include "AliAODRecoDecayHF4Prong.h"
52 #include "AliAODRecoCascadeHF.h"
53 #include "AliRDHFCutsD0toKpi.h"
54 #include "AliRDHFCutsJpsitoee.h"
55 #include "AliRDHFCutsDplustoKpipi.h"
56 #include "AliRDHFCutsDstoKKpi.h"
57 #include "AliRDHFCutsLctopKpi.h"
58 #include "AliRDHFCutsLctoV0.h"
59 #include "AliRDHFCutsD0toKpipipi.h"
60 #include "AliRDHFCutsDStartoKpipi.h"
61 #include "AliAnalysisFilter.h"
62 #include "AliAnalysisVertexingHF.h"
63 #include "AliMixedEvent.h"
67 ClassImp(AliAnalysisVertexingHF)
69 //----------------------------------------------------------------------------
70 AliAnalysisVertexingHF::AliAnalysisVertexingHF():
75 fSecVtxWithKF(kFALSE),
76 fRecoPrimVtxSkippingTrks(kFALSE),
77 fRmTrksFromPrimVtx(kFALSE),
86 fLikeSign3prong(kFALSE),
89 fTrackFilterSoftPi(0x0),
92 fCutsDplustoKpipi(0x0),
96 fCutsD0toKpipipi(0x0),
97 fCutsDStartoKpipi(0x0),
99 fFindVertexForDstar(kTRUE),
100 fFindVertexForCascades(kTRUE),
101 fMassCutBeforeVertexing(kFALSE),
108 // Default constructor
110 Double_t d02[2]={0.,0.};
111 Double_t d03[3]={0.,0.,0.};
112 Double_t d04[4]={0.,0.,0.,0.};
113 fMassCalc2 = new AliAODRecoDecay(0x0,2,0,d02);
114 fMassCalc3 = new AliAODRecoDecay(0x0,3,1,d03);
115 fMassCalc4 = new AliAODRecoDecay(0x0,4,0,d04);
118 //--------------------------------------------------------------------------
119 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
121 fInputAOD(source.fInputAOD),
122 fAODMapSize(source.fAODMapSize),
123 fAODMap(source.fAODMap),
125 fSecVtxWithKF(source.fSecVtxWithKF),
126 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
127 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
129 fD0toKpi(source.fD0toKpi),
130 fJPSItoEle(source.fJPSItoEle),
131 f3Prong(source.f3Prong),
132 f4Prong(source.f4Prong),
133 fDstar(source.fDstar),
134 fCascades(source.fCascades),
135 fLikeSign(source.fLikeSign),
136 fLikeSign3prong(source.fLikeSign3prong),
137 fMixEvent(source.fMixEvent),
138 fTrackFilter(source.fTrackFilter),
139 fTrackFilterSoftPi(source.fTrackFilterSoftPi),
140 fCutsD0toKpi(source.fCutsD0toKpi),
141 fCutsJpsitoee(source.fCutsJpsitoee),
142 fCutsDplustoKpipi(source.fCutsDplustoKpipi),
143 fCutsDstoKKpi(source.fCutsDstoKKpi),
144 fCutsLctopKpi(source.fCutsLctopKpi),
145 fCutsLctoV0(source.fCutsLctoV0),
146 fCutsD0toKpipipi(source.fCutsD0toKpipipi),
147 fCutsDStartoKpipi(source.fCutsDStartoKpipi),
148 fListOfCuts(source.fListOfCuts),
149 fFindVertexForDstar(source.fFindVertexForDstar),
150 fFindVertexForCascades(source.fFindVertexForCascades),
151 fMassCutBeforeVertexing(source.fMassCutBeforeVertexing),
152 fMassCalc2(source.fMassCalc2),
153 fMassCalc3(source.fMassCalc3),
154 fMassCalc4(source.fMassCalc4),
162 //--------------------------------------------------------------------------
163 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
166 // assignment operator
168 if(&source == this) return *this;
169 fInputAOD = source.fInputAOD;
170 fAODMapSize = source.fAODMapSize;
171 fBzkG = source.fBzkG;
172 fSecVtxWithKF = source.fSecVtxWithKF;
173 fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
174 fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
176 fD0toKpi = source.fD0toKpi;
177 fJPSItoEle = source.fJPSItoEle;
178 f3Prong = source.f3Prong;
179 f4Prong = source.f4Prong;
180 fDstar = source.fDstar;
181 fCascades = source.fCascades;
182 fLikeSign = source.fLikeSign;
183 fLikeSign3prong = source.fLikeSign3prong;
184 fMixEvent = source.fMixEvent;
185 fTrackFilter = source.fTrackFilter;
186 fTrackFilterSoftPi = source.fTrackFilterSoftPi;
187 fCutsD0toKpi = source.fCutsD0toKpi;
188 fCutsJpsitoee = source.fCutsJpsitoee;
189 fCutsDplustoKpipi = source.fCutsDplustoKpipi;
190 fCutsDstoKKpi = source.fCutsDstoKKpi;
191 fCutsLctopKpi = source.fCutsLctopKpi;
192 fCutsLctoV0 = source.fCutsLctoV0;
193 fCutsD0toKpipipi = source.fCutsD0toKpipipi;
194 fCutsDStartoKpipi = source.fCutsDStartoKpipi;
195 fListOfCuts = source.fListOfCuts;
196 fFindVertexForDstar = source.fFindVertexForDstar;
197 fFindVertexForCascades = source.fFindVertexForCascades;
198 fMassCutBeforeVertexing = source.fMassCutBeforeVertexing;
199 fMassCalc2 = source.fMassCalc2;
200 fMassCalc3 = source.fMassCalc3;
201 fMassCalc4 = source.fMassCalc4;
205 //----------------------------------------------------------------------------
206 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
208 if(fV1) { delete fV1; fV1=0; }
209 if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
210 if(fTrackFilterSoftPi) { delete fTrackFilterSoftPi; fTrackFilterSoftPi=0; }
211 if(fCutsD0toKpi) { delete fCutsD0toKpi; fCutsD0toKpi=0; }
212 if(fCutsJpsitoee) { delete fCutsJpsitoee; fCutsJpsitoee=0; }
213 if(fCutsDplustoKpipi) { delete fCutsDplustoKpipi; fCutsDplustoKpipi=0; }
214 if(fCutsDstoKKpi) { delete fCutsDstoKKpi; fCutsDstoKKpi=0; }
215 if(fCutsLctopKpi) { delete fCutsLctopKpi; fCutsLctopKpi=0; }
216 if(fCutsLctoV0) { delete fCutsLctoV0; fCutsLctoV0=0; }
217 if(fCutsD0toKpipipi) { delete fCutsD0toKpipipi; fCutsD0toKpipipi=0; }
218 if(fCutsDStartoKpipi) { delete fCutsDStartoKpipi; fCutsDStartoKpipi=0; }
219 if(fAODMap) { delete [] fAODMap; fAODMap=0; }
220 if(fMassCalc2) { delete fMassCalc2; fMassCalc2=0; }
221 if(fMassCalc3) { delete fMassCalc3; fMassCalc3=0; }
222 if(fMassCalc4) { delete fMassCalc4; fMassCalc4=0; }
224 //----------------------------------------------------------------------------
225 TList *AliAnalysisVertexingHF::FillListOfCuts() {
226 // Fill list of analysis cuts
228 TList *list = new TList();
230 list->SetName("ListOfCuts");
233 AliRDHFCutsD0toKpi *cutsD0toKpi = new AliRDHFCutsD0toKpi(*fCutsD0toKpi);
234 list->Add(cutsD0toKpi);
237 AliRDHFCutsJpsitoee *cutsJpsitoee = new AliRDHFCutsJpsitoee(*fCutsJpsitoee);
238 list->Add(cutsJpsitoee);
240 if(fCutsDplustoKpipi) {
241 AliRDHFCutsDplustoKpipi *cutsDplustoKpipi = new AliRDHFCutsDplustoKpipi(*fCutsDplustoKpipi);
242 list->Add(cutsDplustoKpipi);
245 AliRDHFCutsDstoKKpi *cutsDstoKKpi = new AliRDHFCutsDstoKKpi(*fCutsDstoKKpi);
246 list->Add(cutsDstoKKpi);
249 AliRDHFCutsLctopKpi *cutsLctopKpi = new AliRDHFCutsLctopKpi(*fCutsLctopKpi);
250 list->Add(cutsLctopKpi);
253 AliRDHFCutsLctoV0 *cutsLctoV0 = new AliRDHFCutsLctoV0(*fCutsLctoV0);
254 list->Add(cutsLctoV0);
256 if(fCutsD0toKpipipi) {
257 AliRDHFCutsD0toKpipipi *cutsD0toKpipipi = new AliRDHFCutsD0toKpipipi(*fCutsD0toKpipipi);
258 list->Add(cutsD0toKpipipi);
260 if(fCutsDStartoKpipi) {
261 AliRDHFCutsDStartoKpipi *cutsDStartoKpipi = new AliRDHFCutsDStartoKpipi(*fCutsDStartoKpipi);
262 list->Add(cutsDStartoKpipi);
265 // keep a pointer to the list
270 //----------------------------------------------------------------------------
271 void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
272 TClonesArray *aodVerticesHFTClArr,
273 TClonesArray *aodD0toKpiTClArr,
274 TClonesArray *aodJPSItoEleTClArr,
275 TClonesArray *aodCharm3ProngTClArr,
276 TClonesArray *aodCharm4ProngTClArr,
277 TClonesArray *aodDstarTClArr,
278 TClonesArray *aodCascadesTClArr,
279 TClonesArray *aodLikeSign2ProngTClArr,
280 TClonesArray *aodLikeSign3ProngTClArr)
282 // Find heavy-flavour vertex candidates
284 // Output: AOD (additional branches added)
287 TString evtype = event->IsA()->GetName();
288 fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
289 } // if we do mixing AliVEvent is a AliMixedEvent
292 AliDebug(2,"Creating HF candidates from AOD");
294 AliDebug(2,"Creating HF candidates from ESD");
297 if(!aodVerticesHFTClArr) {
298 printf("ERROR: no aodVerticesHFTClArr");
301 if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
302 printf("ERROR: no aodD0toKpiTClArr");
305 if(fJPSItoEle && !aodJPSItoEleTClArr) {
306 printf("ERROR: no aodJPSItoEleTClArr");
309 if(f3Prong && !aodCharm3ProngTClArr) {
310 printf("ERROR: no aodCharm3ProngTClArr");
313 if(f4Prong && !aodCharm4ProngTClArr) {
314 printf("ERROR: no aodCharm4ProngTClArr");
317 if(fDstar && !aodDstarTClArr) {
318 printf("ERROR: no aodDstarTClArr");
321 if(fCascades && !aodCascadesTClArr){
322 printf("ERROR: no aodCascadesTClArr ");
325 if(fLikeSign && !aodLikeSign2ProngTClArr) {
326 printf("ERROR: no aodLikeSign2ProngTClArr");
329 if(fLikeSign3prong && f3Prong && !aodLikeSign3ProngTClArr) {
330 printf("ERROR: no aodLikeSign3ProngTClArr");
334 // delete candidates from previous event and create references
335 Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0,iCascades=0,iLikeSign2Prong=0,iLikeSign3Prong=0;
336 aodVerticesHFTClArr->Delete();
337 iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
338 TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
339 if(fD0toKpi || fDstar) {
340 aodD0toKpiTClArr->Delete();
341 iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
344 aodJPSItoEleTClArr->Delete();
345 iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
348 aodCharm3ProngTClArr->Delete();
349 i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
352 aodCharm4ProngTClArr->Delete();
353 i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
356 aodDstarTClArr->Delete();
357 iDstar = aodDstarTClArr->GetEntriesFast();
360 aodCascadesTClArr->Delete();
361 iCascades = aodCascadesTClArr->GetEntriesFast();
364 aodLikeSign2ProngTClArr->Delete();
365 iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast();
367 if(fLikeSign3prong && f3Prong) {
368 aodLikeSign3ProngTClArr->Delete();
369 iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast();
372 TClonesArray &aodD0toKpiRef = *aodD0toKpiTClArr;
373 TClonesArray &aodJPSItoEleRef = *aodJPSItoEleTClArr;
374 TClonesArray &aodCharm3ProngRef = *aodCharm3ProngTClArr;
375 TClonesArray &aodCharm4ProngRef = *aodCharm4ProngTClArr;
376 TClonesArray &aodDstarRef = *aodDstarTClArr;
377 TClonesArray &aodCascadesRef = *aodCascadesTClArr;
378 TClonesArray &aodLikeSign2ProngRef = *aodLikeSign2ProngTClArr;
379 TClonesArray &aodLikeSign3ProngRef = *aodLikeSign3ProngTClArr;
382 AliAODRecoDecayHF2Prong *io2Prong = 0;
383 AliAODRecoDecayHF3Prong *io3Prong = 0;
384 AliAODRecoDecayHF4Prong *io4Prong = 0;
385 AliAODRecoCascadeHF *ioCascade = 0;
387 Int_t iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries,iv0,nv0;
388 Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaV0,dcaCasc;
389 Bool_t okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
390 Bool_t okDstar=kFALSE,okD0fromDstar=kFALSE;
391 Bool_t okCascades=kFALSE;
392 AliESDtrack *postrack1 = 0;
393 AliESDtrack *postrack2 = 0;
394 AliESDtrack *negtrack1 = 0;
395 AliESDtrack *negtrack2 = 0;
396 AliESDtrack *trackPi = 0;
397 // AliESDtrack *posV0track = 0;
398 // AliESDtrack *negV0track = 0;
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));
415 fnTrksTotal += trkEntries;
417 nv0 = (Int_t)event->GetNumberOfV0s();
418 AliDebug(1,Form(" Number of V0s: %d",nv0));
420 if( trkEntries<2 && (trkEntries<1 || nv0<1) ) {
421 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
426 if(!fCutsD0toKpi->IsEventSelected(event)) return;
428 // call function that applies sigle-track selection,
429 // for displaced tracks and soft pions (both charges) for D*,
430 // and retrieves primary vertex
431 TObjArray seleTrksArray(trkEntries);
432 TObjArray tracksAtVertex(trkEntries);
433 UChar_t *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi
435 Int_t *evtNumber = new Int_t[trkEntries];
436 SelectTracksAndCopyVertex(event,trkEntries,seleTrksArray,tracksAtVertex,nSeleTrks,seleFlags,evtNumber);
438 AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
439 fnSeleTrksTotal += nSeleTrks;
442 TObjArray *twoTrackArray1 = new TObjArray(2);
443 TObjArray *twoTrackArray2 = new TObjArray(2);
444 TObjArray *twoTrackArrayV0 = new TObjArray(2);
445 TObjArray *twoTrackArrayCasc = new TObjArray(2);
446 TObjArray *threeTrackArray = new TObjArray(3);
447 TObjArray *fourTrackArray = new TObjArray(4);
450 Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
452 AliAODRecoDecayHF *rd = 0;
453 AliAODRecoCascadeHF *rc = 0;
457 Bool_t massCutOK=kTRUE;
459 // LOOP ON POSITIVE TRACKS
460 for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
462 //if(iTrkP1%1==0) AliDebug(1,Form(" 1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));
463 //if(iTrkP1%1==0) printf(" 1st loop on pos: track number %d of %d\n",iTrkP1,nSeleTrks);
465 // get track from tracks array
466 postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
468 if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
470 // Make cascades with V0+track
474 for(iv0=0; iv0<nv0; iv0++){
476 //AliDebug(1,Form(" loop on v0s for track number %d and v0 number %d",iTrkP1,iv0));
480 v0 = ((AliAODEvent*)event)->GetV0(iv0);
482 esdV0 = ((AliESDEvent*)event)->GetV0(iv0);
484 if ( (!v0 || !v0->IsA()->InheritsFrom("AliAODv0") ) &&
485 (!esdV0 || !esdV0->IsA()->InheritsFrom("AliESDv0") ) ) continue;
488 // Get the tracks that form the V0
489 // ( parameters at primary vertex )
490 // and define an AliExternalTrackParam out of them
491 AliExternalTrackParam * posV0track;
492 AliExternalTrackParam * negV0track;
495 AliAODTrack *posVV0track = (AliAODTrack*)(v0->GetDaughter(0));
496 AliAODTrack *negVV0track = (AliAODTrack*)(v0->GetDaughter(1));
497 if( !posVV0track || !negVV0track ) continue;
499 // Apply some basic V0 daughter criteria
501 // bachelor must not be a v0-track
502 if (posVV0track->GetID() == postrack1->GetID() ||
503 negVV0track->GetID() == postrack1->GetID()) continue;
504 // reject like-sign v0
505 if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
506 // avoid ghost TPC tracks
507 if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
508 !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
509 // Get AliExternalTrackParam out of the AliAODTracks
510 Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign;
511 posVV0track->PxPyPz(pxpypz); posVV0track->XvYvZv(xyz);
512 posVV0track->GetCovarianceXYZPxPyPz(cv); sign=posVV0track->Charge();
513 posV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
514 negVV0track->PxPyPz(pxpypz); negVV0track->XvYvZv(xyz);
515 negVV0track->GetCovarianceXYZPxPyPz(cv); sign=negVV0track->Charge();
516 negV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
518 AliESDtrack *posVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetPindex() ));
519 AliESDtrack *negVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetNindex() ));
520 if( !posVV0track || !negVV0track ) continue;
522 // Apply some basic V0 daughter criteria
524 // bachelor must not be a v0-track
525 if (posVV0track->GetID() == postrack1->GetID() ||
526 negVV0track->GetID() == postrack1->GetID()) continue;
527 // reject like-sign v0
528 if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
529 // avoid ghost TPC tracks
530 if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
531 !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
532 // reject kinks (only necessary on AliESDtracks)
533 if (posVV0track->GetKinkIndex(0)>0 || negVV0track->GetKinkIndex(0)>0) continue;
534 // Get AliExternalTrackParam out of the AliESDtracks
535 posV0track = new AliExternalTrackParam(*posVV0track);
536 negV0track = new AliExternalTrackParam(*negVV0track);
538 // Define the AODv0 from ESDv0 if reading ESDs
539 v0 = TransformESDv0toAODv0(esdV0,twoTrackArrayV0);
541 if( !posV0track || !negV0track ){
542 AliDebug(1,Form(" Couldn't get the V0 daughters"));
546 // fill in the v0 two-external-track-param array
547 twoTrackArrayV0->AddAt(posV0track,0);
548 twoTrackArrayV0->AddAt(negV0track,1);
551 dcaV0 = v0->DcaV0Daughters();
553 // Define the V0 (neutral) track
554 AliNeutralTrackParam *trackV0;
556 const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
557 if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
559 Double_t xyz[3], pxpypz[3];
561 esdV0->PxPyPz(pxpypz);
562 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
563 trackV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
565 // Fill in the object array to create the cascade
566 twoTrackArrayCasc->AddAt(postrack1,0);
567 twoTrackArrayCasc->AddAt(trackV0,1);
568 // Compute the cascade vertex
569 AliAODVertex *vertexCasc = 0;
570 if(fFindVertexForCascades) {
571 // DCA between the two tracks
572 dcaCasc = postrack1->GetDCA(trackV0,fBzkG,xdummy,ydummy);
574 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
576 // assume Cascade decays at the primary vertex
577 Double_t pos[3],cov[6],chi2perNDF;
579 fV1->GetCovMatrix(cov);
580 chi2perNDF = fV1->GetChi2toNDF();
581 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
585 delete posV0track; posV0track=NULL;
586 delete negV0track; negV0track=NULL;
587 delete trackV0; trackV0=NULL;
588 if(!fInputAOD) {delete v0; v0=NULL;}
589 twoTrackArrayV0->Clear();
590 twoTrackArrayCasc->Clear();
594 // Create and store the Cascade if passed the cuts
595 ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,v0,dcaCasc,okCascades);
596 if(okCascades && ioCascade) {
597 //AliDebug(1,Form("Storing a cascade object... "));
598 // add the vertex and the cascade to the AOD
599 AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
600 rc = new(aodCascadesRef[iCascades++])AliAODRecoCascadeHF(*ioCascade);
601 rc->SetSecondaryVtx(vCasc);
602 vCasc->SetParent(rc);
603 rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
604 if(!fInputAOD) vCasc->AddDaughter(v0); // just to fill ref #0 ??
605 AddRefs(vCasc,rc,event,twoTrackArrayCasc); // add the track (proton)
606 vCasc->AddDaughter(v0); // fill the 2prong V0
610 delete posV0track; posV0track=NULL;
611 delete negV0track; negV0track=NULL;
612 delete trackV0; trackV0=NULL;
613 twoTrackArrayV0->Clear();
614 twoTrackArrayCasc->Clear();
615 if(ioCascade) { delete ioCascade; ioCascade=NULL; }
616 if(vertexCasc) { delete vertexCasc; vertexCasc=NULL; }
617 if(!fInputAOD) {delete v0; v0=NULL;}
619 } // end loop on V0's
622 // If there is less than 2 particles continue
624 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
628 if(postrack1->Charge()<0 && !fLikeSign) continue;
630 // LOOP ON NEGATIVE TRACKS
631 for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
633 //if(iTrkN1%1==0) AliDebug(1,Form(" 1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));
634 //if(iTrkN1%1==0) printf(" 1st loop on neg: track number %d of %d\n",iTrkN1,nSeleTrks);
636 if(iTrkN1==iTrkP1) continue;
638 // get track from tracks array
639 negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
641 if(negtrack1->Charge()>0 && !fLikeSign) continue;
643 if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
646 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
649 if(postrack1->Charge()==negtrack1->Charge()) { // like-sign
650 isLikeSign2Prong=kTRUE;
651 if(!fLikeSign) continue;
652 if(iTrkN1<iTrkP1) continue; // this is needed to avoid double-counting of like-sign
653 } else { // unlike-sign
654 isLikeSign2Prong=kFALSE;
655 if(postrack1->Charge()<0 || negtrack1->Charge()>0) continue; // this is needed to avoid double-counting of unlike-sign
657 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
662 // back to primary vertex
663 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
664 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
665 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
666 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
668 // DCA between the two tracks
669 dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
670 if(dcap1n1>dcaMax) { negtrack1=0; continue; }
673 twoTrackArray1->AddAt(postrack1,0);
674 twoTrackArray1->AddAt(negtrack1,1);
675 AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
677 twoTrackArray1->Clear();
683 if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) {
685 io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
687 if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI) || (isLikeSign2Prong && (okD0 || okJPSI))) {
688 // add the vertex and the decay to the AOD
689 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
690 if(!isLikeSign2Prong) {
692 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
693 rd->SetSecondaryVtx(v2Prong);
694 v2Prong->SetParent(rd);
695 AddRefs(v2Prong,rd,event,twoTrackArray1);
698 rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
699 rd->SetSecondaryVtx(v2Prong);
700 if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
701 AddRefs(v2Prong,rd,event,twoTrackArray1);
703 } else { // isLikeSign2Prong
704 rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
705 rd->SetSecondaryVtx(v2Prong);
706 v2Prong->SetParent(rd);
707 AddRefs(v2Prong,rd,event,twoTrackArray1);
709 // Set selection bit for PID
710 if(okD0) SetSelectionBitForPID(fCutsD0toKpi,rd,AliRDHFCuts::kD0toKpiPID);
713 if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
714 // write references in io2Prong
716 AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
718 vertexp1n1->AddDaughter(postrack1);
719 vertexp1n1->AddDaughter(negtrack1);
721 io2Prong->SetSecondaryVtx(vertexp1n1);
722 //printf("---> %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
723 // create a track from the D0
724 AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
726 // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
727 for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
729 if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
731 if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
734 if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] ||
735 evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] ||
736 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
739 //if(iTrkSoftPi%1==0) AliDebug(1,Form(" 1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));
741 trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
742 if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
744 // get track from tracks array
745 trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
746 // trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
747 SetParametersAtVertex(trackPi,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkSoftPi));
748 twoTrackArrayCasc->AddAt(trackPi,0);
749 twoTrackArrayCasc->AddAt(trackD0,1);
751 AliAODVertex *vertexCasc = 0;
753 if(fFindVertexForDstar) {
754 // DCA between the two tracks
755 dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
757 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
759 // assume Dstar decays at the primary vertex
760 Double_t pos[3],cov[6],chi2perNDF;
762 fV1->GetCovMatrix(cov);
763 chi2perNDF = fV1->GetChi2toNDF();
764 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
768 twoTrackArrayCasc->Clear();
773 ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
775 // add the D0 to the AOD (if not already done)
777 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
778 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
779 rd->SetSecondaryVtx(v2Prong);
780 v2Prong->SetParent(rd);
781 AddRefs(v2Prong,rd,event,twoTrackArray1);
782 okD0=kTRUE; // this is done to add it only once
784 // add the vertex and the cascade to the AOD
785 AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
786 rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
787 rc->SetSecondaryVtx(vCasc);
788 vCasc->SetParent(rc);
789 if(!fInputAOD) vCasc->AddDaughter(rd); // just to fill ref #0
790 AddRefs(vCasc,rc,event,twoTrackArrayCasc);
791 vCasc->AddDaughter(rd); // add the D0 (in ref #1)
792 // Set selection bit for PID
793 SetSelectionBitForPID(fCutsDStartoKpipi,rc,AliRDHFCuts::kDstarPID);
795 twoTrackArrayCasc->Clear();
797 if(ioCascade) {delete ioCascade; ioCascade=NULL;}
798 delete vertexCasc; vertexCasc=NULL;
799 } // end loop on soft pi tracks
801 if(trackD0) {delete trackD0; trackD0=NULL;}
804 if(io2Prong) {delete io2Prong; io2Prong=NULL;}
807 twoTrackArray1->Clear();
808 if( (!f3Prong && !f4Prong) ||
809 (isLikeSign2Prong && !f3Prong) ) {
816 // 2nd LOOP ON POSITIVE TRACKS
817 for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
819 if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
821 //if(iTrkP2%1==0) AliDebug(1,Form(" 2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));
823 // get track from tracks array
824 postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
826 if(postrack2->Charge()<0) continue;
828 if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
831 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
832 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
833 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
836 if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
837 if(!fLikeSign3prong) continue;
838 if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
839 isLikeSign3Prong=kTRUE;
843 } else { // normal triplet (+-+)
844 isLikeSign3Prong=kFALSE;
846 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
847 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
848 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
852 // back to primary vertex
853 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
854 // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
855 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
856 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
857 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
858 SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
860 //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
862 dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
863 if(dcap2n1>dcaMax) { postrack2=0; continue; }
864 dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
865 if(dcap1p2>dcaMax) { postrack2=0; continue; }
867 // check invariant mass cuts for D+,Ds,Lc
870 if(postrack2->Charge()>0) {
871 threeTrackArray->AddAt(postrack1,0);
872 threeTrackArray->AddAt(negtrack1,1);
873 threeTrackArray->AddAt(postrack2,2);
875 threeTrackArray->AddAt(negtrack1,0);
876 threeTrackArray->AddAt(postrack1,1);
877 threeTrackArray->AddAt(postrack2,2);
879 if(fMassCutBeforeVertexing)
880 massCutOK = SelectInvMassAndPt(threeTrackArray);
883 if(f3Prong && !massCutOK) {
884 threeTrackArray->Clear();
892 twoTrackArray2->AddAt(postrack2,0);
893 twoTrackArray2->AddAt(negtrack1,1);
894 AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
896 twoTrackArray2->Clear();
901 // 3 prong candidates
902 if(f3Prong && massCutOK) {
904 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
905 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
907 AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
908 if(!isLikeSign3Prong) {
909 rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
910 rd->SetSecondaryVtx(v3Prong);
911 v3Prong->SetParent(rd);
912 AddRefs(v3Prong,rd,event,threeTrackArray);
913 } else { // isLikeSign3Prong
915 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
916 rd->SetSecondaryVtx(v3Prong);
917 v3Prong->SetParent(rd);
918 AddRefs(v3Prong,rd,event,threeTrackArray);
921 // Set selection bit for PID
922 SetSelectionBitForPID(fCutsDplustoKpipi,rd,AliRDHFCuts::kDplusPID);
923 SetSelectionBitForPID(fCutsDstoKKpi,rd,AliRDHFCuts::kDsPID);
924 SetSelectionBitForPID(fCutsLctopKpi,rd,AliRDHFCuts::kLcPID);
926 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
927 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
930 // 4 prong candidates
932 // don't make 4 prong with like-sign pairs and triplets
933 && !isLikeSign2Prong && !isLikeSign3Prong
934 // track-to-track dca cuts already now
935 && dcap1n1 < fCutsD0toKpipipi->GetDCACut()
936 && dcap2n1 < fCutsD0toKpipipi->GetDCACut()) {
938 // back to primary vertex
939 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
940 // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
941 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
942 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
943 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
944 SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
946 // Vertexing for these 3 (can be taken from above?)
947 threeTrackArray->AddAt(postrack1,0);
948 threeTrackArray->AddAt(negtrack1,1);
949 threeTrackArray->AddAt(postrack2,2);
950 AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
952 // 3rd LOOP ON NEGATIVE TRACKS (for 4 prong)
953 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
955 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
957 //if(iTrkN2%1==0) AliDebug(1,Form(" 3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
959 // get track from tracks array
960 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
962 if(negtrack2->Charge()>0) continue;
964 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
966 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
967 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
968 evtNumber[iTrkP2]==evtNumber[iTrkN2] ||
969 evtNumber[iTrkP1]==evtNumber[iTrkN1] ||
970 evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
971 evtNumber[iTrkN1]==evtNumber[iTrkP2]) continue;
974 // back to primary vertex
975 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
976 // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
977 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
978 // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
979 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
980 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
981 SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
982 SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
984 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
985 if(dcap1n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
986 dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
987 if(dcap2n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
990 fourTrackArray->AddAt(postrack1,0);
991 fourTrackArray->AddAt(negtrack1,1);
992 fourTrackArray->AddAt(postrack2,2);
993 fourTrackArray->AddAt(negtrack2,3);
995 // check invariant mass cuts for D0
997 if(fMassCutBeforeVertexing)
998 massCutOK = SelectInvMassAndPt(fourTrackArray);
1001 fourTrackArray->Clear();
1007 AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
1008 io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
1010 AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
1011 rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
1012 rd->SetSecondaryVtx(v4Prong);
1013 v4Prong->SetParent(rd);
1014 AddRefs(v4Prong,rd,event,fourTrackArray);
1017 if(io4Prong) {delete io4Prong; io4Prong=NULL;}
1018 if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;}
1019 fourTrackArray->Clear();
1022 } // end loop on negative tracks
1024 threeTrackArray->Clear();
1025 delete vertexp1n1p2;
1032 } // end 2nd loop on positive tracks
1034 twoTrackArray2->Clear();
1036 // 2nd LOOP ON NEGATIVE TRACKS (for 3 prong -+-)
1037 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
1039 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1041 //if(iTrkN2%1==0) AliDebug(1,Form(" 2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
1043 // get track from tracks array
1044 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1046 if(negtrack2->Charge()>0) continue;
1048 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1051 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
1052 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
1053 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1056 if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
1057 if(!fLikeSign3prong) continue;
1058 if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
1059 isLikeSign3Prong=kTRUE;
1063 } else { // normal triplet (-+-)
1064 isLikeSign3Prong=kFALSE;
1067 // back to primary vertex
1068 // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1069 // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1070 // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1071 SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
1072 SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
1073 SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
1074 //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
1076 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1077 if(dcap1n2>dcaMax) { negtrack2=0; continue; }
1078 dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1079 if(dcan1n2>dcaMax) { negtrack2=0; continue; }
1081 threeTrackArray->AddAt(negtrack1,0);
1082 threeTrackArray->AddAt(postrack1,1);
1083 threeTrackArray->AddAt(negtrack2,2);
1085 // check invariant mass cuts for D+,Ds,Lc
1087 if(fMassCutBeforeVertexing && f3Prong)
1088 massCutOK = SelectInvMassAndPt(threeTrackArray);
1091 threeTrackArray->Clear();
1097 twoTrackArray2->AddAt(postrack1,0);
1098 twoTrackArray2->AddAt(negtrack2,1);
1100 AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1102 twoTrackArray2->Clear();
1108 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1109 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
1111 AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1112 if(!isLikeSign3Prong) {
1113 rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1114 rd->SetSecondaryVtx(v3Prong);
1115 v3Prong->SetParent(rd);
1116 AddRefs(v3Prong,rd,event,threeTrackArray);
1117 } else { // isLikeSign3Prong
1118 if(fLikeSign3prong){
1119 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1120 rd->SetSecondaryVtx(v3Prong);
1121 v3Prong->SetParent(rd);
1122 AddRefs(v3Prong,rd,event,threeTrackArray);
1125 // Set selection bit for PID
1126 SetSelectionBitForPID(fCutsDplustoKpipi,rd,AliRDHFCuts::kDplusPID);
1127 SetSelectionBitForPID(fCutsDstoKKpi,rd,AliRDHFCuts::kDsPID);
1128 SetSelectionBitForPID(fCutsLctopKpi,rd,AliRDHFCuts::kLcPID);
1130 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1131 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1133 threeTrackArray->Clear();
1137 } // end 2nd loop on negative tracks
1139 twoTrackArray2->Clear();
1143 } // end 1st loop on negative tracks
1146 } // end 1st loop on positive tracks
1149 AliDebug(1,Form(" Total HF vertices in event = %d;",
1150 (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
1152 AliDebug(1,Form(" D0->Kpi in event = %d;",
1153 (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
1156 AliDebug(1,Form(" JPSI->ee in event = %d;",
1157 (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
1160 AliDebug(1,Form(" Charm->3Prong in event = %d;",
1161 (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
1164 AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
1165 (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
1168 AliDebug(1,Form(" D*->D0pi in event = %d;\n",
1169 (Int_t)aodDstarTClArr->GetEntriesFast()));
1172 AliDebug(1,Form(" cascades -> v0 + track in event = %d;\n",
1173 (Int_t)aodCascadesTClArr->GetEntriesFast()));
1176 AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
1177 (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
1179 if(fLikeSign3prong && f3Prong) {
1180 AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
1181 (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
1185 twoTrackArray1->Delete(); delete twoTrackArray1;
1186 twoTrackArray2->Delete(); delete twoTrackArray2;
1187 twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1188 twoTrackArrayV0->Delete(); delete twoTrackArrayV0;
1189 threeTrackArray->Clear();
1190 threeTrackArray->Delete(); delete threeTrackArray;
1191 fourTrackArray->Delete(); delete fourTrackArray;
1192 delete [] seleFlags; seleFlags=NULL;
1193 if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
1194 tracksAtVertex.Delete();
1197 seleTrksArray.Delete();
1198 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
1202 //printf("Trks: total %d sele %d\n",fnTrksTotal,fnSeleTrksTotal);
1206 //----------------------------------------------------------------------------
1207 void AliAnalysisVertexingHF::AddRefs(AliAODVertex *v,AliAODRecoDecayHF *rd,
1208 const AliVEvent *event,
1209 const TObjArray *trkArray) const
1211 // Add the AOD tracks as daughters of the vertex (TRef)
1212 // Also add the references to the primary vertex and to the cuts
1215 AddDaughterRefs(v,event,trkArray);
1216 rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1220 rd->SetListOfCutsRef((TList*)fListOfCuts);
1221 //fListOfCuts->Print();
1222 cout<<fListOfCuts<<endl;
1223 TList *l=(TList*)rd->GetListOfCuts();
1225 if(l) {l->Print(); }else{printf("error\n");}
1230 //----------------------------------------------------------------------------
1231 void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,
1232 const AliVEvent *event,
1233 const TObjArray *trkArray) const
1235 // Add the AOD tracks as daughters of the vertex (TRef)
1237 Int_t nDg = v->GetNDaughters();
1239 if(nDg) dg = v->GetDaughter(0);
1241 if(dg) return; // daughters already added
1243 Int_t nTrks = trkArray->GetEntriesFast();
1245 AliExternalTrackParam *track = 0;
1246 AliAODTrack *aodTrack = 0;
1249 for(Int_t i=0; i<nTrks; i++) {
1250 track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1251 id = (Int_t)track->GetID();
1252 //printf("---> %d\n",id);
1253 if(id<0) continue; // this track is a AliAODRecoDecay
1254 aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
1255 v->AddDaughter(aodTrack);
1260 //---------------------------------------------------------------------------
1261 void AliAnalysisVertexingHF::FixReferences(AliAODEvent *aod)
1263 // Checks that the references to the daughter tracks are properly
1264 // assigned and reassigns them if needed
1268 TClonesArray *inputArray=(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
1269 if(!inputArray) return;
1271 AliAODTrack *track = 0;
1272 AliAODVertex *vertex = 0;
1274 Bool_t needtofix=kFALSE;
1275 for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1276 vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1277 for(Int_t id=0; id<vertex->GetNDaughters(); id++) {
1278 track = (AliAODTrack*)vertex->GetDaughter(id);
1279 if(!track->GetStatus()) needtofix=kTRUE;
1281 if(needtofix) break;
1284 if(!needtofix) return;
1287 printf("Fixing references\n");
1289 fAODMapSize = 100000;
1290 fAODMap = new Int_t[fAODMapSize];
1292 for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
1293 track = aod->GetTrack(i);
1295 // skip pure ITS SA tracks
1296 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
1298 // skip tracks without ITS
1299 if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
1301 // TEMPORARY: check that the cov matrix is there
1302 Double_t covtest[21];
1303 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1306 fAODMap[(Int_t)track->GetID()] = i;
1310 Int_t ids[4]={-1,-1,-1,-1};
1311 for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
1312 Bool_t cascade=kFALSE;
1313 vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
1315 Int_t nDgs = vertex->GetNDaughters();
1316 for(id=0; id<nDgs; id++) {
1317 track = (AliAODTrack*)vertex->GetDaughter(id);
1318 if(track->Charge()==0) {cascade=kTRUE; continue;} // cascade
1319 ids[id]=(Int_t)track->GetID();
1320 vertex->RemoveDaughter(track);
1322 if(cascade) continue;
1323 for(id=0; id<nDgs; id++) {
1324 track = aod->GetTrack(fAODMap[ids[id]]);
1325 vertex->AddDaughter(track);
1332 //----------------------------------------------------------------------------
1333 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1334 TObjArray *twoTrackArray,AliVEvent *event,
1335 AliAODVertex *secVert,
1336 AliAODRecoDecayHF2Prong *rd2Prong,
1338 Bool_t &okDstar) const
1340 // Make the cascade as a 2Prong decay and check if it passes Dstar
1341 // reconstruction cuts
1345 Bool_t dummy1,dummy2,dummy3;
1347 // We use Make2Prong to construct the AliAODRecoCascadeHF
1348 // (which inherits from AliAODRecoDecayHF2Prong)
1349 AliAODRecoCascadeHF *theCascade =
1350 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1351 dummy1,dummy2,dummy3);
1352 if(!theCascade) return 0x0;
1355 AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1356 theCascade->SetCharge(trackPi->Charge());
1358 //--- selection cuts
1360 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1362 Int_t idSoftPi=(Int_t)trackPi->GetID();
1363 AliAODTrack* trackPiAOD=(AliAODTrack*)event->GetTrack(fAODMap[idSoftPi]);
1364 tmpCascade->GetSecondaryVtx()->AddDaughter(trackPiAOD);
1366 tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
1368 tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
1370 AliAODVertex *primVertexAOD=0;
1371 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1372 // take event primary vertex
1373 primVertexAOD = PrimaryVertex();
1374 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1375 rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
1379 okDstar = (Bool_t)fCutsDStartoKpipi->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1380 if(okDstar) theCascade->SetSelectionBit(AliRDHFCuts::kDstarCuts);
1382 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1383 tmpCascade->UnsetOwnPrimaryVtx();
1384 delete tmpCascade; tmpCascade=NULL;
1385 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1386 rd2Prong->UnsetOwnPrimaryVtx();
1388 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1396 //----------------------------------------------------------------------------
1397 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1398 TObjArray *twoTrackArray,AliVEvent *event,
1399 AliAODVertex *secVert,
1402 Bool_t &okCascades) const
1405 // Make the cascade as a 2Prong decay and check if it passes
1406 // cascades reconstruction cuts
1408 // AliDebug(2,Form(" building the cascade"));
1410 Bool_t dummy1,dummy2,dummy3;
1412 // We use Make2Prong to construct the AliAODRecoCascadeHF
1413 // (which inherits from AliAODRecoDecayHF2Prong)
1414 AliAODRecoCascadeHF *theCascade =
1415 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1416 dummy1,dummy2,dummy3);
1417 if(!theCascade) return 0x0;
1419 // bachelor track and charge
1420 AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1421 theCascade->SetCharge(trackBachelor->Charge());
1423 //--- selection cuts
1426 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1428 Int_t idBachelor=(Int_t)trackBachelor->GetID();
1429 AliAODTrack* trackBachelorAOD=(AliAODTrack*)event->GetTrack(fAODMap[idBachelor]);
1430 tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelorAOD);
1432 tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
1434 tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
1436 AliAODVertex *primVertexAOD=0;
1437 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1438 // take event primary vertex
1439 primVertexAOD = PrimaryVertex();
1440 if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
1441 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1445 if(fCascades && fInputAOD){
1446 okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1449 //AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied"));
1451 }// no cuts implemented from ESDs
1452 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1453 tmpCascade->UnsetOwnPrimaryVtx();
1454 delete tmpCascade; tmpCascade=NULL;
1455 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1461 //-----------------------------------------------------------------------------
1462 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
1463 TObjArray *twoTrackArray,AliVEvent *event,
1464 AliAODVertex *secVert,Double_t dca,
1465 Bool_t &okD0,Bool_t &okJPSI,
1466 Bool_t &okD0fromDstar) const
1468 // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
1469 // reconstruction cuts
1470 // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
1472 okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
1474 Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1475 AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1476 AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
1478 // propagate tracks to secondary vertex, to compute inv. mass
1479 postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1480 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1482 Double_t momentum[3];
1483 postrack->GetPxPyPz(momentum);
1484 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1485 negtrack->GetPxPyPz(momentum);
1486 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1489 // invariant mass cut (try to improve coding here..)
1490 Bool_t okMassCut=kFALSE;
1491 if(!okMassCut && fD0toKpi) if(SelectInvMassAndPt(0,2,px,py,pz)) okMassCut=kTRUE;
1492 if(!okMassCut && fJPSItoEle) if(SelectInvMassAndPt(1,2,px,py,pz)) okMassCut=kTRUE;
1493 if(!okMassCut && fDstar) if(SelectInvMassAndPt(3,2,px,py,pz)) okMassCut=kTRUE;
1494 if(!okMassCut && fCascades) if(SelectInvMassAndPt(5,2,px,py,pz)) okMassCut=kTRUE;
1496 //AliDebug(2," candidate didn't pass mass cut");
1499 // primary vertex to be used by this candidate
1500 AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
1501 if(!primVertexAOD) return 0x0;
1504 Double_t d0z0[2],covd0z0[3];
1505 postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1507 d0err[0] = TMath::Sqrt(covd0z0[0]);
1508 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1510 d0err[1] = TMath::Sqrt(covd0z0[0]);
1512 // create the object AliAODRecoDecayHF2Prong
1513 AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
1514 the2Prong->SetOwnPrimaryVtx(primVertexAOD);
1515 UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1516 the2Prong->SetProngIDs(2,id);
1517 delete primVertexAOD; primVertexAOD=NULL;
1520 if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
1522 // Add daughter references already here
1523 if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
1527 okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event);
1528 if(okD0) the2Prong->SetSelectionBit(AliRDHFCuts::kD0toKpiCuts);
1530 //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
1531 // select J/psi from B
1533 okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1535 //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
1536 // select D0->Kpi from Dstar
1538 okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
1539 if(okD0fromDstar) the2Prong->SetSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
1541 //if(fDebug && fDstar) printf(" %d\n",(Int_t)okD0fromDstar);
1545 // remove the primary vertex (was used only for selection)
1546 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1547 the2Prong->UnsetOwnPrimaryVtx();
1550 // get PID info from ESD
1551 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1552 if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1553 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1554 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1555 Double_t esdpid[10];
1556 for(Int_t i=0;i<5;i++) {
1557 esdpid[i] = esdpid0[i];
1558 esdpid[5+i] = esdpid1[i];
1560 the2Prong->SetPID(2,esdpid);
1564 //----------------------------------------------------------------------------
1565 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
1566 TObjArray *threeTrackArray,AliVEvent *event,
1567 AliAODVertex *secVert,Double_t dispersion,
1568 const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
1569 Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1570 Bool_t &ok3Prong) const
1572 // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1573 // reconstruction cuts
1578 if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
1580 Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1581 Double_t momentum[3];
1584 AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1585 AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1586 AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1588 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1589 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1590 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1591 postrack1->GetPxPyPz(momentum);
1592 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1593 negtrack->GetPxPyPz(momentum);
1594 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1595 postrack2->GetPxPyPz(momentum);
1596 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1598 // invariant mass cut for D+, Ds, Lc
1599 Bool_t okMassCut=kFALSE;
1600 if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1601 if(!okMassCut && f3Prong) if(SelectInvMassAndPt(2,3,px,py,pz)) okMassCut=kTRUE;
1603 //AliDebug(2," candidate didn't pass mass cut");
1607 // primary vertex to be used by this candidate
1608 AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
1609 if(!primVertexAOD) return 0x0;
1611 Double_t d0z0[2],covd0z0[3];
1612 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1614 d0err[0] = TMath::Sqrt(covd0z0[0]);
1615 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1617 d0err[1] = TMath::Sqrt(covd0z0[0]);
1618 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1620 d0err[2] = TMath::Sqrt(covd0z0[0]);
1623 // create the object AliAODRecoDecayHF3Prong
1624 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1625 Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1626 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]));
1627 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]));
1628 Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
1629 AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1630 the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1631 UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1632 the3Prong->SetProngIDs(3,id);
1634 delete primVertexAOD; primVertexAOD=NULL;
1636 // Add daughter references already here
1637 if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
1639 // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1643 if(fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1645 the3Prong->SetSelectionBit(AliRDHFCuts::kDplusCuts);
1647 if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1649 the3Prong->SetSelectionBit(AliRDHFCuts::kDsCuts);
1651 if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
1653 the3Prong->SetSelectionBit(AliRDHFCuts::kLcCuts);
1656 //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1658 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1659 the3Prong->UnsetOwnPrimaryVtx();
1662 // get PID info from ESD
1663 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1664 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1665 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1666 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1667 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1668 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1670 Double_t esdpid[15];
1671 for(Int_t i=0;i<5;i++) {
1672 esdpid[i] = esdpid0[i];
1673 esdpid[5+i] = esdpid1[i];
1674 esdpid[10+i] = esdpid2[i];
1676 the3Prong->SetPID(3,esdpid);
1680 //----------------------------------------------------------------------------
1681 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1682 TObjArray *fourTrackArray,AliVEvent *event,
1683 AliAODVertex *secVert,
1684 const AliAODVertex *vertexp1n1,
1685 const AliAODVertex *vertexp1n1p2,
1686 Double_t dcap1n1,Double_t dcap1n2,
1687 Double_t dcap2n1,Double_t dcap2n2,
1688 Bool_t &ok4Prong) const
1690 // Make 4Prong candidates and check if they pass D0toKpipipi
1691 // reconstruction cuts
1692 // G.E.Bruno, R.Romita
1695 if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
1697 Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1699 AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1700 AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1701 AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1702 AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1704 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1705 negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1706 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1707 negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1709 Double_t momentum[3];
1710 postrack1->GetPxPyPz(momentum);
1711 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1712 negtrack1->GetPxPyPz(momentum);
1713 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1714 postrack2->GetPxPyPz(momentum);
1715 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1716 negtrack2->GetPxPyPz(momentum);
1717 px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1719 // invariant mass cut for rho or D0 (try to improve coding here..)
1720 Bool_t okMassCut=kFALSE;
1721 if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1722 if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) { //no PID, to be implemented with PID
1723 if(SelectInvMassAndPt(4,4,px,py,pz)) okMassCut=kTRUE;
1726 //if(fDebug) printf(" candidate didn't pass mass cut\n");
1727 //printf(" candidate didn't pass mass cut\n");
1731 // primary vertex to be used by this candidate
1732 AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
1733 if(!primVertexAOD) return 0x0;
1735 Double_t d0z0[2],covd0z0[3];
1736 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1738 d0err[0] = TMath::Sqrt(covd0z0[0]);
1739 negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1741 d0err[1] = TMath::Sqrt(covd0z0[0]);
1742 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1744 d0err[2] = TMath::Sqrt(covd0z0[0]);
1745 negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1747 d0err[3] = TMath::Sqrt(covd0z0[0]);
1750 // create the object AliAODRecoDecayHF4Prong
1751 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1752 Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
1753 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]));
1754 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]));
1755 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]));
1757 AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
1758 the4Prong->SetOwnPrimaryVtx(primVertexAOD);
1759 UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
1760 the4Prong->SetProngIDs(4,id);
1762 delete primVertexAOD; primVertexAOD=NULL;
1764 ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
1767 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1768 the4Prong->UnsetOwnPrimaryVtx();
1772 // get PID info from ESD
1773 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1774 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1775 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1776 if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
1777 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1778 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1779 Double_t esdpid3[5]={0.,0.,0.,0.,0.};
1780 if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
1782 Double_t esdpid[20];
1783 for(Int_t i=0;i<5;i++) {
1784 esdpid[i] = esdpid0[i];
1785 esdpid[5+i] = esdpid1[i];
1786 esdpid[10+i] = esdpid2[i];
1787 esdpid[15+i] = esdpid3[i];
1789 the4Prong->SetPID(4,esdpid);
1793 //-----------------------------------------------------------------------------
1794 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
1795 AliVEvent *event) const
1797 // Returns primary vertex to be used for this candidate
1799 AliESDVertex *vertexESD = 0;
1800 AliAODVertex *vertexAOD = 0;
1803 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1804 // primary vertex from the input event
1806 vertexESD = new AliESDVertex(*fV1);
1809 // primary vertex specific to this candidate
1811 Int_t nTrks = trkArray->GetEntriesFast();
1812 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1814 if(fRecoPrimVtxSkippingTrks) {
1815 // recalculating the vertex
1817 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1818 Float_t diamondcovxy[3];
1819 event->GetDiamondCovXY(diamondcovxy);
1820 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1821 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1822 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1823 vertexer->SetVtxStart(diamond);
1824 delete diamond; diamond=NULL;
1825 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1826 vertexer->SetOnlyFitter();
1828 Int_t skipped[1000];
1829 Int_t nTrksToSkip=0,id;
1830 AliExternalTrackParam *t = 0;
1831 for(Int_t i=0; i<nTrks; i++) {
1832 t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1833 id = (Int_t)t->GetID();
1835 skipped[nTrksToSkip++] = id;
1838 // For AOD, skip also tracks without covariance matrix
1840 Double_t covtest[21];
1841 for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1842 AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1843 if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1844 id = (Int_t)vtrack->GetID();
1846 skipped[nTrksToSkip++] = id;
1850 for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
1852 vertexer->SetSkipTracks(nTrksToSkip,skipped);
1853 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1855 } else if(fRmTrksFromPrimVtx && nTrks>0) {
1856 // removing the prongs tracks
1858 TObjArray rmArray(nTrks);
1859 UShort_t *rmId = new UShort_t[nTrks];
1860 AliESDtrack *esdTrack = 0;
1862 for(Int_t i=0; i<nTrks; i++) {
1863 t = (AliESDtrack*)trkArray->UncheckedAt(i);
1864 esdTrack = new AliESDtrack(*t);
1865 rmArray.AddLast(esdTrack);
1866 if(esdTrack->GetID()>=0) {
1867 rmId[i]=(UShort_t)esdTrack->GetID();
1872 Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
1873 vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1874 delete [] rmId; rmId=NULL;
1879 if(!vertexESD) return vertexAOD;
1880 if(vertexESD->GetNContributors()<=0) {
1881 //AliDebug(2,"vertexing failed");
1882 delete vertexESD; vertexESD=NULL;
1886 delete vertexer; vertexer=NULL;
1890 // convert to AliAODVertex
1891 Double_t pos[3],cov[6],chi2perNDF;
1892 vertexESD->GetXYZ(pos); // position
1893 vertexESD->GetCovMatrix(cov); //covariance matrix
1894 chi2perNDF = vertexESD->GetChi2toNDF();
1895 delete vertexESD; vertexESD=NULL;
1897 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1901 //-----------------------------------------------------------------------------
1902 void AliAnalysisVertexingHF::PrintStatus() const {
1903 // Print parameters being used
1905 //printf("Preselections:\n");
1906 // fTrackFilter->Dump();
1908 printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1910 printf("Secondary vertex with AliVertexerTracks\n");
1912 if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1913 if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1915 printf("Reconstruct D0->Kpi candidates with cuts:\n");
1916 if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
1919 printf("Reconstruct D*->D0pi candidates with cuts:\n");
1920 if(fFindVertexForDstar) {
1921 printf(" Reconstruct a secondary vertex for the D*\n");
1923 printf(" Assume the D* comes from the primary vertex\n");
1925 if(fCutsDStartoKpipi) fCutsDStartoKpipi->PrintAll();
1928 printf("Reconstruct J/psi from B candidates with cuts:\n");
1929 if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
1932 printf("Reconstruct 3 prong candidates.\n");
1933 printf(" D+->Kpipi cuts:\n");
1934 if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
1935 printf(" Ds->KKpi cuts:\n");
1936 if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
1937 printf(" Lc->pKpi cuts:\n");
1938 if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
1941 printf("Reconstruct 4 prong candidates.\n");
1942 printf(" D0->Kpipipi cuts:\n");
1943 if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
1946 printf("Reconstruct cascades candidates formed with v0s.\n");
1947 printf(" Lc -> k0s P & Lc -> L Pi cuts:\n");
1948 if(fCutsLctoV0) fCutsLctoV0->PrintAll();
1953 //-----------------------------------------------------------------------------
1954 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
1955 Double_t &dispersion,Bool_t useTRefArray) const
1957 // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1959 AliESDVertex *vertexESD = 0;
1960 AliAODVertex *vertexAOD = 0;
1962 if(!fSecVtxWithKF) { // AliVertexerTracks
1964 AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1965 vertexer->SetVtxStart(fV1);
1966 vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1967 delete vertexer; vertexer=NULL;
1969 if(!vertexESD) return vertexAOD;
1971 if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
1972 //AliDebug(2,"vertexing failed");
1973 delete vertexESD; vertexESD=NULL;
1977 Double_t vertRadius2=vertexESD->GetXv()*vertexESD->GetXv()+vertexESD->GetYv()*vertexESD->GetYv();
1979 // vertex outside beam pipe, reject candidate to avoid propagation through material
1980 delete vertexESD; vertexESD=NULL;
1984 } else { // Kalman Filter vertexer (AliKFParticle)
1986 AliKFParticle::SetField(fBzkG);
1988 AliKFVertex vertexKF;
1990 Int_t nTrks = trkArray->GetEntriesFast();
1991 for(Int_t i=0; i<nTrks; i++) {
1992 AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1993 AliKFParticle daughterKF(*esdTrack,211);
1994 vertexKF.AddDaughter(daughterKF);
1996 vertexESD = new AliESDVertex(vertexKF.Parameters(),
1997 vertexKF.CovarianceMatrix(),
1999 vertexKF.GetNContributors());
2003 // convert to AliAODVertex
2004 Double_t pos[3],cov[6],chi2perNDF;
2005 vertexESD->GetXYZ(pos); // position
2006 vertexESD->GetCovMatrix(cov); //covariance matrix
2007 chi2perNDF = vertexESD->GetChi2toNDF();
2008 dispersion = vertexESD->GetDispersion();
2009 delete vertexESD; vertexESD=NULL;
2011 Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
2012 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
2016 //-----------------------------------------------------------------------------
2017 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt(TObjArray *trkArray) const {
2018 // Invariant mass cut on tracks
2020 Int_t nProngs=trkArray->GetEntriesFast();
2021 Int_t retval=kFALSE;
2023 Double_t momentum[3];
2024 Double_t px3[3],py3[3],pz3[3];
2025 Double_t px4[4],py4[4],pz4[4];
2026 AliESDtrack *postrack1=0;
2027 AliESDtrack *postrack2=0;
2028 AliESDtrack *negtrack1=0;
2029 AliESDtrack *negtrack2=0;
2033 // invariant mass cut for D+, Ds, Lc
2034 postrack1 = (AliESDtrack*)trkArray->UncheckedAt(0);
2035 negtrack1 = (AliESDtrack*)trkArray->UncheckedAt(1);
2036 postrack2 = (AliESDtrack*)trkArray->UncheckedAt(2);
2037 postrack1->GetPxPyPz(momentum);
2038 px3[0] = momentum[0]; py3[0] = momentum[1]; pz3[0] = momentum[2];
2039 negtrack1->GetPxPyPz(momentum);
2040 px3[1] = momentum[0]; py3[1] = momentum[1]; pz3[1] = momentum[2];
2041 postrack2->GetPxPyPz(momentum);
2042 px3[2] = momentum[0]; py3[2] = momentum[1]; pz3[2] = momentum[2];
2043 retval = SelectInvMassAndPt(2,3,px3,py3,pz3);
2046 // invariant mass cut for D0->4prong
2047 postrack1 = (AliESDtrack*)trkArray->UncheckedAt(0);
2048 negtrack1 = (AliESDtrack*)trkArray->UncheckedAt(1);
2049 postrack2 = (AliESDtrack*)trkArray->UncheckedAt(2);
2050 negtrack2 = (AliESDtrack*)trkArray->UncheckedAt(3);
2051 postrack1->GetPxPyPz(momentum);
2052 px4[0] = momentum[0]; py4[0] = momentum[1]; pz4[0] = momentum[2];
2053 negtrack1->GetPxPyPz(momentum);
2054 px4[1] = momentum[0]; py4[1] = momentum[1]; pz4[1] = momentum[2];
2055 postrack2->GetPxPyPz(momentum);
2056 px4[2] = momentum[0]; py4[2] = momentum[1]; pz4[2] = momentum[2];
2057 negtrack2->GetPxPyPz(momentum);
2058 px4[3] = momentum[0]; py4[3] = momentum[1]; pz4[3] = momentum[2];
2059 retval = SelectInvMassAndPt(4,4,px4,py4,pz4);
2062 printf("SelectInvMassAndPt(): wrong decay selection\n");
2068 //-----------------------------------------------------------------------------
2069 Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt(Int_t decay,
2073 Double_t *pz) const {
2074 // Check invariant mass cut and pt candidate cut
2076 UInt_t pdg2[2],pdg3[3],pdg4[4];
2080 Bool_t retval=kFALSE;
2084 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2086 minPt=fCutsD0toKpi->GetMinPtCandidate();
2088 if(fMassCalc2->Pt2() < minPt*minPt) break;
2090 pdg2[0]=211; pdg2[1]=321;
2091 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2092 minv = fMassCalc2->InvMass(nprongs,pdg2);
2093 if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
2094 pdg2[0]=321; pdg2[1]=211;
2095 minv = fMassCalc2->InvMass(nprongs,pdg2);
2096 if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
2099 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2101 minPt=fCutsJpsitoee->GetMinPtCandidate();
2103 if(fMassCalc2->Pt2() < minPt*minPt) break;
2105 pdg2[0]=11; pdg2[1]=11;
2106 mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
2107 minv = fMassCalc2->InvMass(nprongs,pdg2);
2108 if(TMath::Abs(minv-mPDG)<fCutsJpsitoee->GetMassCut()) retval=kTRUE;
2110 case 2: // D+->Kpipi
2111 fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
2113 minPt=TMath::Min(fCutsDplustoKpipi->GetMinPtCandidate(),fCutsDstoKKpi->GetMinPtCandidate());
2114 minPt=TMath::Min(minPt,fCutsLctopKpi->GetMinPtCandidate());
2116 if(fMassCalc3->Pt2() < minPt*minPt) break;
2118 pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
2119 mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
2120 minv = fMassCalc3->InvMass(nprongs,pdg3);
2121 if(TMath::Abs(minv-mPDG)<fCutsDplustoKpipi->GetMassCut()) retval=kTRUE;
2123 pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
2124 mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
2125 minv = fMassCalc3->InvMass(nprongs,pdg3);
2126 if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
2127 pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
2128 minv = fMassCalc3->InvMass(nprongs,pdg3);
2129 if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
2131 pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
2132 mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
2133 minv = fMassCalc3->InvMass(nprongs,pdg3);
2134 if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
2135 pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
2136 minv = fMassCalc3->InvMass(nprongs,pdg3);
2137 if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
2140 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2142 minPt=fCutsDStartoKpipi->GetMinPtCandidate();
2144 if(fMassCalc2->Pt2() < minPt*minPt) break;
2146 pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
2147 mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
2148 minv = fMassCalc2->InvMass(nprongs,pdg2);
2149 if(TMath::Abs(minv-mPDG)<fCutsDStartoKpipi->GetMassCut()) retval=kTRUE;
2151 case 4: // D0->Kpipipi without PID
2152 fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
2154 minPt=fCutsD0toKpipipi->GetMinPtCandidate();
2156 if(fMassCalc4->Pt2() < minPt*minPt) break;
2158 pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
2159 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2160 minv = fMassCalc4->InvMass(nprongs,pdg4);
2161 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
2162 pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
2163 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2164 minv = fMassCalc4->InvMass(nprongs,pdg4);
2165 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
2166 pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
2167 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2168 minv = fMassCalc4->InvMass(nprongs,pdg4);
2169 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
2170 pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
2171 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2172 minv = fMassCalc4->InvMass(nprongs,pdg4);
2173 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
2176 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
2177 mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
2178 pdg2[0]=2212;pdg2[1]=310;
2179 minv=fMassCalc2->InvMass(2,pdg2);
2180 if(TMath::Abs(minv-mPDG)<fCutsLctoV0->GetMassCut()) retval=kTRUE;
2181 pdg2[0]=211;pdg2[1]=3122;
2182 minv=fMassCalc2->InvMass(2,pdg2);
2183 if(TMath::Abs(minv-mPDG)<fCutsLctoV0->GetMassCut()) retval=kTRUE;
2186 printf("SelectInvMassAndPt(): wrong decay selection\n");
2192 //-----------------------------------------------------------------------------
2193 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
2195 TObjArray &seleTrksArray,TObjArray &tracksAtVertex,
2197 UChar_t *seleFlags,Int_t *evtNumber)
2199 // Apply single-track preselection.
2200 // Fill a TObjArray with selected tracks (for displaced vertices or
2201 // soft pion from D*). Selection flag stored in seleFlags.
2202 // Create the AliESDVertex object (convert from AliAODVertex if necessary)
2203 // In case of AOD input, also fill fAODMap for track index<->ID
2205 const AliVVertex *vprimary = event->GetPrimaryVertex();
2207 if(fV1) { delete fV1; fV1=NULL; }
2208 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
2211 UShort_t *indices = 0;
2212 Double_t pos[3],cov[6];
2214 if(!fInputAOD) { // ESD
2215 fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
2217 vprimary->GetXYZ(pos);
2218 vprimary->GetCovarianceMatrix(cov);
2219 fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2220 indices = new UShort_t[event->GetNumberOfTracks()];
2221 for(Int_t ijk=0; ijk<event->GetNumberOfTracks(); ijk++) indices[ijk]=0;
2222 fAODMapSize = 100000;
2223 fAODMap = new Int_t[fAODMapSize];
2227 Int_t entries = (Int_t)event->GetNumberOfTracks();
2228 Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
2231 // transfer ITS tracks from event to arrays
2232 for(Int_t i=0; i<entries; i++) {
2234 track = (AliVTrack*)event->GetTrack(i);
2236 // skip pure ITS SA tracks
2237 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2239 // skip tracks without ITS
2240 if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
2242 // TEMPORARY: check that the cov matrix is there
2243 Double_t covtest[21];
2244 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2248 AliAODTrack *aodt = (AliAODTrack*)track;
2249 if(aodt->GetUsedForPrimVtxFit()) {
2250 indices[nindices]=aodt->GetID(); nindices++;
2252 fAODMap[(Int_t)aodt->GetID()] = i;
2255 AliESDtrack *esdt = 0;
2258 esdt = (AliESDtrack*)track;
2260 esdt = new AliESDtrack(track);
2263 // single track selection
2264 okDisplaced=kFALSE; okSoftPi=kFALSE;
2265 if(fMixEvent && i<trkEntries){
2266 evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
2267 const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
2268 Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
2269 eventVtx->GetXYZ(vtxPos);
2270 vprimary->GetXYZ(primPos);
2271 eventVtx->GetCovarianceMatrix(primCov);
2272 for(Int_t ind=0;ind<3;ind++){
2273 trasl[ind]=vtxPos[ind]-primPos[ind];
2276 Bool_t isTransl=esdt->Translate(trasl,primCov);
2284 if(SingleTrkCuts(esdt,okDisplaced,okSoftPi) && nSeleTrks<trkEntries) {
2285 esdt->PropagateToDCA(fV1,fBzkG,kVeryBig);
2286 seleTrksArray.AddLast(esdt);
2287 tracksAtVertex.AddLast(new AliExternalTrackParam(*esdt));
2289 seleFlags[nSeleTrks]=0;
2290 if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2291 if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2294 if(fInputAOD) delete esdt;
2299 } // end loop on tracks
2301 // primary vertex from AOD
2303 delete fV1; fV1=NULL;
2304 vprimary->GetXYZ(pos);
2305 vprimary->GetCovarianceMatrix(cov);
2306 Double_t chi2toNDF = vprimary->GetChi2perNDF();
2307 Int_t ncontr=nindices;
2308 if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
2309 Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
2310 fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2311 fV1->SetTitle(vprimary->GetTitle());
2312 fV1->SetIndices(nindices,indices);
2314 if(indices) { delete [] indices; indices=NULL; }
2319 //-----------------------------------------------------------------------------
2320 void AliAnalysisVertexingHF::SetSelectionBitForPID(AliRDHFCuts *cuts,AliAODRecoDecayHF *rd,Int_t bit) {
2322 // Set the selection bit for PID
2324 if(cuts->GetPidHF()) {
2325 Bool_t usepid=cuts->GetIsUsePID();
2326 cuts->SetUsePID(kTRUE);
2327 if(cuts->IsSelectedPID(rd))
2328 rd->SetSelectionBit(bit);
2329 cuts->SetUsePID(usepid);
2333 //-----------------------------------------------------------------------------
2334 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2335 Bool_t &okDisplaced,Bool_t &okSoftPi) const
2337 // Check if track passes some kinematical cuts
2339 // this is needed to store the impact parameters
2340 trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2344 // Track selection, displaced tracks
2347 selectInfo = fTrackFilter->IsSelected(trk);
2349 if(selectInfo) okDisplaced=kTRUE;
2351 // Track selection, soft pions
2353 if(fDstar && fTrackFilterSoftPi) {
2354 selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2356 if(selectInfo) okSoftPi=kTRUE;
2358 if(okDisplaced || okSoftPi) return kTRUE;
2364 //-----------------------------------------------------------------------------
2365 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2367 // Transform ESDv0 to AODv0
2369 // this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2370 // and creates an AODv0 out of them
2372 Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
2373 AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2375 // create the v0 neutral track to compute the DCA to the primary vertex
2376 Double_t xyz[3], pxpypz[3];
2378 esdV0->PxPyPz(pxpypz);
2379 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2380 AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2385 Double_t d0z0[2],covd0z0[3];
2386 AliAODVertex *primVertexAOD = PrimaryVertex();
2387 trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2388 Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2389 // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2390 Double_t dcaV0DaughterToPrimVertex[2];
2391 AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2392 AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
2393 if( !posV0track || !negV0track) {
2394 if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2396 delete primVertexAOD;
2399 posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2400 // if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2402 dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
2403 negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2404 // if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2406 dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
2407 Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
2408 Double_t pmom[3],nmom[3];
2409 esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2410 esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2412 AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2413 aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2416 delete primVertexAOD;
2420 //-----------------------------------------------------------------------------
2421 void AliAnalysisVertexingHF::SetParametersAtVertex(AliESDtrack* esdt, AliExternalTrackParam* extpar) const{
2422 // Set the stored track parameters at primary vertex into AliESDtrack
2423 Double_t xyz[3], pxpypz[3], cv[21];
2424 extpar->PxPyPz(pxpypz);
2425 extpar->XvYvZv(xyz);
2426 extpar->GetCovarianceXYZPxPyPz(cv);
2427 Short_t sign=extpar->Charge();
2428 esdt->Set(xyz,pxpypz,cv,sign);