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),
98 fMassCutBeforeVertexing(kFALSE),
103 // Default constructor
105 Double_t d02[2],d03[3],d04[4];
106 fMassCalc2 = new AliAODRecoDecay(0x0,2,0,d02);
107 fMassCalc3 = new AliAODRecoDecay(0x0,3,1,d03);
108 fMassCalc4 = new AliAODRecoDecay(0x0,4,0,d04);
111 //--------------------------------------------------------------------------
112 AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
114 fInputAOD(source.fInputAOD),
115 fAODMapSize(source.fAODMapSize),
116 fAODMap(source.fAODMap),
118 fSecVtxWithKF(source.fSecVtxWithKF),
119 fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
120 fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
122 fD0toKpi(source.fD0toKpi),
123 fJPSItoEle(source.fJPSItoEle),
124 f3Prong(source.f3Prong),
125 f4Prong(source.f4Prong),
126 fDstar(source.fDstar),
127 fCascades(source.fCascades),
128 fLikeSign(source.fLikeSign),
129 fMixEvent(source.fMixEvent),
130 fTrackFilter(source.fTrackFilter),
131 fTrackFilterSoftPi(source.fTrackFilterSoftPi),
132 fCutsD0toKpi(source.fCutsD0toKpi),
133 fCutsJpsitoee(source.fCutsJpsitoee),
134 fCutsDplustoKpipi(source.fCutsDplustoKpipi),
135 fCutsDstoKKpi(source.fCutsDstoKKpi),
136 fCutsLctopKpi(source.fCutsLctopKpi),
137 fCutsLctoV0(source.fCutsLctoV0),
138 fCutsD0toKpipipi(source.fCutsD0toKpipipi),
139 fCutsDStartoKpipi(source.fCutsDStartoKpipi),
140 fListOfCuts(source.fListOfCuts),
141 fFindVertexForDstar(source.fFindVertexForDstar),
142 fFindVertexForCascades(source.fFindVertexForCascades),
143 fMassCutBeforeVertexing(source.fMassCutBeforeVertexing),
144 fMassCalc2(source.fMassCalc2),
145 fMassCalc3(source.fMassCalc3),
146 fMassCalc4(source.fMassCalc4)
152 //--------------------------------------------------------------------------
153 AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
156 // assignment operator
158 if(&source == this) return *this;
159 fInputAOD = source.fInputAOD;
160 fAODMapSize = source.fAODMapSize;
161 fBzkG = source.fBzkG;
162 fSecVtxWithKF = source.fSecVtxWithKF;
163 fRecoPrimVtxSkippingTrks = source.fRecoPrimVtxSkippingTrks;
164 fRmTrksFromPrimVtx = source.fRmTrksFromPrimVtx;
166 fD0toKpi = source.fD0toKpi;
167 fJPSItoEle = source.fJPSItoEle;
168 f3Prong = source.f3Prong;
169 f4Prong = source.f4Prong;
170 fDstar = source.fDstar;
171 fCascades = source.fCascades;
172 fLikeSign = source.fLikeSign;
173 fMixEvent = source.fMixEvent;
174 fTrackFilter = source.fTrackFilter;
175 fTrackFilterSoftPi = source.fTrackFilterSoftPi;
176 fCutsD0toKpi = source.fCutsD0toKpi;
177 fCutsJpsitoee = source.fCutsJpsitoee;
178 fCutsDplustoKpipi = source.fCutsDplustoKpipi;
179 fCutsDstoKKpi = source.fCutsDstoKKpi;
180 fCutsLctopKpi = source.fCutsLctopKpi;
181 fCutsLctoV0 = source.fCutsLctoV0;
182 fCutsD0toKpipipi = source.fCutsD0toKpipipi;
183 fCutsDStartoKpipi = source.fCutsDStartoKpipi;
184 fListOfCuts = source.fListOfCuts;
185 fFindVertexForDstar = source.fFindVertexForDstar;
186 fFindVertexForCascades = source.fFindVertexForCascades;
187 fMassCutBeforeVertexing = source.fMassCutBeforeVertexing;
188 fMassCalc2 = source.fMassCalc2;
189 fMassCalc3 = source.fMassCalc3;
190 fMassCalc4 = source.fMassCalc4;
194 //----------------------------------------------------------------------------
195 AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
197 if(fV1) { delete fV1; fV1=0; }
198 if(fTrackFilter) { delete fTrackFilter; fTrackFilter=0; }
199 if(fTrackFilterSoftPi) { delete fTrackFilterSoftPi; fTrackFilterSoftPi=0; }
200 if(fCutsD0toKpi) { delete fCutsD0toKpi; fCutsD0toKpi=0; }
201 if(fCutsJpsitoee) { delete fCutsJpsitoee; fCutsJpsitoee=0; }
202 if(fCutsDplustoKpipi) { delete fCutsDplustoKpipi; fCutsDplustoKpipi=0; }
203 if(fCutsDstoKKpi) { delete fCutsDstoKKpi; fCutsDstoKKpi=0; }
204 if(fCutsLctopKpi) { delete fCutsLctopKpi; fCutsLctopKpi=0; }
205 if(fCutsLctoV0) { delete fCutsLctoV0; fCutsLctoV0=0; }
206 if(fCutsD0toKpipipi) { delete fCutsD0toKpipipi; fCutsD0toKpipipi=0; }
207 if(fCutsDStartoKpipi) { delete fCutsDStartoKpipi; fCutsDStartoKpipi=0; }
208 if(fAODMap) { delete fAODMap; fAODMap=0; }
209 if(fMassCalc2) { delete fMassCalc2; fMassCalc2=0; }
210 if(fMassCalc3) { delete fMassCalc3; fMassCalc3=0; }
211 if(fMassCalc4) { delete fMassCalc4; fMassCalc4=0; }
213 //----------------------------------------------------------------------------
214 TList *AliAnalysisVertexingHF::FillListOfCuts() {
215 // Fill list of analysis cuts
217 TList *list = new TList();
219 list->SetName("ListOfCuts");
222 AliRDHFCutsD0toKpi *cutsD0toKpi = new AliRDHFCutsD0toKpi(*fCutsD0toKpi);
223 list->Add(cutsD0toKpi);
226 AliRDHFCutsJpsitoee *cutsJpsitoee = new AliRDHFCutsJpsitoee(*fCutsJpsitoee);
227 list->Add(cutsJpsitoee);
229 if(fCutsDplustoKpipi) {
230 AliRDHFCutsDplustoKpipi *cutsDplustoKpipi = new AliRDHFCutsDplustoKpipi(*fCutsDplustoKpipi);
231 list->Add(cutsDplustoKpipi);
234 AliRDHFCutsDstoKKpi *cutsDstoKKpi = new AliRDHFCutsDstoKKpi(*fCutsDstoKKpi);
235 list->Add(cutsDstoKKpi);
238 AliRDHFCutsLctopKpi *cutsLctopKpi = new AliRDHFCutsLctopKpi(*fCutsLctopKpi);
239 list->Add(cutsLctopKpi);
242 AliRDHFCutsLctoV0 *cutsLctoV0 = new AliRDHFCutsLctoV0(*fCutsLctoV0);
243 list->Add(cutsLctoV0);
245 if(fCutsD0toKpipipi) {
246 AliRDHFCutsD0toKpipipi *cutsD0toKpipipi = new AliRDHFCutsD0toKpipipi(*fCutsD0toKpipipi);
247 list->Add(cutsD0toKpipipi);
249 if(fCutsDStartoKpipi) {
250 AliRDHFCutsDStartoKpipi *cutsDStartoKpipi = new AliRDHFCutsDStartoKpipi(*fCutsDStartoKpipi);
251 list->Add(cutsDStartoKpipi);
254 // keep a pointer to the list
259 //----------------------------------------------------------------------------
260 void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
261 TClonesArray *aodVerticesHFTClArr,
262 TClonesArray *aodD0toKpiTClArr,
263 TClonesArray *aodJPSItoEleTClArr,
264 TClonesArray *aodCharm3ProngTClArr,
265 TClonesArray *aodCharm4ProngTClArr,
266 TClonesArray *aodDstarTClArr,
267 TClonesArray *aodCascadesTClArr,
268 TClonesArray *aodLikeSign2ProngTClArr,
269 TClonesArray *aodLikeSign3ProngTClArr)
271 // Find heavy-flavour vertex candidates
273 // Output: AOD (additional branches added)
276 TString evtype = event->IsA()->GetName();
277 fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
278 } // if we do mixing AliVEvent is a AliMixedEvent
281 AliDebug(2,"Creating HF candidates from AOD");
283 AliDebug(2,"Creating HF candidates from ESD");
286 if(!aodVerticesHFTClArr) {
287 printf("ERROR: no aodVerticesHFTClArr");
290 if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
291 printf("ERROR: no aodD0toKpiTClArr");
294 if(fJPSItoEle && !aodJPSItoEleTClArr) {
295 printf("ERROR: no aodJPSItoEleTClArr");
298 if(f3Prong && !aodCharm3ProngTClArr) {
299 printf("ERROR: no aodCharm3ProngTClArr");
302 if(f4Prong && !aodCharm4ProngTClArr) {
303 printf("ERROR: no aodCharm4ProngTClArr");
306 if(fDstar && !aodDstarTClArr) {
307 printf("ERROR: no aodDstarTClArr");
310 if(fCascades && !aodCascadesTClArr){
311 printf("ERROR: no aodCascadesTClArr ");
314 if(fLikeSign && !aodLikeSign2ProngTClArr) {
315 printf("ERROR: no aodLikeSign2ProngTClArr");
318 if(fLikeSign && f3Prong && !aodLikeSign3ProngTClArr) {
319 printf("ERROR: no aodLikeSign2ProngTClArr");
323 // delete candidates from previous event and create references
324 Int_t iVerticesHF=0,iD0toKpi=0,iJPSItoEle=0,i3Prong=0,i4Prong=0,iDstar=0,iCascades=0,iLikeSign2Prong=0,iLikeSign3Prong=0;
325 aodVerticesHFTClArr->Delete();
326 iVerticesHF = aodVerticesHFTClArr->GetEntriesFast();
327 TClonesArray &verticesHFRef = *aodVerticesHFTClArr;
328 if(fD0toKpi || fDstar) {
329 aodD0toKpiTClArr->Delete();
330 iD0toKpi = aodD0toKpiTClArr->GetEntriesFast();
333 aodJPSItoEleTClArr->Delete();
334 iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
337 aodCharm3ProngTClArr->Delete();
338 i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
341 aodCharm4ProngTClArr->Delete();
342 i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
345 aodDstarTClArr->Delete();
346 iDstar = aodDstarTClArr->GetEntriesFast();
349 aodCascadesTClArr->Delete();
350 iCascades = aodCascadesTClArr->GetEntriesFast();
353 aodLikeSign2ProngTClArr->Delete();
354 iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast();
356 if(fLikeSign && f3Prong) {
357 aodLikeSign3ProngTClArr->Delete();
358 iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast();
361 TClonesArray &aodD0toKpiRef = *aodD0toKpiTClArr;
362 TClonesArray &aodJPSItoEleRef = *aodJPSItoEleTClArr;
363 TClonesArray &aodCharm3ProngRef = *aodCharm3ProngTClArr;
364 TClonesArray &aodCharm4ProngRef = *aodCharm4ProngTClArr;
365 TClonesArray &aodDstarRef = *aodDstarTClArr;
366 TClonesArray &aodCascadesRef = *aodCascadesTClArr;
367 TClonesArray &aodLikeSign2ProngRef = *aodLikeSign2ProngTClArr;
368 TClonesArray &aodLikeSign3ProngRef = *aodLikeSign3ProngTClArr;
371 AliAODRecoDecayHF2Prong *io2Prong = 0;
372 AliAODRecoDecayHF3Prong *io3Prong = 0;
373 AliAODRecoDecayHF4Prong *io4Prong = 0;
374 AliAODRecoCascadeHF *ioCascade = 0;
376 Int_t iTrkP1,iTrkP2,iTrkN1,iTrkN2,iTrkSoftPi,trkEntries,iv0,nv0;
377 Double_t xdummy,ydummy,dcap1n1,dcap1n2,dcap2n1,dcap1p2,dcan1n2,dcap2n2,dcaV0,dcaCasc;
378 Bool_t okD0=kFALSE,okJPSI=kFALSE,ok3Prong=kFALSE,ok4Prong=kFALSE;
379 Bool_t okDstar=kFALSE,okD0fromDstar=kFALSE;
380 Bool_t okCascades=kFALSE;
381 AliESDtrack *postrack1 = 0;
382 AliESDtrack *postrack2 = 0;
383 AliESDtrack *negtrack1 = 0;
384 AliESDtrack *negtrack2 = 0;
385 AliESDtrack *trackPi = 0;
386 // AliESDtrack *posV0track = 0;
387 // AliESDtrack *negV0track = 0;
388 Float_t dcaMax = fCutsD0toKpi->GetDCACut();
389 if(fCutsJpsitoee) dcaMax=TMath::Max(dcaMax,fCutsJpsitoee->GetDCACut());
390 if(fCutsDplustoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDplustoKpipi->GetDCACut());
391 if(fCutsDstoKKpi) dcaMax=TMath::Max(dcaMax,fCutsDstoKKpi->GetDCACut());
392 if(fCutsLctopKpi) dcaMax=TMath::Max(dcaMax,fCutsLctopKpi->GetDCACut());
393 if(fCutsD0toKpipipi) dcaMax=TMath::Max(dcaMax,fCutsD0toKpipipi->GetDCACut());
394 if(fCutsDStartoKpipi) dcaMax=TMath::Max(dcaMax,fCutsDStartoKpipi->GetDCACut());
396 AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
400 fBzkG = (Double_t)event->GetMagneticField();
402 trkEntries = (Int_t)event->GetNumberOfTracks();
403 AliDebug(1,Form(" Number of tracks: %d",trkEntries));
405 nv0 = (Int_t)event->GetNumberOfV0s();
406 AliDebug(1,Form(" Number of V0s: %d",nv0));
408 if( trkEntries<2 && (trkEntries<1 || nv0<1) ) {
409 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
414 if(!fCutsD0toKpi->IsEventSelected(event)) return;
416 // call function that applies sigle-track selection,
417 // for displaced tracks and soft pions (both charges) for D*,
418 // and retrieves primary vertex
419 TObjArray seleTrksArray(trkEntries);
420 UChar_t *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi
422 Int_t *evtNumber = new Int_t[trkEntries];
423 SelectTracksAndCopyVertex(event,trkEntries,seleTrksArray,nSeleTrks,seleFlags,evtNumber);
425 //AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
426 printf(" Selected tracks: %d\n",nSeleTrks);
428 TObjArray *twoTrackArray1 = new TObjArray(2);
429 TObjArray *twoTrackArray2 = new TObjArray(2);
430 TObjArray *twoTrackArrayV0 = new TObjArray(2);
431 TObjArray *twoTrackArrayCasc = new TObjArray(2);
432 TObjArray *threeTrackArray = new TObjArray(3);
433 TObjArray *fourTrackArray = new TObjArray(4);
436 Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
438 AliAODRecoDecayHF *rd = 0;
439 AliAODRecoCascadeHF *rc = 0;
443 Bool_t massCutOK=kTRUE;
445 // LOOP ON POSITIVE TRACKS
446 for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
448 if(iTrkP1%1==0) AliDebug(1,Form(" 1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));
449 //if(iTrkP1%1==0) printf(" 1st loop on pos: track number %d of %d\n",iTrkP1,nSeleTrks);
451 // get track from tracks array
452 postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
454 if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
456 // Make cascades with V0+track
460 for(iv0=0; iv0<nv0; iv0++){
462 AliDebug(1,Form(" loop on v0s for track number %d and v0 number %d",iTrkP1,iv0));
466 v0 = ((AliAODEvent*)event)->GetV0(iv0);
468 esdV0 = ((AliESDEvent*)event)->GetV0(iv0);
470 if ( (!v0 || !v0->IsA()->InheritsFrom("AliAODv0") ) &&
471 (!esdV0 || !esdV0->IsA()->InheritsFrom("AliESDv0") ) ) continue;
474 // Get the tracks that form the V0
475 // ( parameters at primary vertex )
476 // and define an AliExternalTrackParam out of them
477 AliExternalTrackParam * posV0track;
478 AliExternalTrackParam * negV0track;
481 AliAODTrack *posVV0track = (AliAODTrack*)(v0->GetDaughter(0));
482 AliAODTrack *negVV0track = (AliAODTrack*)(v0->GetDaughter(1));
483 if( !posVV0track || !negVV0track ) continue;
485 // Apply some basic V0 daughter criteria
487 // bachelor must not be a v0-track
488 if (posVV0track->GetID() == postrack1->GetID() ||
489 negVV0track->GetID() == postrack1->GetID()) continue;
490 // reject like-sign v0
491 if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
492 // avoid ghost TPC tracks
493 if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
494 !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
495 // Get AliExternalTrackParam out of the AliAODTracks
496 Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign;
497 posVV0track->PxPyPz(pxpypz); posVV0track->XvYvZv(xyz);
498 posVV0track->GetCovarianceXYZPxPyPz(cv); sign=posVV0track->Charge();
499 posV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
500 negVV0track->PxPyPz(pxpypz); negVV0track->XvYvZv(xyz);
501 negVV0track->GetCovarianceXYZPxPyPz(cv); sign=negVV0track->Charge();
502 negV0track = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
504 AliESDtrack *posVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetPindex() ));
505 AliESDtrack *negVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetNindex() ));
506 if( !posVV0track || !negVV0track ) continue;
508 // Apply some basic V0 daughter criteria
510 // bachelor must not be a v0-track
511 if (posVV0track->GetID() == postrack1->GetID() ||
512 negVV0track->GetID() == postrack1->GetID()) continue;
513 // reject like-sign v0
514 if ( posVV0track->Charge() == negVV0track->Charge() ) continue;
515 // avoid ghost TPC tracks
516 if(!(posVV0track->GetStatus() & AliESDtrack::kTPCrefit) ||
517 !(negVV0track->GetStatus() & AliESDtrack::kTPCrefit)) continue;
518 // reject kinks (only necessary on AliESDtracks)
519 if (posVV0track->GetKinkIndex(0)>0 || negVV0track->GetKinkIndex(0)>0) continue;
520 // Get AliExternalTrackParam out of the AliESDtracks
521 posV0track = new AliExternalTrackParam(*posVV0track);
522 negV0track = new AliExternalTrackParam(*negVV0track);
524 // Define the AODv0 from ESDv0 if reading ESDs
525 v0 = TransformESDv0toAODv0(esdV0,twoTrackArrayV0);
527 if( !posV0track || !negV0track ){
528 AliDebug(1,Form(" Couldn't get the V0 daughters"));
532 // fill in the v0 two-external-track-param array
533 twoTrackArrayV0->AddAt(posV0track,0);
534 twoTrackArrayV0->AddAt(negV0track,1);
537 dcaV0 = v0->DcaV0Daughters();
539 // Define the V0 (neutral) track
540 AliNeutralTrackParam *trackV0;
542 const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
543 if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
545 Double_t xyz[3], pxpypz[3];
547 esdV0->PxPyPz(pxpypz);
548 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
549 trackV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
551 // Fill in the object array to create the cascade
552 twoTrackArrayCasc->AddAt(postrack1,0);
553 twoTrackArrayCasc->AddAt(trackV0,1);
555 // Compute the cascade vertex
556 AliAODVertex *vertexCasc = 0;
557 if(fFindVertexForCascades) {
558 // DCA between the two tracks
559 dcaCasc = postrack1->GetDCA(trackV0,fBzkG,xdummy,ydummy);
561 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
563 // assume Cascade decays at the primary vertex
564 Double_t pos[3],cov[6],chi2perNDF;
566 fV1->GetCovMatrix(cov);
567 chi2perNDF = fV1->GetChi2toNDF();
568 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
572 delete posV0track; posV0track=NULL;
573 delete negV0track; negV0track=NULL;
574 delete trackV0; trackV0=NULL;
575 if(!fInputAOD) {delete v0; v0=NULL;}
576 twoTrackArrayCasc->Clear();
580 // Create and store the Cascade if passed the cuts
581 ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,v0,dcaCasc,okCascades);
582 if(okCascades && ioCascade) {
583 AliDebug(1,Form("Storing a cascade object... "));
584 // add the vertex and the cascade to the AOD
585 AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
586 rc = new(aodCascadesRef[iCascades++])AliAODRecoCascadeHF(*ioCascade);
587 rc->SetSecondaryVtx(vCasc);
588 vCasc->SetParent(rc);
589 rc->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
590 if(!fInputAOD) vCasc->AddDaughter(v0); // just to fill ref #0 ??
591 AddRefs(vCasc,rc,event,twoTrackArrayCasc); // add the track (proton)
592 vCasc->AddDaughter(v0); // fill the 2prong V0
596 delete posV0track; posV0track=NULL;
597 delete negV0track; negV0track=NULL;
598 delete trackV0; trackV0=NULL;
599 twoTrackArrayV0->Clear();
600 twoTrackArrayCasc->Clear();
601 if(ioCascade) { delete ioCascade; ioCascade=NULL; }
602 if(vertexCasc) { delete vertexCasc; vertexCasc=NULL; }
603 if(!fInputAOD) {delete v0; v0=NULL;}
605 } // end loop on V0's
608 // If there is less than 2 particles continue
610 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
614 if(postrack1->Charge()<0 && !fLikeSign) continue;
616 // LOOP ON NEGATIVE TRACKS
617 for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
619 if(iTrkN1%1==0) AliDebug(1,Form(" 1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));
620 //if(iTrkN1%1==0) printf(" 1st loop on neg: track number %d of %d\n",iTrkN1,nSeleTrks);
622 if(iTrkN1==iTrkP1) continue;
624 // get track from tracks array
625 negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
627 if(negtrack1->Charge()>0 && !fLikeSign) continue;
629 if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
632 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
635 if(postrack1->Charge()==negtrack1->Charge()) { // like-sign
636 isLikeSign2Prong=kTRUE;
637 if(!fLikeSign) continue;
638 if(iTrkN1<iTrkP1) continue; // this is needed to avoid double-counting of like-sign
639 } else { // unlike-sign
640 isLikeSign2Prong=kFALSE;
641 if(postrack1->Charge()<0 || negtrack1->Charge()>0) continue; // this is needed to avoid double-counting of unlike-sign
643 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
648 // back to primary vertex
649 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
650 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
652 // DCA between the two tracks
653 dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
654 if(dcap1n1>dcaMax) { negtrack1=0; continue; }
657 twoTrackArray1->AddAt(postrack1,0);
658 twoTrackArray1->AddAt(negtrack1,1);
659 AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
661 twoTrackArray1->Clear();
667 if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) {
669 io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
671 if((fD0toKpi && okD0) || (fJPSItoEle && okJPSI) || (isLikeSign2Prong && (okD0 || okJPSI))) {
672 // add the vertex and the decay to the AOD
673 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
674 if(!isLikeSign2Prong) {
676 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
677 rd->SetSecondaryVtx(v2Prong);
678 v2Prong->SetParent(rd);
679 AddRefs(v2Prong,rd,event,twoTrackArray1);
682 rd = new(aodJPSItoEleRef[iJPSItoEle++])AliAODRecoDecayHF2Prong(*io2Prong);
683 rd->SetSecondaryVtx(v2Prong);
684 if(!okD0) v2Prong->SetParent(rd); // it cannot have two mothers ...
685 AddRefs(v2Prong,rd,event,twoTrackArray1);
687 } else { // isLikeSign2Prong
688 rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
689 rd->SetSecondaryVtx(v2Prong);
690 v2Prong->SetParent(rd);
691 AddRefs(v2Prong,rd,event,twoTrackArray1);
695 if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
696 // write references in io2Prong
698 AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
700 vertexp1n1->AddDaughter(postrack1);
701 vertexp1n1->AddDaughter(negtrack1);
703 io2Prong->SetSecondaryVtx(vertexp1n1);
704 //printf("---> %d %d %d %d %d\n",vertexp1n1->GetNDaughters(),iTrkP1,iTrkN1,postrack1->Charge(),negtrack1->Charge());
705 // create a track from the D0
706 AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(io2Prong);
708 // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
709 for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
711 if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
713 if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
716 if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] ||
717 evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] ||
718 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
721 if(iTrkSoftPi%1==0) AliDebug(1,Form(" 1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));
723 trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
724 if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
726 // get track from tracks array
727 trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
728 trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
730 twoTrackArrayCasc->AddAt(trackPi,0);
731 twoTrackArrayCasc->AddAt(trackD0,1);
733 AliAODVertex *vertexCasc = 0;
735 if(fFindVertexForDstar) {
736 // DCA between the two tracks
737 dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
739 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
741 // assume Dstar decays at the primary vertex
742 Double_t pos[3],cov[6],chi2perNDF;
744 fV1->GetCovMatrix(cov);
745 chi2perNDF = fV1->GetChi2toNDF();
746 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
750 twoTrackArrayCasc->Clear();
755 ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
757 // add the D0 to the AOD (if not already done)
759 AliAODVertex *v2Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexp1n1);
760 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
761 rd->SetSecondaryVtx(v2Prong);
762 v2Prong->SetParent(rd);
763 AddRefs(v2Prong,rd,event,twoTrackArray1);
764 okD0=kTRUE; // this is done to add it only once
766 // add the vertex and the cascade to the AOD
767 AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
768 rc = new(aodDstarRef[iDstar++])AliAODRecoCascadeHF(*ioCascade);
769 rc->SetSecondaryVtx(vCasc);
770 vCasc->SetParent(rc);
771 if(!fInputAOD) vCasc->AddDaughter(rd); // just to fill ref #0
772 AddRefs(vCasc,rc,event,twoTrackArrayCasc);
773 vCasc->AddDaughter(rd); // add the D0 (in ref #1)
775 twoTrackArrayCasc->Clear();
777 if(ioCascade) {delete ioCascade; ioCascade=NULL;}
778 delete vertexCasc; vertexCasc=NULL;
779 } // end loop on soft pi tracks
781 if(trackD0) {delete trackD0; trackD0=NULL;}
784 if(io2Prong) {delete io2Prong; io2Prong=NULL;}
787 twoTrackArray1->Clear();
788 if( (!f3Prong && !f4Prong) ||
789 (isLikeSign2Prong && !f3Prong) ) {
796 // 2nd LOOP ON POSITIVE TRACKS
797 for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
799 if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
801 if(iTrkP2%1==0) AliDebug(1,Form(" 2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));
803 // get track from tracks array
804 postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
806 if(postrack2->Charge()<0) continue;
808 if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
811 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
812 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
813 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
816 if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
817 if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
818 isLikeSign3Prong=kTRUE;
822 } else { // normal triplet (+-+)
823 isLikeSign3Prong=kFALSE;
825 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
826 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
827 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
831 // back to primary vertex
832 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
833 postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
834 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
835 //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
837 dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
838 if(dcap2n1>dcaMax) { postrack2=0; continue; }
839 dcap1p2 = postrack2->GetDCA(postrack1,fBzkG,xdummy,ydummy);
840 if(dcap1p2>dcaMax) { postrack2=0; continue; }
842 // check invariant mass cuts for D+,Ds,Lc
845 if(postrack2->Charge()>0) {
846 threeTrackArray->AddAt(postrack1,0);
847 threeTrackArray->AddAt(negtrack1,1);
848 threeTrackArray->AddAt(postrack2,2);
850 threeTrackArray->AddAt(negtrack1,0);
851 threeTrackArray->AddAt(postrack1,1);
852 threeTrackArray->AddAt(postrack2,2);
854 if(fMassCutBeforeVertexing)
855 massCutOK = SelectInvMass(threeTrackArray);
858 if(f3Prong && !massCutOK) {
859 threeTrackArray->Clear();
867 twoTrackArray2->AddAt(postrack2,0);
868 twoTrackArray2->AddAt(negtrack1,1);
869 AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
871 twoTrackArray2->Clear();
876 // 3 prong candidates
877 if(f3Prong && massCutOK) {
879 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
880 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
882 AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
883 if(!isLikeSign3Prong) {
884 rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
885 rd->SetSecondaryVtx(v3Prong);
886 v3Prong->SetParent(rd);
887 AddRefs(v3Prong,rd,event,threeTrackArray);
888 } else { // isLikeSign3Prong
889 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
890 rd->SetSecondaryVtx(v3Prong);
891 v3Prong->SetParent(rd);
892 AddRefs(v3Prong,rd,event,threeTrackArray);
895 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
896 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
899 // 4 prong candidates
901 // don't make 4 prong with like-sign pairs and triplets
902 && !isLikeSign2Prong && !isLikeSign3Prong
903 // track-to-track dca cuts already now
904 && dcap1n1 < fCutsD0toKpipipi->GetDCACut()
905 && dcap2n1 < fCutsD0toKpipipi->GetDCACut()) {
907 // back to primary vertex
908 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
909 postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
910 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
911 // Vertexing for these 3 (can be taken from above?)
912 threeTrackArray->AddAt(postrack1,0);
913 threeTrackArray->AddAt(negtrack1,1);
914 threeTrackArray->AddAt(postrack2,2);
915 AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
917 // 3rd LOOP ON NEGATIVE TRACKS (for 4 prong)
918 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
920 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
922 if(iTrkN2%1==0) AliDebug(1,Form(" 3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
924 // get track from tracks array
925 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
927 if(negtrack2->Charge()>0) continue;
929 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
931 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
932 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
933 evtNumber[iTrkP2]==evtNumber[iTrkN2] ||
934 evtNumber[iTrkP1]==evtNumber[iTrkN1] ||
935 evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
936 evtNumber[iTrkN1]==evtNumber[iTrkP2]) continue;
939 // back to primary vertex
940 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
941 postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
942 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
943 negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
944 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
945 if(dcap1n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
946 dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
947 if(dcap2n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
950 fourTrackArray->AddAt(postrack1,0);
951 fourTrackArray->AddAt(negtrack1,1);
952 fourTrackArray->AddAt(postrack2,2);
953 fourTrackArray->AddAt(negtrack2,3);
955 // check invariant mass cuts for D0
957 if(fMassCutBeforeVertexing)
958 massCutOK = SelectInvMass(fourTrackArray);
961 fourTrackArray->Clear();
967 AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
968 io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
970 AliAODVertex *v4Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert4PrAOD);
971 rd = new(aodCharm4ProngRef[i4Prong++])AliAODRecoDecayHF4Prong(*io4Prong);
972 rd->SetSecondaryVtx(v4Prong);
973 v4Prong->SetParent(rd);
974 AddRefs(v4Prong,rd,event,fourTrackArray);
977 if(io4Prong) {delete io4Prong; io4Prong=NULL;}
978 if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;}
979 fourTrackArray->Clear();
982 } // end loop on negative tracks
984 threeTrackArray->Clear();
992 } // end 2nd loop on positive tracks
994 twoTrackArray2->Clear();
996 // 2nd LOOP ON NEGATIVE TRACKS (for 3 prong -+-)
997 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
999 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1001 if(iTrkN2%1==0) AliDebug(1,Form(" 2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
1003 // get track from tracks array
1004 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1006 if(negtrack2->Charge()>0) continue;
1008 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1011 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
1012 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
1013 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1016 if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet
1017 if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
1018 isLikeSign3Prong=kTRUE;
1022 } else { // normal triplet (-+-)
1023 isLikeSign3Prong=kFALSE;
1026 // back to primary vertex
1027 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1028 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
1029 negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
1030 //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
1032 dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1033 if(dcap1n2>dcaMax) { negtrack2=0; continue; }
1034 dcan1n2 = negtrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
1035 if(dcan1n2>dcaMax) { negtrack2=0; continue; }
1037 threeTrackArray->AddAt(negtrack1,0);
1038 threeTrackArray->AddAt(postrack1,1);
1039 threeTrackArray->AddAt(negtrack2,2);
1041 // check invariant mass cuts for D+,Ds,Lc
1043 if(fMassCutBeforeVertexing && f3Prong)
1044 massCutOK = SelectInvMass(threeTrackArray);
1047 threeTrackArray->Clear();
1053 twoTrackArray2->AddAt(postrack1,0);
1054 twoTrackArray2->AddAt(negtrack2,1);
1056 AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1058 twoTrackArray2->Clear();
1064 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1065 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
1067 AliAODVertex *v3Prong = new(verticesHFRef[iVerticesHF++])AliAODVertex(*secVert3PrAOD);
1068 if(!isLikeSign3Prong) {
1069 rd = new(aodCharm3ProngRef[i3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1070 rd->SetSecondaryVtx(v3Prong);
1071 v3Prong->SetParent(rd);
1072 AddRefs(v3Prong,rd,event,threeTrackArray);
1073 } else { // isLikeSign3Prong
1074 rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
1075 rd->SetSecondaryVtx(v3Prong);
1076 v3Prong->SetParent(rd);
1077 AddRefs(v3Prong,rd,event,threeTrackArray);
1080 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1081 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1083 threeTrackArray->Clear();
1087 } // end 2nd loop on negative tracks
1089 twoTrackArray2->Clear();
1093 } // end 1st loop on negative tracks
1096 } // end 1st loop on positive tracks
1099 AliDebug(1,Form(" Total HF vertices in event = %d;",
1100 (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
1102 AliDebug(1,Form(" D0->Kpi in event = %d;",
1103 (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
1106 AliDebug(1,Form(" JPSI->ee in event = %d;",
1107 (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
1110 AliDebug(1,Form(" Charm->3Prong in event = %d;",
1111 (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
1114 AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
1115 (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
1118 AliDebug(1,Form(" D*->D0pi in event = %d;\n",
1119 (Int_t)aodDstarTClArr->GetEntriesFast()));
1122 AliDebug(1,Form(" cascades -> v0 + track in event = %d;\n",
1123 (Int_t)aodCascadesTClArr->GetEntriesFast()));
1126 AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
1127 (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
1129 if(fLikeSign && f3Prong) {
1130 AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
1131 (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
1135 twoTrackArray1->Delete(); delete twoTrackArray1;
1136 twoTrackArray2->Delete(); delete twoTrackArray2;
1137 twoTrackArrayCasc->Delete(); delete twoTrackArrayCasc;
1138 twoTrackArrayV0->Delete(); delete twoTrackArrayV0;
1139 threeTrackArray->Clear();
1140 threeTrackArray->Delete(); delete threeTrackArray;
1141 fourTrackArray->Delete(); delete fourTrackArray;
1142 delete [] seleFlags; seleFlags=NULL;
1143 if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
1146 seleTrksArray.Delete();
1147 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
1152 //----------------------------------------------------------------------------
1153 void AliAnalysisVertexingHF::AddRefs(AliAODVertex *v,AliAODRecoDecayHF *rd,
1154 const AliVEvent *event,
1155 const TObjArray *trkArray) const
1157 // Add the AOD tracks as daughters of the vertex (TRef)
1158 // Also add the references to the primary vertex and to the cuts
1161 AddDaughterRefs(v,event,trkArray);
1162 rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1166 rd->SetListOfCutsRef((TList*)fListOfCuts);
1167 //fListOfCuts->Print();
1168 cout<<fListOfCuts<<endl;
1169 TList *l=(TList*)rd->GetListOfCuts();
1171 if(l) {l->Print(); }else{printf("error\n");}
1176 //----------------------------------------------------------------------------
1177 void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,
1178 const AliVEvent *event,
1179 const TObjArray *trkArray) const
1181 // Add the AOD tracks as daughters of the vertex (TRef)
1183 Int_t nDg = v->GetNDaughters();
1185 if(nDg) dg = v->GetDaughter(0);
1187 if(dg) return; // daughters already added
1189 Int_t nTrks = trkArray->GetEntriesFast();
1191 AliExternalTrackParam *track = 0;
1192 AliAODTrack *aodTrack = 0;
1195 for(Int_t i=0; i<nTrks; i++) {
1196 track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1197 id = (Int_t)track->GetID();
1198 //printf("---> %d\n",id);
1199 if(id<0) continue; // this track is a AliAODRecoDecay
1200 aodTrack = (AliAODTrack*)event->GetTrack(fAODMap[id]);
1201 v->AddDaughter(aodTrack);
1206 //----------------------------------------------------------------------------
1207 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1208 TObjArray *twoTrackArray,AliVEvent *event,
1209 AliAODVertex *secVert,
1210 AliAODRecoDecayHF2Prong *rd2Prong,
1212 Bool_t &okDstar) const
1214 // Make the cascade as a 2Prong decay and check if it passes Dstar
1215 // reconstruction cuts
1219 Bool_t dummy1,dummy2,dummy3;
1221 // We use Make2Prong to construct the AliAODRecoCascadeHF
1222 // (which inherits from AliAODRecoDecayHF2Prong)
1223 AliAODRecoCascadeHF *theCascade =
1224 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1225 dummy1,dummy2,dummy3);
1226 if(!theCascade) return 0x0;
1229 AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1230 theCascade->SetCharge(trackPi->Charge());
1232 //--- selection cuts
1234 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1235 tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
1236 tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
1237 AliAODVertex *primVertexAOD=0;
1238 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1239 // take event primary vertex
1240 primVertexAOD = PrimaryVertex();
1241 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1242 rd2Prong->SetOwnPrimaryVtx(primVertexAOD);
1246 okDstar = (Bool_t)fCutsDStartoKpipi->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1248 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1249 tmpCascade->UnsetOwnPrimaryVtx();
1250 delete tmpCascade; tmpCascade=NULL;
1251 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1252 rd2Prong->UnsetOwnPrimaryVtx();
1254 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1261 //----------------------------------------------------------------------------
1262 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1263 TObjArray *twoTrackArray,AliVEvent *event,
1264 AliAODVertex *secVert,
1267 Bool_t &okCascades) const
1270 // Make the cascade as a 2Prong decay and check if it passes
1271 // cascades reconstruction cuts
1273 // AliDebug(2,Form(" building the cascade"));
1275 Bool_t dummy1,dummy2,dummy3;
1277 // We use Make2Prong to construct the AliAODRecoCascadeHF
1278 // (which inherits from AliAODRecoDecayHF2Prong)
1279 AliAODRecoCascadeHF *theCascade =
1280 (AliAODRecoCascadeHF*)Make2Prong(twoTrackArray,event,secVert,dca,
1281 dummy1,dummy2,dummy3);
1282 if(!theCascade) return 0x0;
1284 // bachelor track and charge
1285 AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1286 theCascade->SetCharge(trackBachelor->Charge());
1288 //--- selection cuts
1290 AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
1291 tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
1292 tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
1293 AliAODVertex *primVertexAOD=0;
1294 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1295 // take event primary vertex
1296 primVertexAOD = PrimaryVertex();
1297 if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex();
1298 tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
1302 Bool_t okLcksp=0, okLcLpi=0;
1303 if(fCascades && fInputAOD){
1304 okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1305 if(okCascades==1) okLcksp=1;
1306 if(okCascades==2) okLcLpi=1;
1307 if(okCascades==3) { okLcksp=1; okLcLpi=1;}
1309 else { AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied")); okCascades=kTRUE; }// no cuts implemented from ESDs
1310 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1311 tmpCascade->UnsetOwnPrimaryVtx();
1312 delete tmpCascade; tmpCascade=NULL;
1313 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1319 //-----------------------------------------------------------------------------
1320 AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
1321 TObjArray *twoTrackArray,AliVEvent *event,
1322 AliAODVertex *secVert,Double_t dca,
1323 Bool_t &okD0,Bool_t &okJPSI,
1324 Bool_t &okD0fromDstar) const
1326 // Make 2Prong candidates and check if they pass D0toKpi or BtoJPSI
1327 // reconstruction cuts
1328 // G.E.Bruno (J/psi), A.Dainese (D0->Kpi)
1330 okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
1332 Double_t px[2],py[2],pz[2],d0[2],d0err[2];
1333 AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1334 AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
1336 // propagate tracks to secondary vertex, to compute inv. mass
1337 postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1338 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1340 Double_t momentum[3];
1341 postrack->GetPxPyPz(momentum);
1342 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1343 negtrack->GetPxPyPz(momentum);
1344 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1347 // invariant mass cut (try to improve coding here..)
1348 Bool_t okMassCut=kFALSE;
1349 if(!okMassCut && fD0toKpi) if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
1350 if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
1351 if(!okMassCut && fDstar) if(SelectInvMass(3,2,px,py,pz)) okMassCut=kTRUE;
1353 AliDebug(2," candidate didn't pass mass cut");
1356 // primary vertex to be used by this candidate
1357 AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
1358 if(!primVertexAOD) return 0x0;
1361 Double_t d0z0[2],covd0z0[3];
1362 postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1364 d0err[0] = TMath::Sqrt(covd0z0[0]);
1365 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1367 d0err[1] = TMath::Sqrt(covd0z0[0]);
1369 // create the object AliAODRecoDecayHF2Prong
1370 AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
1371 the2Prong->SetOwnPrimaryVtx(primVertexAOD);
1372 UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
1373 the2Prong->SetProngIDs(2,id);
1374 delete primVertexAOD; primVertexAOD=NULL;
1377 if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
1379 if(fD0toKpi) okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1380 //if(fDebug && fD0toKpi) printf(" %d\n",(Int_t)okD0);
1381 // select J/psi from B
1382 if(fJPSItoEle) okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
1383 //if(fDebug && fJPSItoEle) printf(" %d\n",(Int_t)okJPSI);
1384 // select D0->Kpi from Dstar
1385 if(fDstar) okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
1386 //if(fDebug && fDstar) printf(" %d\n",(Int_t)okD0fromDstar);
1390 // remove the primary vertex (was used only for selection)
1391 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1392 the2Prong->UnsetOwnPrimaryVtx();
1395 // get PID info from ESD
1396 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1397 if(postrack->GetStatus()&AliESDtrack::kESDpid) postrack->GetESDpid(esdpid0);
1398 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1399 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1400 Double_t esdpid[10];
1401 for(Int_t i=0;i<5;i++) {
1402 esdpid[i] = esdpid0[i];
1403 esdpid[5+i] = esdpid1[i];
1405 the2Prong->SetPID(2,esdpid);
1409 //----------------------------------------------------------------------------
1410 AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
1411 TObjArray *threeTrackArray,AliVEvent *event,
1412 AliAODVertex *secVert,Double_t dispersion,
1413 const AliAODVertex *vertexp1n1,const AliAODVertex *vertexp2n1,
1414 Double_t dcap1n1,Double_t dcap2n1,Double_t dcap1p2,
1415 Bool_t &ok3Prong) const
1417 // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1418 // reconstruction cuts
1423 if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
1425 Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1426 Double_t momentum[3];
1429 AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1430 AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1431 AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1433 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1434 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1435 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1436 postrack1->GetPxPyPz(momentum);
1437 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1438 negtrack->GetPxPyPz(momentum);
1439 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1440 postrack2->GetPxPyPz(momentum);
1441 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1443 // invariant mass cut for D+, Ds, Lc
1444 Bool_t okMassCut=kFALSE;
1445 if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1446 if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
1448 AliDebug(2," candidate didn't pass mass cut");
1452 // primary vertex to be used by this candidate
1453 AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
1454 if(!primVertexAOD) return 0x0;
1456 Double_t d0z0[2],covd0z0[3];
1457 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1459 d0err[0] = TMath::Sqrt(covd0z0[0]);
1460 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1462 d0err[1] = TMath::Sqrt(covd0z0[0]);
1463 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1465 d0err[2] = TMath::Sqrt(covd0z0[0]);
1468 // create the object AliAODRecoDecayHF3Prong
1469 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1470 Double_t dca[3]={dcap1n1,dcap2n1,dcap1p2};
1471 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]));
1472 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]));
1473 Short_t charge=(Short_t)(postrack1->Charge()+postrack2->Charge()+negtrack->Charge());
1474 AliAODRecoDecayHF3Prong *the3Prong = new AliAODRecoDecayHF3Prong(secVert,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23,charge);
1475 the3Prong->SetOwnPrimaryVtx(primVertexAOD);
1476 UShort_t id[3]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack->GetID(),(UShort_t)postrack2->GetID()};
1477 the3Prong->SetProngIDs(3,id);
1479 delete primVertexAOD; primVertexAOD=NULL;
1481 // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1485 if(fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1486 if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1487 if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
1490 //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1492 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1493 the3Prong->UnsetOwnPrimaryVtx();
1496 // get PID info from ESD
1497 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1498 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1499 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1500 if(negtrack->GetStatus()&AliESDtrack::kESDpid) negtrack->GetESDpid(esdpid1);
1501 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1502 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1504 Double_t esdpid[15];
1505 for(Int_t i=0;i<5;i++) {
1506 esdpid[i] = esdpid0[i];
1507 esdpid[5+i] = esdpid1[i];
1508 esdpid[10+i] = esdpid2[i];
1510 the3Prong->SetPID(3,esdpid);
1514 //----------------------------------------------------------------------------
1515 AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
1516 TObjArray *fourTrackArray,AliVEvent *event,
1517 AliAODVertex *secVert,
1518 const AliAODVertex *vertexp1n1,
1519 const AliAODVertex *vertexp1n1p2,
1520 Double_t dcap1n1,Double_t dcap1n2,
1521 Double_t dcap2n1,Double_t dcap2n2,
1522 Bool_t &ok4Prong) const
1524 // Make 4Prong candidates and check if they pass D0toKpipipi
1525 // reconstruction cuts
1526 // G.E.Bruno, R.Romita
1529 if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
1531 Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1533 AliESDtrack *postrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(0);
1534 AliESDtrack *negtrack1 = (AliESDtrack*)fourTrackArray->UncheckedAt(1);
1535 AliESDtrack *postrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(2);
1536 AliESDtrack *negtrack2 = (AliESDtrack*)fourTrackArray->UncheckedAt(3);
1538 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1539 negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1540 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1541 negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1543 Double_t momentum[3];
1544 postrack1->GetPxPyPz(momentum);
1545 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
1546 negtrack1->GetPxPyPz(momentum);
1547 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
1548 postrack2->GetPxPyPz(momentum);
1549 px[2] = momentum[0]; py[2] = momentum[1]; pz[2] = momentum[2];
1550 negtrack2->GetPxPyPz(momentum);
1551 px[3] = momentum[0]; py[3] = momentum[1]; pz[3] = momentum[2];
1553 // invariant mass cut for rho or D0 (try to improve coding here..)
1554 Bool_t okMassCut=kFALSE;
1555 if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed
1556 if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) { //no PID, to be implemented with PID
1557 if(SelectInvMass(4,4,px,py,pz)) okMassCut=kTRUE;
1560 //if(fDebug) printf(" candidate didn't pass mass cut\n");
1561 //printf(" candidate didn't pass mass cut\n");
1565 // primary vertex to be used by this candidate
1566 AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
1567 if(!primVertexAOD) return 0x0;
1569 Double_t d0z0[2],covd0z0[3];
1570 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1572 d0err[0] = TMath::Sqrt(covd0z0[0]);
1573 negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1575 d0err[1] = TMath::Sqrt(covd0z0[0]);
1576 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1578 d0err[2] = TMath::Sqrt(covd0z0[0]);
1579 negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1581 d0err[3] = TMath::Sqrt(covd0z0[0]);
1584 // create the object AliAODRecoDecayHF4Prong
1585 Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1586 Double_t dca[6]={dcap1n1,0.,dcap1n2,dcap2n1,0.,dcap2n2};
1587 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]));
1588 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]));
1589 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]));
1591 AliAODRecoDecayHF4Prong *the4Prong = new AliAODRecoDecayHF4Prong(secVert,px,py,pz,d0,d0err,dca,dist12,dist3,dist4,charge);
1592 the4Prong->SetOwnPrimaryVtx(primVertexAOD);
1593 UShort_t id[4]={(UShort_t)postrack1->GetID(),(UShort_t)negtrack1->GetID(),(UShort_t)postrack2->GetID(),(UShort_t)negtrack2->GetID()};
1594 the4Prong->SetProngIDs(4,id);
1596 delete primVertexAOD; primVertexAOD=NULL;
1598 ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
1601 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1602 the4Prong->UnsetOwnPrimaryVtx();
1606 // get PID info from ESD
1607 Double_t esdpid0[5]={0.,0.,0.,0.,0.};
1608 if(postrack1->GetStatus()&AliESDtrack::kESDpid) postrack1->GetESDpid(esdpid0);
1609 Double_t esdpid1[5]={0.,0.,0.,0.,0.};
1610 if(negtrack1->GetStatus()&AliESDtrack::kESDpid) negtrack1->GetESDpid(esdpid1);
1611 Double_t esdpid2[5]={0.,0.,0.,0.,0.};
1612 if(postrack2->GetStatus()&AliESDtrack::kESDpid) postrack2->GetESDpid(esdpid2);
1613 Double_t esdpid3[5]={0.,0.,0.,0.,0.};
1614 if(negtrack2->GetStatus()&AliESDtrack::kESDpid) negtrack2->GetESDpid(esdpid3);
1616 Double_t esdpid[20];
1617 for(Int_t i=0;i<5;i++) {
1618 esdpid[i] = esdpid0[i];
1619 esdpid[5+i] = esdpid1[i];
1620 esdpid[10+i] = esdpid2[i];
1621 esdpid[15+i] = esdpid3[i];
1623 the4Prong->SetPID(4,esdpid);
1627 //-----------------------------------------------------------------------------
1628 AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
1629 AliVEvent *event) const
1631 // Returns primary vertex to be used for this candidate
1633 AliESDVertex *vertexESD = 0;
1634 AliAODVertex *vertexAOD = 0;
1637 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1638 // primary vertex from the input event
1640 vertexESD = new AliESDVertex(*fV1);
1643 // primary vertex specific to this candidate
1645 Int_t nTrks = trkArray->GetEntriesFast();
1646 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1648 if(fRecoPrimVtxSkippingTrks) {
1649 // recalculating the vertex
1651 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1652 Float_t diamondcovxy[3];
1653 event->GetDiamondCovXY(diamondcovxy);
1654 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1655 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1656 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1657 vertexer->SetVtxStart(diamond);
1658 delete diamond; diamond=NULL;
1659 if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1660 vertexer->SetOnlyFitter();
1662 Int_t skipped[1000];
1663 Int_t nTrksToSkip=0,id;
1664 AliExternalTrackParam *t = 0;
1665 for(Int_t i=0; i<nTrks; i++) {
1666 t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1667 id = (Int_t)t->GetID();
1669 skipped[nTrksToSkip++] = id;
1672 // For AOD, skip also tracks without covariance matrix
1674 Double_t covtest[21];
1675 for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1676 AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1677 if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1678 id = (Int_t)vtrack->GetID();
1680 skipped[nTrksToSkip++] = id;
1685 vertexer->SetSkipTracks(nTrksToSkip,skipped);
1686 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1688 } else if(fRmTrksFromPrimVtx) {
1689 // removing the prongs tracks
1691 TObjArray rmArray(nTrks);
1692 UShort_t *rmId = new UShort_t[nTrks];
1693 AliESDtrack *esdTrack = 0;
1695 for(Int_t i=0; i<nTrks; i++) {
1696 t = (AliESDtrack*)trkArray->UncheckedAt(i);
1697 esdTrack = new AliESDtrack(*t);
1698 rmArray.AddLast(esdTrack);
1699 if(esdTrack->GetID()>=0) {
1700 rmId[i]=(UShort_t)esdTrack->GetID();
1705 Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
1706 vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1707 delete [] rmId; rmId=NULL;
1712 if(!vertexESD) return vertexAOD;
1713 if(vertexESD->GetNContributors()<=0) {
1714 AliDebug(2,"vertexing failed");
1715 delete vertexESD; vertexESD=NULL;
1719 delete vertexer; vertexer=NULL;
1723 // convert to AliAODVertex
1724 Double_t pos[3],cov[6],chi2perNDF;
1725 vertexESD->GetXYZ(pos); // position
1726 vertexESD->GetCovMatrix(cov); //covariance matrix
1727 chi2perNDF = vertexESD->GetChi2toNDF();
1728 delete vertexESD; vertexESD=NULL;
1730 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1734 //-----------------------------------------------------------------------------
1735 void AliAnalysisVertexingHF::PrintStatus() const {
1736 // Print parameters being used
1738 //printf("Preselections:\n");
1739 // fTrackFilter->Dump();
1741 printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1743 printf("Secondary vertex with AliVertexerTracks\n");
1745 if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1746 if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1748 printf("Reconstruct D0->Kpi candidates with cuts:\n");
1749 if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
1752 printf("Reconstruct D*->D0pi candidates with cuts:\n");
1753 if(fFindVertexForDstar) {
1754 printf(" Reconstruct a secondary vertex for the D*\n");
1756 printf(" Assume the D* comes from the primary vertex\n");
1758 if(fCutsDStartoKpipi) fCutsDStartoKpipi->PrintAll();
1761 printf("Reconstruct J/psi from B candidates with cuts:\n");
1762 if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
1765 printf("Reconstruct 3 prong candidates.\n");
1766 printf(" D+->Kpipi cuts:\n");
1767 if(fCutsDplustoKpipi) fCutsDplustoKpipi->PrintAll();
1768 printf(" Ds->KKpi cuts:\n");
1769 if(fCutsDstoKKpi) fCutsDstoKKpi->PrintAll();
1770 printf(" Lc->pKpi cuts:\n");
1771 if(fCutsLctopKpi) fCutsLctopKpi->PrintAll();
1774 printf("Reconstruct 4 prong candidates.\n");
1775 printf(" D0->Kpipipi cuts:\n");
1776 if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
1779 printf("Reconstruct cascades candidates formed with v0s.\n");
1780 printf(" Lc -> k0s P & Lc -> L Pi cuts:\n");
1781 if(fCutsLctoV0) fCutsLctoV0->PrintAll();
1786 //-----------------------------------------------------------------------------
1787 AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
1788 Double_t &dispersion,Bool_t useTRefArray) const
1790 // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1792 AliESDVertex *vertexESD = 0;
1793 AliAODVertex *vertexAOD = 0;
1795 if(!fSecVtxWithKF) { // AliVertexerTracks
1797 AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1798 vertexer->SetVtxStart(fV1);
1799 vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1800 delete vertexer; vertexer=NULL;
1802 if(!vertexESD) return vertexAOD;
1804 if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
1805 AliDebug(2,"vertexing failed");
1806 delete vertexESD; vertexESD=NULL;
1810 } else { // Kalman Filter vertexer (AliKFParticle)
1812 AliKFParticle::SetField(fBzkG);
1814 AliKFVertex vertexKF;
1816 Int_t nTrks = trkArray->GetEntriesFast();
1817 for(Int_t i=0; i<nTrks; i++) {
1818 AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
1819 AliKFParticle daughterKF(*esdTrack,211);
1820 vertexKF.AddDaughter(daughterKF);
1822 vertexESD = new AliESDVertex(vertexKF.Parameters(),
1823 vertexKF.CovarianceMatrix(),
1825 vertexKF.GetNContributors());
1829 // convert to AliAODVertex
1830 Double_t pos[3],cov[6],chi2perNDF;
1831 vertexESD->GetXYZ(pos); // position
1832 vertexESD->GetCovMatrix(cov); //covariance matrix
1833 chi2perNDF = vertexESD->GetChi2toNDF();
1834 dispersion = vertexESD->GetDispersion();
1835 delete vertexESD; vertexESD=NULL;
1837 Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
1838 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
1842 //-----------------------------------------------------------------------------
1843 Bool_t AliAnalysisVertexingHF::SelectInvMass(TObjArray *trkArray) const {
1844 // Invariant mass cut on tracks
1846 Int_t nProngs=trkArray->GetEntriesFast();
1847 Int_t retval=kFALSE;
1849 Double_t momentum[3];
1850 Double_t px3[3],py3[3],pz3[3];
1851 Double_t px4[4],py4[4],pz4[4];
1852 AliESDtrack *postrack1=0;
1853 AliESDtrack *postrack2=0;
1854 AliESDtrack *negtrack1=0;
1855 AliESDtrack *negtrack2=0;
1859 // invariant mass cut for D+, Ds, Lc
1860 postrack1 = (AliESDtrack*)trkArray->UncheckedAt(0);
1861 negtrack1 = (AliESDtrack*)trkArray->UncheckedAt(1);
1862 postrack2 = (AliESDtrack*)trkArray->UncheckedAt(2);
1863 postrack1->GetPxPyPz(momentum);
1864 px3[0] = momentum[0]; py3[0] = momentum[1]; pz3[0] = momentum[2];
1865 negtrack1->GetPxPyPz(momentum);
1866 px3[1] = momentum[0]; py3[1] = momentum[1]; pz3[1] = momentum[2];
1867 postrack2->GetPxPyPz(momentum);
1868 px3[2] = momentum[0]; py3[2] = momentum[1]; pz3[2] = momentum[2];
1869 retval = SelectInvMass(2,3,px3,py3,pz3);
1872 // invariant mass cut for D0->4prong
1873 postrack1 = (AliESDtrack*)trkArray->UncheckedAt(0);
1874 negtrack1 = (AliESDtrack*)trkArray->UncheckedAt(1);
1875 postrack2 = (AliESDtrack*)trkArray->UncheckedAt(2);
1876 negtrack2 = (AliESDtrack*)trkArray->UncheckedAt(3);
1877 postrack1->GetPxPyPz(momentum);
1878 px4[0] = momentum[0]; py4[0] = momentum[1]; pz4[0] = momentum[2];
1879 negtrack1->GetPxPyPz(momentum);
1880 px4[1] = momentum[0]; py4[1] = momentum[1]; pz4[1] = momentum[2];
1881 postrack2->GetPxPyPz(momentum);
1882 px4[2] = momentum[0]; py4[2] = momentum[1]; pz4[2] = momentum[2];
1883 negtrack2->GetPxPyPz(momentum);
1884 px4[3] = momentum[0]; py4[3] = momentum[1]; pz4[3] = momentum[2];
1885 retval = SelectInvMass(4,4,px4,py4,pz4);
1888 printf("SelectInvMass(): wrong decay selection\n");
1894 //-----------------------------------------------------------------------------
1895 Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1899 Double_t *pz) const {
1900 // Check invariant mass cut
1902 UInt_t pdg2[2],pdg3[3],pdg4[4];
1905 Bool_t retval=kFALSE;
1909 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
1910 pdg2[0]=211; pdg2[1]=321;
1911 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1912 minv = fMassCalc2->InvMass(nprongs,pdg2);
1913 if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
1914 pdg2[0]=321; pdg2[1]=211;
1915 minv = fMassCalc2->InvMass(nprongs,pdg2);
1916 if(TMath::Abs(minv-mPDG)<fCutsD0toKpi->GetMassCut()) retval=kTRUE;
1919 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
1920 pdg2[0]=11; pdg2[1]=11;
1921 mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
1922 minv = fMassCalc2->InvMass(nprongs,pdg2);
1923 if(TMath::Abs(minv-mPDG)<fCutsJpsitoee->GetMassCut()) retval=kTRUE;
1925 case 2: // D+->Kpipi
1926 fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
1927 pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
1928 mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
1929 minv = fMassCalc3->InvMass(nprongs,pdg3);
1930 if(TMath::Abs(minv-mPDG)<fCutsDplustoKpipi->GetMassCut()) retval=kTRUE;
1932 pdg3[0]=321; pdg3[1]=321; pdg3[2]=211;
1933 mPDG=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1934 minv = fMassCalc3->InvMass(nprongs,pdg3);
1935 if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
1936 pdg3[0]=211; pdg3[1]=321; pdg3[2]=321;
1937 minv = fMassCalc3->InvMass(nprongs,pdg3);
1938 if(TMath::Abs(minv-mPDG)<fCutsDstoKKpi->GetMassCut()) retval=kTRUE;
1940 pdg3[0]=2212; pdg3[1]=321; pdg3[2]=211;
1941 mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1942 minv = fMassCalc3->InvMass(nprongs,pdg3);
1943 if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
1944 pdg3[0]=211; pdg3[1]=321; pdg3[2]=2212;
1945 minv = fMassCalc3->InvMass(nprongs,pdg3);
1946 if(TMath::Abs(minv-mPDG)<fCutsLctopKpi->GetMassCut()) retval=kTRUE;
1949 fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
1950 pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
1951 mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
1952 minv = fMassCalc2->InvMass(nprongs,pdg2);
1953 if(TMath::Abs(minv-mPDG)<fCutsDStartoKpipi->GetMassCut()) retval=kTRUE;
1955 case 4: // D0->Kpipipi without PID
1956 fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
1957 pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
1958 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1959 minv = fMassCalc4->InvMass(nprongs,pdg4);
1960 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1961 pdg4[0]=211; pdg4[1]=321; pdg4[2]=211; pdg4[3]=211;
1962 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1963 minv = fMassCalc4->InvMass(nprongs,pdg4);
1964 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1965 pdg4[0]=211; pdg4[1]=211; pdg4[2]=321; pdg4[3]=211;
1966 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1967 minv = fMassCalc4->InvMass(nprongs,pdg4);
1968 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1969 pdg4[0]=211; pdg4[1]=211; pdg4[2]=211; pdg4[3]=321;
1970 mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1971 minv = fMassCalc4->InvMass(nprongs,pdg4);
1972 if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
1975 printf("SelectInvMass(): wrong decay selection\n");
1981 //-----------------------------------------------------------------------------
1982 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
1984 TObjArray &seleTrksArray,Int_t &nSeleTrks,
1985 UChar_t *seleFlags,Int_t *evtNumber)
1987 // Apply single-track preselection.
1988 // Fill a TObjArray with selected tracks (for displaced vertices or
1989 // soft pion from D*). Selection flag stored in seleFlags.
1990 // Create the AliESDVertex object (convert from AliAODVertex if necessary)
1991 // In case of AOD input, also fill fAODMap for track index<->ID
1993 const AliVVertex *vprimary = event->GetPrimaryVertex();
1995 if(fV1) { delete fV1; fV1=NULL; }
1996 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
1999 UShort_t *indices = 0;
2000 Double_t pos[3],cov[6];
2002 if(!fInputAOD) { // ESD
2003 fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
2005 vprimary->GetXYZ(pos);
2006 vprimary->GetCovarianceMatrix(cov);
2007 fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2008 indices = new UShort_t[event->GetNumberOfTracks()];
2009 fAODMapSize = 100000;
2010 fAODMap = new Int_t[fAODMapSize];
2014 Int_t entries = (Int_t)event->GetNumberOfTracks();
2015 Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
2018 // transfer ITS tracks from event to arrays
2019 for(Int_t i=0; i<entries; i++) {
2021 track = (AliVTrack*)event->GetTrack(i);
2023 // skip pure ITS SA tracks
2024 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2026 // TEMPORARY: check that the cov matrix is there
2027 Double_t covtest[21];
2028 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2032 AliAODTrack *aodt = (AliAODTrack*)track;
2033 if(aodt->GetUsedForPrimVtxFit()) {
2034 indices[nindices]=aodt->GetID(); nindices++;
2036 fAODMap[(Int_t)aodt->GetID()] = i;
2039 AliESDtrack *esdt = 0;
2042 esdt = (AliESDtrack*)track;
2044 esdt = new AliESDtrack(track);
2047 // single track selection
2048 okDisplaced=kFALSE; okSoftPi=kFALSE;
2049 if(fMixEvent && i<trkEntries){
2050 evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
2051 const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
2052 Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
2053 eventVtx->GetXYZ(vtxPos);
2054 vprimary->GetXYZ(primPos);
2055 eventVtx->GetCovarianceMatrix(primCov);
2056 for(Int_t ind=0;ind<3;ind++){
2057 trasl[ind]=vtxPos[ind]-primPos[ind];
2060 Bool_t isTransl=esdt->Translate(trasl,primCov);
2068 if(SingleTrkCuts(esdt,okDisplaced,okSoftPi) && nSeleTrks<trkEntries) {
2069 seleTrksArray.AddLast(esdt);
2070 seleFlags[nSeleTrks]=0;
2071 if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2072 if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2075 if(fInputAOD) delete esdt;
2080 } // end loop on tracks
2082 // primary vertex from AOD
2084 delete fV1; fV1=NULL;
2085 vprimary->GetXYZ(pos);
2086 vprimary->GetCovarianceMatrix(cov);
2087 Double_t chi2toNDF = vprimary->GetChi2perNDF();
2088 Int_t ncontr=nindices;
2089 if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
2090 Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
2091 fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2092 fV1->SetTitle(vprimary->GetTitle());
2093 fV1->SetIndices(nindices,indices);
2095 if(indices) { delete [] indices; indices=NULL; }
2100 //-----------------------------------------------------------------------------
2101 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2102 Bool_t &okDisplaced,Bool_t &okSoftPi) const
2104 // Check if track passes some kinematical cuts
2106 // this is needed to store the impact parameters
2107 trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2111 // Track selection, displaced tracks
2114 selectInfo = fTrackFilter->IsSelected(trk);
2116 if(selectInfo) okDisplaced=kTRUE;
2118 // Track selection, soft pions
2120 if(fDstar && fTrackFilterSoftPi) {
2121 selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2123 if(selectInfo) okSoftPi=kTRUE;
2125 if(okDisplaced || okSoftPi) return kTRUE;
2131 //-----------------------------------------------------------------------------
2132 AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2134 // Transform ESDv0 to AODv0
2136 // this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2137 // and creates an AODv0 out of them
2139 Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
2140 AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2142 // create the v0 neutral track to compute the DCA to the primary vertex
2143 Double_t xyz[3], pxpypz[3];
2145 esdV0->PxPyPz(pxpypz);
2146 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2147 AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2148 if(!trackesdV0) return 0;
2149 Double_t d0z0[2],covd0z0[3];
2150 AliAODVertex *primVertexAOD = PrimaryVertex();
2151 trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2152 Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2153 // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2154 Double_t dcaV0DaughterToPrimVertex[2];
2155 AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2156 AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
2157 if( !posV0track || !negV0track) {
2158 if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2161 posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2162 // if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2164 dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
2165 negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2166 // if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2168 dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
2169 Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
2170 Double_t pmom[3],nmom[3];
2171 esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2172 esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2174 AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2175 aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2177 if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2178 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2182 //-----------------------------------------------------------------------------