Added flag to separate prompt and secondary charm in MC (Giacomo, Chiara)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisVertexingHF.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//----------------------------------------------------------------------------
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.
22//
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>
28#include <TFile.h>
29#include <TDatabasePDG.h>
30#include <TString.h>
31#include <TList.h>
32#include "AliLog.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"
62#include "AliESDv0.h"
63#include "AliAODv0.h"
64
65ClassImp(AliAnalysisVertexingHF)
66
67//----------------------------------------------------------------------------
68AliAnalysisVertexingHF::AliAnalysisVertexingHF():
69fInputAOD(kFALSE),
70fAODMapSize(0),
71fAODMap(0),
72fBzkG(0.),
73fSecVtxWithKF(kFALSE),
74fRecoPrimVtxSkippingTrks(kFALSE),
75fRmTrksFromPrimVtx(kFALSE),
76fV1(0x0),
77fD0toKpi(kTRUE),
78fJPSItoEle(kTRUE),
79f3Prong(kTRUE),
80f4Prong(kTRUE),
81fDstar(kTRUE),
82fCascades(kTRUE),
83fLikeSign(kFALSE),
84fMixEvent(kFALSE),
85fTrackFilter(0x0),
86fTrackFilterSoftPi(0x0),
87fCutsD0toKpi(0x0),
88fCutsJpsitoee(0x0),
89fCutsDplustoKpipi(0x0),
90fCutsDstoKKpi(0x0),
91fCutsLctopKpi(0x0),
92fCutsLctoV0(0x0),
93fCutsD0toKpipipi(0x0),
94fCutsDStartoKpipi(0x0),
95fListOfCuts(0x0),
96fFindVertexForDstar(kTRUE),
97fFindVertexForCascades(kTRUE),
98fMassCutBeforeVertexing(kFALSE),
99fMassCalc2(0),
100fMassCalc3(0),
101fMassCalc4(0)
102{
103 // Default constructor
104
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);
109
110}
111//--------------------------------------------------------------------------
112AliAnalysisVertexingHF::AliAnalysisVertexingHF(const AliAnalysisVertexingHF &source) :
113TNamed(source),
114fInputAOD(source.fInputAOD),
115fAODMapSize(source.fAODMapSize),
116fAODMap(source.fAODMap),
117fBzkG(source.fBzkG),
118fSecVtxWithKF(source.fSecVtxWithKF),
119fRecoPrimVtxSkippingTrks(source.fRecoPrimVtxSkippingTrks),
120fRmTrksFromPrimVtx(source.fRmTrksFromPrimVtx),
121fV1(source.fV1),
122fD0toKpi(source.fD0toKpi),
123fJPSItoEle(source.fJPSItoEle),
124f3Prong(source.f3Prong),
125f4Prong(source.f4Prong),
126fDstar(source.fDstar),
127fCascades(source.fCascades),
128fLikeSign(source.fLikeSign),
129fMixEvent(source.fMixEvent),
130fTrackFilter(source.fTrackFilter),
131fTrackFilterSoftPi(source.fTrackFilterSoftPi),
132fCutsD0toKpi(source.fCutsD0toKpi),
133fCutsJpsitoee(source.fCutsJpsitoee),
134fCutsDplustoKpipi(source.fCutsDplustoKpipi),
135fCutsDstoKKpi(source.fCutsDstoKKpi),
136fCutsLctopKpi(source.fCutsLctopKpi),
137fCutsLctoV0(source.fCutsLctoV0),
138fCutsD0toKpipipi(source.fCutsD0toKpipipi),
139fCutsDStartoKpipi(source.fCutsDStartoKpipi),
140fListOfCuts(source.fListOfCuts),
141fFindVertexForDstar(source.fFindVertexForDstar),
142fFindVertexForCascades(source.fFindVertexForCascades),
143fMassCutBeforeVertexing(source.fMassCutBeforeVertexing),
144fMassCalc2(source.fMassCalc2),
145fMassCalc3(source.fMassCalc3),
146fMassCalc4(source.fMassCalc4)
147{
148 //
149 // Copy constructor
150 //
151}
152//--------------------------------------------------------------------------
153AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVertexingHF &source)
154{
155 //
156 // assignment operator
157 //
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;
165 fV1 = source.fV1;
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;
191
192 return *this;
193}
194//----------------------------------------------------------------------------
195AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
196 // Destructor
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; }
212}
213//----------------------------------------------------------------------------
214TList *AliAnalysisVertexingHF::FillListOfCuts() {
215 // Fill list of analysis cuts
216
217 TList *list = new TList();
218 list->SetOwner();
219 list->SetName("ListOfCuts");
220
221 if(fCutsD0toKpi) {
222 AliRDHFCutsD0toKpi *cutsD0toKpi = new AliRDHFCutsD0toKpi(*fCutsD0toKpi);
223 list->Add(cutsD0toKpi);
224 }
225 if(fCutsJpsitoee) {
226 AliRDHFCutsJpsitoee *cutsJpsitoee = new AliRDHFCutsJpsitoee(*fCutsJpsitoee);
227 list->Add(cutsJpsitoee);
228 }
229 if(fCutsDplustoKpipi) {
230 AliRDHFCutsDplustoKpipi *cutsDplustoKpipi = new AliRDHFCutsDplustoKpipi(*fCutsDplustoKpipi);
231 list->Add(cutsDplustoKpipi);
232 }
233 if(fCutsDstoKKpi) {
234 AliRDHFCutsDstoKKpi *cutsDstoKKpi = new AliRDHFCutsDstoKKpi(*fCutsDstoKKpi);
235 list->Add(cutsDstoKKpi);
236 }
237 if(fCutsLctopKpi) {
238 AliRDHFCutsLctopKpi *cutsLctopKpi = new AliRDHFCutsLctopKpi(*fCutsLctopKpi);
239 list->Add(cutsLctopKpi);
240 }
241 if(fCutsLctoV0){
242 AliRDHFCutsLctoV0 *cutsLctoV0 = new AliRDHFCutsLctoV0(*fCutsLctoV0);
243 list->Add(cutsLctoV0);
244 }
245 if(fCutsD0toKpipipi) {
246 AliRDHFCutsD0toKpipipi *cutsD0toKpipipi = new AliRDHFCutsD0toKpipipi(*fCutsD0toKpipipi);
247 list->Add(cutsD0toKpipipi);
248 }
249 if(fCutsDStartoKpipi) {
250 AliRDHFCutsDStartoKpipi *cutsDStartoKpipi = new AliRDHFCutsDStartoKpipi(*fCutsDStartoKpipi);
251 list->Add(cutsDStartoKpipi);
252 }
253
254 // keep a pointer to the list
255 fListOfCuts = list;
256
257 return list;
258}
259//----------------------------------------------------------------------------
260void 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)
270{
271 // Find heavy-flavour vertex candidates
272 // Input: ESD or AOD
273 // Output: AOD (additional branches added)
274
275 if(!fMixEvent){
276 TString evtype = event->IsA()->GetName();
277 fInputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
278 } // if we do mixing AliVEvent is a AliMixedEvent
279
280 if(fInputAOD) {
281 AliDebug(2,"Creating HF candidates from AOD");
282 } else {
283 AliDebug(2,"Creating HF candidates from ESD");
284 }
285
286 if(!aodVerticesHFTClArr) {
287 printf("ERROR: no aodVerticesHFTClArr");
288 return;
289 }
290 if((fD0toKpi || fDstar) && !aodD0toKpiTClArr) {
291 printf("ERROR: no aodD0toKpiTClArr");
292 return;
293 }
294 if(fJPSItoEle && !aodJPSItoEleTClArr) {
295 printf("ERROR: no aodJPSItoEleTClArr");
296 return;
297 }
298 if(f3Prong && !aodCharm3ProngTClArr) {
299 printf("ERROR: no aodCharm3ProngTClArr");
300 return;
301 }
302 if(f4Prong && !aodCharm4ProngTClArr) {
303 printf("ERROR: no aodCharm4ProngTClArr");
304 return;
305 }
306 if(fDstar && !aodDstarTClArr) {
307 printf("ERROR: no aodDstarTClArr");
308 return;
309 }
310 if(fCascades && !aodCascadesTClArr){
311 printf("ERROR: no aodCascadesTClArr ");
312 return;
313 }
314 if(fLikeSign && !aodLikeSign2ProngTClArr) {
315 printf("ERROR: no aodLikeSign2ProngTClArr");
316 return;
317 }
318 if(fLikeSign && f3Prong && !aodLikeSign3ProngTClArr) {
319 printf("ERROR: no aodLikeSign2ProngTClArr");
320 return;
321 }
322
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();
331 }
332 if(fJPSItoEle) {
333 aodJPSItoEleTClArr->Delete();
334 iJPSItoEle = aodJPSItoEleTClArr->GetEntriesFast();
335 }
336 if(f3Prong) {
337 aodCharm3ProngTClArr->Delete();
338 i3Prong = aodCharm3ProngTClArr->GetEntriesFast();
339 }
340 if(f4Prong) {
341 aodCharm4ProngTClArr->Delete();
342 i4Prong = aodCharm4ProngTClArr->GetEntriesFast();
343 }
344 if(fDstar) {
345 aodDstarTClArr->Delete();
346 iDstar = aodDstarTClArr->GetEntriesFast();
347 }
348 if(fCascades) {
349 aodCascadesTClArr->Delete();
350 iCascades = aodCascadesTClArr->GetEntriesFast();
351 }
352 if(fLikeSign) {
353 aodLikeSign2ProngTClArr->Delete();
354 iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast();
355 }
356 if(fLikeSign && f3Prong) {
357 aodLikeSign3ProngTClArr->Delete();
358 iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast();
359 }
360
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;
369
370
371 AliAODRecoDecayHF2Prong *io2Prong = 0;
372 AliAODRecoDecayHF3Prong *io3Prong = 0;
373 AliAODRecoDecayHF4Prong *io4Prong = 0;
374 AliAODRecoCascadeHF *ioCascade = 0;
375
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());
395
396 AliDebug(2,Form(" dca cut set to %f cm",dcaMax));
397
398
399 // get Bz
400 fBzkG = (Double_t)event->GetMagneticField();
401
402 trkEntries = (Int_t)event->GetNumberOfTracks();
403 AliDebug(1,Form(" Number of tracks: %d",trkEntries));
404
405 nv0 = (Int_t)event->GetNumberOfV0s();
406 AliDebug(1,Form(" Number of V0s: %d",nv0));
407
408 if( trkEntries<2 && (trkEntries<1 || nv0<1) ) {
409 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
410 return;
411 }
412
413 // event selection
414 if(!fCutsD0toKpi->IsEventSelected(event)) return;
415
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
421 Int_t nSeleTrks=0;
422 Int_t *evtNumber = new Int_t[trkEntries];
423 SelectTracksAndCopyVertex(event,seleTrksArray,nSeleTrks,seleFlags,evtNumber);
424
425 //AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
426 printf(" Selected tracks: %d\n",nSeleTrks);
427
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);
434
435 Double_t dispersion;
436 Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
437
438 AliAODRecoDecayHF *rd = 0;
439 AliAODRecoCascadeHF *rc = 0;
440 AliAODv0 *v0 = 0;
441 AliESDv0 *esdV0 = 0;
442
443 Bool_t massCutOK=kTRUE;
444
445 // LOOP ON POSITIVE TRACKS
446 for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
447
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);
450
451 // get track from tracks array
452 postrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP1);
453
454 if(!TESTBIT(seleFlags[iTrkP1],kBitDispl)) continue;
455
456 // Make cascades with V0+track
457 //
458 if(fCascades) {
459 // loop on V0's
460 for(iv0=0; iv0<nv0; iv0++){
461
462 AliDebug(1,Form(" loop on v0s for track number %d and v0 number %d",iTrkP1,iv0));
463
464 // Get the V0
465 if(fInputAOD) {
466 v0 = ((AliAODEvent*)event)->GetV0(iv0);
467 } else {
468 esdV0 = ((AliESDEvent*)event)->GetV0(iv0);
469 }
470 if ( (!v0 || !v0->IsA()->InheritsFrom("AliAODv0") ) &&
471 (!esdV0 || !esdV0->IsA()->InheritsFrom("AliESDv0") ) ) continue;
472
473
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;
479
480 if(fInputAOD){
481 AliAODTrack *posVV0track = (AliAODTrack*)(v0->GetDaughter(0));
482 AliAODTrack *negVV0track = (AliAODTrack*)(v0->GetDaughter(1));
483 if( !posVV0track || !negVV0track ) continue;
484 //
485 // Apply some basic V0 daughter criteria
486 //
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);
503 } else {
504 AliESDtrack *posVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetPindex() ));
505 AliESDtrack *negVV0track = (AliESDtrack*)(event->GetTrack( esdV0->GetNindex() ));
506 if( !posVV0track || !negVV0track ) continue;
507 //
508 // Apply some basic V0 daughter criteria
509 //
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);
523
524 // Define the AODv0 from ESDv0 if reading ESDs
525 v0 = TransformESDv0toAODv0(esdV0,twoTrackArrayV0);
526 }
527 if( !posV0track || !negV0track ){
528 AliDebug(1,Form(" Couldn't get the V0 daughters"));
529 continue;
530 }
531
532 // fill in the v0 two-external-track-param array
533 twoTrackArrayV0->AddAt(posV0track,0);
534 twoTrackArrayV0->AddAt(negV0track,1);
535
536 // Get the V0 dca
537 dcaV0 = v0->DcaV0Daughters();
538
539 // Define the V0 (neutral) track
540 AliNeutralTrackParam *trackV0;
541 if(fInputAOD) {
542 const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
543 if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
544 } else {
545 Double_t xyz[3], pxpypz[3];
546 esdV0->XvYvZv(xyz);
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);
550 }
551 // Fill in the object array to create the cascade
552 twoTrackArrayCasc->AddAt(postrack1,0);
553 twoTrackArrayCasc->AddAt(trackV0,1);
554
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);
560 // Vertexing+
561 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
562 } else {
563 // assume Cascade decays at the primary vertex
564 Double_t pos[3],cov[6],chi2perNDF;
565 fV1->GetXYZ(pos);
566 fV1->GetCovMatrix(cov);
567 chi2perNDF = fV1->GetChi2toNDF();
568 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
569 dcaCasc = 0.;
570 }
571 if(!vertexCasc) {
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();
577 continue;
578 }
579
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
593 }
594
595 // Clean up
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;}
604
605 } // end loop on V0's
606 }
607
608 // If there is less than 2 particles exit
609 if(trkEntries<2) {
610 AliDebug(1,Form(" Not enough tracks: %d",trkEntries));
611 return;
612 }
613
614 if(postrack1->Charge()<0 && !fLikeSign) continue;
615
616 // LOOP ON NEGATIVE TRACKS
617 for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
618
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);
621
622 if(iTrkN1==iTrkP1) continue;
623
624 // get track from tracks array
625 negtrack1 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN1);
626
627 if(negtrack1->Charge()>0 && !fLikeSign) continue;
628
629 if(!TESTBIT(seleFlags[iTrkN1],kBitDispl)) continue;
630
631 if(fMixEvent) {
632 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
633 }
634
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
642 if(fMixEvent) {
643 if(evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
644 }
645
646 }
647
648 // back to primary vertex
649 postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
650 negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
651
652 // DCA between the two tracks
653 dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
654 if(dcap1n1>dcaMax) { negtrack1=0; continue; }
655
656 // Vertexing
657 twoTrackArray1->AddAt(postrack1,0);
658 twoTrackArray1->AddAt(negtrack1,1);
659 AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
660 if(!vertexp1n1) {
661 twoTrackArray1->Clear();
662 negtrack1=0;
663 continue;
664 }
665
666 // 2 prong candidate
667 if(fD0toKpi || fJPSItoEle || fDstar || fLikeSign) {
668
669 io2Prong = Make2Prong(twoTrackArray1,event,vertexp1n1,dcap1n1,okD0,okJPSI,okD0fromDstar);
670
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) {
675 if(okD0) {
676 rd = new(aodD0toKpiRef[iD0toKpi++])AliAODRecoDecayHF2Prong(*io2Prong);
677 rd->SetSecondaryVtx(v2Prong);
678 v2Prong->SetParent(rd);
679 AddRefs(v2Prong,rd,event,twoTrackArray1);
680 }
681 if(okJPSI) {
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);
686 }
687 } else { // isLikeSign2Prong
688 rd = new(aodLikeSign2ProngRef[iLikeSign2Prong++])AliAODRecoDecayHF2Prong(*io2Prong);
689 rd->SetSecondaryVtx(v2Prong);
690 v2Prong->SetParent(rd);
691 AddRefs(v2Prong,rd,event,twoTrackArray1);
692 }
693 }
694 // D* candidates
695 if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
696 // write references in io2Prong
697 if(fInputAOD) {
698 AddDaughterRefs(vertexp1n1,event,twoTrackArray1);
699 } else {
700 vertexp1n1->AddDaughter(postrack1);
701 vertexp1n1->AddDaughter(negtrack1);
702 }
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);
707
708 // LOOP ON TRACKS THAT PASSED THE SOFT PION CUTS
709 for(iTrkSoftPi=0; iTrkSoftPi<nSeleTrks; iTrkSoftPi++) {
710
711 if(iTrkSoftPi==iTrkP1 || iTrkSoftPi==iTrkN1) continue;
712
713 if(!TESTBIT(seleFlags[iTrkSoftPi],kBitSoftPi)) continue;
714
715 if(fMixEvent) {
716 if(evtNumber[iTrkP1]==evtNumber[iTrkSoftPi] ||
717 evtNumber[iTrkN1]==evtNumber[iTrkSoftPi] ||
718 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
719 }
720
721 if(iTrkSoftPi%1==0) AliDebug(1,Form(" 1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));
722
723 trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
724 if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
725
726 // get track from tracks array
727 trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
728 trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
729
730 twoTrackArrayCasc->AddAt(trackPi,0);
731 twoTrackArrayCasc->AddAt(trackD0,1);
732
733 AliAODVertex *vertexCasc = 0;
734
735 if(fFindVertexForDstar) {
736 // DCA between the two tracks
737 dcaCasc = trackPi->GetDCA(trackD0,fBzkG,xdummy,ydummy);
738 // Vertexing
739 vertexCasc = ReconstructSecondaryVertex(twoTrackArrayCasc,dispersion,kFALSE);
740 } else {
741 // assume Dstar decays at the primary vertex
742 Double_t pos[3],cov[6],chi2perNDF;
743 fV1->GetXYZ(pos);
744 fV1->GetCovMatrix(cov);
745 chi2perNDF = fV1->GetChi2toNDF();
746 vertexCasc = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,2);
747 dcaCasc = 0.;
748 }
749 if(!vertexCasc) {
750 twoTrackArrayCasc->Clear();
751 trackPi=0;
752 continue;
753 }
754
755 ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,io2Prong,dcaCasc,okDstar);
756 if(okDstar) {
757 // add the D0 to the AOD (if not already done)
758 if(!okD0) {
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
765 }
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)
774 }
775 twoTrackArrayCasc->Clear();
776 trackPi=0;
777 if(ioCascade) {delete ioCascade; ioCascade=NULL;}
778 delete vertexCasc; vertexCasc=NULL;
779 } // end loop on soft pi tracks
780
781 if(trackD0) {delete trackD0; trackD0=NULL;}
782
783 }
784 if(io2Prong) {delete io2Prong; io2Prong=NULL;}
785 }
786
787 twoTrackArray1->Clear();
788 if( (!f3Prong && !f4Prong) ||
789 (isLikeSign2Prong && !f3Prong) ) {
790 negtrack1=0;
791 delete vertexp1n1;
792 continue;
793 }
794
795
796 // 2nd LOOP ON POSITIVE TRACKS
797 for(iTrkP2=iTrkP1+1; iTrkP2<nSeleTrks; iTrkP2++) {
798
799 if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
800
801 if(iTrkP2%1==0) AliDebug(1,Form(" 2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));
802
803 // get track from tracks array
804 postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
805
806 if(postrack2->Charge()<0) continue;
807
808 if(!TESTBIT(seleFlags[iTrkP2],kBitDispl)) continue;
809
810 if(fMixEvent) {
811 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
812 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
813 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
814 }
815
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;
819 } else { // not ok
820 continue;
821 }
822 } else { // normal triplet (+-+)
823 isLikeSign3Prong=kFALSE;
824 if(fMixEvent) {
825 if(evtNumber[iTrkP1]==evtNumber[iTrkP2] ||
826 evtNumber[iTrkN1]==evtNumber[iTrkP2] ||
827 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
828 }
829 }
830
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());
836
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; }
841
842 // check invariant mass cuts for D+,Ds,Lc
843 massCutOK=kTRUE;
844 if(f3Prong) {
845 if(postrack2->Charge()>0) {
846 threeTrackArray->AddAt(postrack1,0);
847 threeTrackArray->AddAt(negtrack1,1);
848 threeTrackArray->AddAt(postrack2,2);
849 } else {
850 threeTrackArray->AddAt(negtrack1,0);
851 threeTrackArray->AddAt(postrack1,1);
852 threeTrackArray->AddAt(postrack2,2);
853 }
854 if(fMassCutBeforeVertexing)
855 massCutOK = SelectInvMass(threeTrackArray);
856 }
857
858 if(f3Prong && !massCutOK) {
859 threeTrackArray->Clear();
860 if(!f4Prong) {
861 postrack2=0;
862 continue;
863 }
864 }
865
866 // Vertexing
867 twoTrackArray2->AddAt(postrack2,0);
868 twoTrackArray2->AddAt(negtrack1,1);
869 AliAODVertex *vertexp2n1 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
870 if(!vertexp2n1) {
871 twoTrackArray2->Clear();
872 postrack2=0;
873 continue;
874 }
875
876 // 3 prong candidates
877 if(f3Prong && massCutOK) {
878
879 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
880 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp2n1,dcap1n1,dcap2n1,dcap1p2,ok3Prong);
881 if(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);
893 }
894 }
895 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
896 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
897 }
898
899 // 4 prong candidates
900 if(f4Prong
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()) {
906
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);
916
917 // 3rd LOOP ON NEGATIVE TRACKS (for 4 prong)
918 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
919
920 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
921
922 if(iTrkN2%1==0) AliDebug(1,Form(" 3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
923
924 // get track from tracks array
925 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
926
927 if(negtrack2->Charge()>0) continue;
928
929 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
930 if(fMixEvent){
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;
937 }
938
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; }
948
949
950 fourTrackArray->AddAt(postrack1,0);
951 fourTrackArray->AddAt(negtrack1,1);
952 fourTrackArray->AddAt(postrack2,2);
953 fourTrackArray->AddAt(negtrack2,3);
954
955 // check invariant mass cuts for D0
956 massCutOK=kTRUE;
957 if(fMassCutBeforeVertexing)
958 massCutOK = SelectInvMass(fourTrackArray);
959
960 if(!massCutOK) {
961 fourTrackArray->Clear();
962 negtrack2=0;
963 continue;
964 }
965
966 // Vertexing
967 AliAODVertex* secVert4PrAOD = ReconstructSecondaryVertex(fourTrackArray,dispersion);
968 io4Prong = Make4Prong(fourTrackArray,event,secVert4PrAOD,vertexp1n1,vertexp1n1p2,dcap1n1,dcap1n2,dcap2n1,dcap2n2,ok4Prong);
969 if(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);
975 }
976
977 if(io4Prong) {delete io4Prong; io4Prong=NULL;}
978 if(secVert4PrAOD) {delete secVert4PrAOD; secVert4PrAOD=NULL;}
979 fourTrackArray->Clear();
980 negtrack2 = 0;
981
982 } // end loop on negative tracks
983
984 threeTrackArray->Clear();
985 delete vertexp1n1p2;
986
987 }
988
989 postrack2 = 0;
990 delete vertexp2n1;
991
992 } // end 2nd loop on positive tracks
993
994 twoTrackArray2->Clear();
995
996 // 2nd LOOP ON NEGATIVE TRACKS (for 3 prong -+-)
997 for(iTrkN2=iTrkN1+1; iTrkN2<nSeleTrks; iTrkN2++) {
998
999 if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
1000
1001 if(iTrkN2%1==0) AliDebug(1,Form(" 2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));
1002
1003 // get track from tracks array
1004 negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
1005
1006 if(negtrack2->Charge()>0) continue;
1007
1008 if(!TESTBIT(seleFlags[iTrkN2],kBitDispl)) continue;
1009
1010 if(fMixEvent) {
1011 if(evtNumber[iTrkP1]==evtNumber[iTrkN2] ||
1012 evtNumber[iTrkN1]==evtNumber[iTrkN2] ||
1013 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
1014 }
1015
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;
1019 } else { // not ok
1020 continue;
1021 }
1022 } else { // normal triplet (-+-)
1023 isLikeSign3Prong=kFALSE;
1024 }
1025
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());
1031
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; }
1036
1037 threeTrackArray->AddAt(negtrack1,0);
1038 threeTrackArray->AddAt(postrack1,1);
1039 threeTrackArray->AddAt(negtrack2,2);
1040
1041 // check invariant mass cuts for D+,Ds,Lc
1042 massCutOK=kTRUE;
1043 if(fMassCutBeforeVertexing && f3Prong)
1044 massCutOK = SelectInvMass(threeTrackArray);
1045
1046 if(!massCutOK) {
1047 threeTrackArray->Clear();
1048 negtrack2=0;
1049 continue;
1050 }
1051
1052 // Vertexing
1053 twoTrackArray2->AddAt(postrack1,0);
1054 twoTrackArray2->AddAt(negtrack2,1);
1055
1056 AliAODVertex *vertexp1n2 = ReconstructSecondaryVertex(twoTrackArray2,dispersion);
1057 if(!vertexp1n2) {
1058 twoTrackArray2->Clear();
1059 negtrack2=0;
1060 continue;
1061 }
1062
1063 if(f3Prong) {
1064 AliAODVertex* secVert3PrAOD = ReconstructSecondaryVertex(threeTrackArray,dispersion);
1065 io3Prong = Make3Prong(threeTrackArray,event,secVert3PrAOD,dispersion,vertexp1n1,vertexp1n2,dcap1n1,dcap1n2,dcan1n2,ok3Prong);
1066 if(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);
1078 }
1079 }
1080 if(io3Prong) {delete io3Prong; io3Prong=NULL;}
1081 if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
1082 }
1083 threeTrackArray->Clear();
1084 negtrack2 = 0;
1085 delete vertexp1n2;
1086
1087 } // end 2nd loop on negative tracks
1088
1089 twoTrackArray2->Clear();
1090
1091 negtrack1 = 0;
1092 delete vertexp1n1;
1093 } // end 1st loop on negative tracks
1094
1095 postrack1 = 0;
1096 } // end 1st loop on positive tracks
1097
1098
1099 AliDebug(1,Form(" Total HF vertices in event = %d;",
1100 (Int_t)aodVerticesHFTClArr->GetEntriesFast()));
1101 if(fD0toKpi) {
1102 AliDebug(1,Form(" D0->Kpi in event = %d;",
1103 (Int_t)aodD0toKpiTClArr->GetEntriesFast()));
1104 }
1105 if(fJPSItoEle) {
1106 AliDebug(1,Form(" JPSI->ee in event = %d;",
1107 (Int_t)aodJPSItoEleTClArr->GetEntriesFast()));
1108 }
1109 if(f3Prong) {
1110 AliDebug(1,Form(" Charm->3Prong in event = %d;",
1111 (Int_t)aodCharm3ProngTClArr->GetEntriesFast()));
1112 }
1113 if(f4Prong) {
1114 AliDebug(1,Form(" Charm->4Prong in event = %d;\n",
1115 (Int_t)aodCharm4ProngTClArr->GetEntriesFast()));
1116 }
1117 if(fDstar) {
1118 AliDebug(1,Form(" D*->D0pi in event = %d;\n",
1119 (Int_t)aodDstarTClArr->GetEntriesFast()));
1120 }
1121 if(fCascades){
1122 AliDebug(1,Form(" cascades -> v0 + track in event = %d;\n",
1123 (Int_t)aodCascadesTClArr->GetEntriesFast()));
1124 }
1125 if(fLikeSign) {
1126 AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
1127 (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
1128 }
1129 if(fLikeSign && f3Prong) {
1130 AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
1131 (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
1132 }
1133
1134
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;}
1144
1145 if(fInputAOD) {
1146 seleTrksArray.Delete();
1147 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
1148 }
1149
1150 return;
1151}
1152//----------------------------------------------------------------------------
1153void AliAnalysisVertexingHF::AddRefs(AliAODVertex *v,AliAODRecoDecayHF *rd,
1154 const AliVEvent *event,
1155 const TObjArray *trkArray) const
1156{
1157 // Add the AOD tracks as daughters of the vertex (TRef)
1158 // Also add the references to the primary vertex and to the cuts
1159
1160 if(fInputAOD) {
1161 AddDaughterRefs(v,event,trkArray);
1162 rd->SetPrimaryVtxRef((AliAODVertex*)event->GetPrimaryVertex());
1163 }
1164
1165 /*
1166 rd->SetListOfCutsRef((TList*)fListOfCuts);
1167 //fListOfCuts->Print();
1168 cout<<fListOfCuts<<endl;
1169 TList *l=(TList*)rd->GetListOfCuts();
1170 cout<<l<<endl;
1171 if(l) {l->Print(); }else{printf("error\n");}
1172 */
1173
1174 return;
1175}
1176//----------------------------------------------------------------------------
1177void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,
1178 const AliVEvent *event,
1179 const TObjArray *trkArray) const
1180{
1181 // Add the AOD tracks as daughters of the vertex (TRef)
1182
1183 Int_t nDg = v->GetNDaughters();
1184 TObject *dg = 0;
1185 if(nDg) dg = v->GetDaughter(0);
1186
1187 if(dg) return; // daughters already added
1188
1189 Int_t nTrks = trkArray->GetEntriesFast();
1190
1191 AliExternalTrackParam *track = 0;
1192 AliAODTrack *aodTrack = 0;
1193 Int_t id;
1194
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);
1202 }
1203
1204 return;
1205}
1206//----------------------------------------------------------------------------
1207AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1208 TObjArray *twoTrackArray,AliVEvent *event,
1209 AliAODVertex *secVert,
1210 AliAODRecoDecayHF2Prong *rd2Prong,
1211 Double_t dca,
1212 Bool_t &okDstar) const
1213{
1214 // Make the cascade as a 2Prong decay and check if it passes Dstar
1215 // reconstruction cuts
1216
1217 okDstar = kFALSE;
1218
1219 Bool_t dummy1,dummy2,dummy3;
1220
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;
1227
1228 // charge
1229 AliESDtrack *trackPi = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1230 theCascade->SetCharge(trackPi->Charge());
1231
1232 //--- selection cuts
1233 //
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);
1243 }
1244 // select D*->D0pi
1245 if(fDstar) {
1246 okDstar = (Bool_t)fCutsDStartoKpipi->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
1247 }
1248 tmpCascade->GetSecondaryVtx()->RemoveDaughters();
1249 tmpCascade->UnsetOwnPrimaryVtx();
1250 delete tmpCascade; tmpCascade=NULL;
1251 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1252 rd2Prong->UnsetOwnPrimaryVtx();
1253 }
1254 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
1255 //---
1256
1257 return theCascade;
1258}
1259
1260
1261//----------------------------------------------------------------------------
1262AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
1263 TObjArray *twoTrackArray,AliVEvent *event,
1264 AliAODVertex *secVert,
1265 AliAODv0 *v0,
1266 Double_t dca,
1267 Bool_t &okCascades) const
1268{
1269 //
1270 // Make the cascade as a 2Prong decay and check if it passes
1271 // cascades reconstruction cuts
1272
1273 // AliDebug(2,Form(" building the cascade"));
1274 okCascades= kFALSE;
1275 Bool_t dummy1,dummy2,dummy3;
1276
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;
1283
1284 // bachelor track and charge
1285 AliESDtrack *trackBachelor = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
1286 theCascade->SetCharge(trackBachelor->Charge());
1287
1288 //--- selection cuts
1289 //
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);
1299 }
1300
1301 // select Cascades
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;}
1308 }
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;}
1314 //---
1315
1316 return theCascade;
1317}
1318
1319//-----------------------------------------------------------------------------
1320AliAODRecoDecayHF2Prong *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
1325{
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)
1329
1330 okD0=kFALSE; okJPSI=kFALSE; okD0fromDstar=kFALSE;
1331
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);
1335
1336 // propagate tracks to secondary vertex, to compute inv. mass
1337 postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1338 negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
1339
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];
1345
1346
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;
1352 if(!okMassCut) {
1353 AliDebug(2," candidate didn't pass mass cut");
1354 return 0x0;
1355 }
1356 // primary vertex to be used by this candidate
1357 AliAODVertex *primVertexAOD = PrimaryVertex(twoTrackArray,event);
1358 if(!primVertexAOD) return 0x0;
1359
1360
1361 Double_t d0z0[2],covd0z0[3];
1362 postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1363 d0[0] = d0z0[0];
1364 d0err[0] = TMath::Sqrt(covd0z0[0]);
1365 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1366 d0[1] = d0z0[0];
1367 d0err[1] = TMath::Sqrt(covd0z0[0]);
1368
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;
1375
1376
1377 if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar
1378 // select D0->Kpi
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);
1387 }
1388
1389
1390 // remove the primary vertex (was used only for selection)
1391 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1392 the2Prong->UnsetOwnPrimaryVtx();
1393 }
1394
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];
1404 }
1405 the2Prong->SetPID(2,esdpid);
1406
1407 return the2Prong;
1408}
1409//----------------------------------------------------------------------------
1410AliAODRecoDecayHF3Prong* 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
1416{
1417 // Make 3Prong candidates and check if they pass Dplus or Ds or Lambdac
1418 // reconstruction cuts
1419 // E.Bruna, F.Prino
1420
1421
1422 ok3Prong=kFALSE;
1423 if(!secVert || !vertexp1n1 || !vertexp2n1) return 0x0;
1424
1425 Double_t px[3],py[3],pz[3],d0[3],d0err[3];
1426 Double_t momentum[3];
1427
1428
1429 AliESDtrack *postrack1 = (AliESDtrack*)threeTrackArray->UncheckedAt(0);
1430 AliESDtrack *negtrack = (AliESDtrack*)threeTrackArray->UncheckedAt(1);
1431 AliESDtrack *postrack2 = (AliESDtrack*)threeTrackArray->UncheckedAt(2);
1432
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];
1442
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;
1447 if(!okMassCut) {
1448 AliDebug(2," candidate didn't pass mass cut");
1449 return 0x0;
1450 }
1451
1452 // primary vertex to be used by this candidate
1453 AliAODVertex *primVertexAOD = PrimaryVertex(threeTrackArray,event);
1454 if(!primVertexAOD) return 0x0;
1455
1456 Double_t d0z0[2],covd0z0[3];
1457 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1458 d0[0]=d0z0[0];
1459 d0err[0] = TMath::Sqrt(covd0z0[0]);
1460 negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1461 d0[1]=d0z0[0];
1462 d0err[1] = TMath::Sqrt(covd0z0[0]);
1463 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1464 d0[2]=d0z0[0];
1465 d0err[2] = TMath::Sqrt(covd0z0[0]);
1466
1467
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);
1478
1479 delete primVertexAOD; primVertexAOD=NULL;
1480
1481 // select D+->Kpipi, Ds->KKpi, Lc->pKpi
1482 if(f3Prong) {
1483 ok3Prong = kFALSE;
1484
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;
1488
1489 }
1490 //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
1491
1492 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1493 the3Prong->UnsetOwnPrimaryVtx();
1494 }
1495
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);
1503
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];
1509 }
1510 the3Prong->SetPID(3,esdpid);
1511
1512 return the3Prong;
1513}
1514//----------------------------------------------------------------------------
1515AliAODRecoDecayHF4Prong* 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
1523{
1524 // Make 4Prong candidates and check if they pass D0toKpipipi
1525 // reconstruction cuts
1526 // G.E.Bruno, R.Romita
1527
1528 ok4Prong=kFALSE;
1529 if(!secVert || !vertexp1n1 || !vertexp1n1p2) return 0x0;
1530
1531 Double_t px[4],py[4],pz[4],d0[4],d0err[4];//d0z[3];
1532
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);
1537
1538 postrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1539 negtrack1->PropagateToDCA(secVert,fBzkG,kVeryBig);
1540 postrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1541 negtrack2->PropagateToDCA(secVert,fBzkG,kVeryBig);
1542
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];
1552
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;
1558 }
1559 if(!okMassCut) {
1560 //if(fDebug) printf(" candidate didn't pass mass cut\n");
1561 //printf(" candidate didn't pass mass cut\n");
1562 return 0x0;
1563 }
1564
1565 // primary vertex to be used by this candidate
1566 AliAODVertex *primVertexAOD = PrimaryVertex(fourTrackArray,event);
1567 if(!primVertexAOD) return 0x0;
1568
1569 Double_t d0z0[2],covd0z0[3];
1570 postrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1571 d0[0]=d0z0[0];
1572 d0err[0] = TMath::Sqrt(covd0z0[0]);
1573 negtrack1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1574 d0[1]=d0z0[0];
1575 d0err[1] = TMath::Sqrt(covd0z0[0]);
1576 postrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1577 d0[2]=d0z0[0];
1578 d0err[2] = TMath::Sqrt(covd0z0[0]);
1579 negtrack2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1580 d0[3]=d0z0[0];
1581 d0err[3] = TMath::Sqrt(covd0z0[0]);
1582
1583
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]));
1590 Short_t charge=0;
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);
1595
1596 delete primVertexAOD; primVertexAOD=NULL;
1597
1598 ok4Prong=(Bool_t)fCutsD0toKpipipi->IsSelected(the4Prong,AliRDHFCuts::kCandidate);
1599
1600
1601 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx && !fMixEvent) {
1602 the4Prong->UnsetOwnPrimaryVtx();
1603 }
1604
1605
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);
1615
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];
1622 }
1623 the4Prong->SetPID(4,esdpid);
1624
1625 return the4Prong;
1626}
1627//-----------------------------------------------------------------------------
1628AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
1629 AliVEvent *event) const
1630{
1631 // Returns primary vertex to be used for this candidate
1632
1633 AliESDVertex *vertexESD = 0;
1634 AliAODVertex *vertexAOD = 0;
1635
1636
1637 if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1638 // primary vertex from the input event
1639
1640 vertexESD = new AliESDVertex(*fV1);
1641
1642 } else {
1643 // primary vertex specific to this candidate
1644
1645 Int_t nTrks = trkArray->GetEntriesFast();
1646 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1647
1648 if(fRecoPrimVtxSkippingTrks) {
1649 // recalculating the vertex
1650
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();
1661 }
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();
1668 if(id<0) continue;
1669 skipped[nTrksToSkip++] = id;
1670 }
1671 // TEMPORARY FIX
1672 // For AOD, skip also tracks without covariance matrix
1673 if(fInputAOD) {
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();
1679 if(id<0) continue;
1680 skipped[nTrksToSkip++] = id;
1681 }
1682 }
1683 }
1684 //
1685 vertexer->SetSkipTracks(nTrksToSkip,skipped);
1686 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1687
1688 } else if(fRmTrksFromPrimVtx) {
1689 // removing the prongs tracks
1690
1691 TObjArray rmArray(nTrks);
1692 UShort_t *rmId = new UShort_t[nTrks];
1693 AliESDtrack *esdTrack = 0;
1694 AliESDtrack *t = 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();
1701 } else {
1702 rmId[i]=9999;
1703 }
1704 }
1705 Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
1706 vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1707 delete [] rmId; rmId=NULL;
1708 rmArray.Delete();
1709
1710 }
1711
1712 if(!vertexESD) return vertexAOD;
1713 if(vertexESD->GetNContributors()<=0) {
1714 AliDebug(2,"vertexing failed");
1715 delete vertexESD; vertexESD=NULL;
1716 return vertexAOD;
1717 }
1718
1719 delete vertexer; vertexer=NULL;
1720
1721 }
1722
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;
1729
1730 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1731
1732 return vertexAOD;
1733}
1734//-----------------------------------------------------------------------------
1735void AliAnalysisVertexingHF::PrintStatus() const {
1736 // Print parameters being used
1737
1738 //printf("Preselections:\n");
1739 // fTrackFilter->Dump();
1740 if(fSecVtxWithKF) {
1741 printf("Secondary vertex with Kalman filter package (AliKFParticle)\n");
1742 } else {
1743 printf("Secondary vertex with AliVertexerTracks\n");
1744 }
1745 if(fRecoPrimVtxSkippingTrks) printf("RecoPrimVtxSkippingTrks\n");
1746 if(fRmTrksFromPrimVtx) printf("RmTrksFromPrimVtx\n");
1747 if(fD0toKpi) {
1748 printf("Reconstruct D0->Kpi candidates with cuts:\n");
1749 if(fCutsD0toKpi) fCutsD0toKpi->PrintAll();
1750 }
1751 if(fDstar) {
1752 printf("Reconstruct D*->D0pi candidates with cuts:\n");
1753 if(fFindVertexForDstar) {
1754 printf(" Reconstruct a secondary vertex for the D*\n");
1755 } else {
1756 printf(" Assume the D* comes from the primary vertex\n");
1757 }
1758 if(fCutsDStartoKpipi) fCutsDStartoKpipi->PrintAll();
1759 }
1760 if(fJPSItoEle) {
1761 printf("Reconstruct J/psi from B candidates with cuts:\n");
1762 if(fCutsJpsitoee) fCutsJpsitoee->PrintAll();
1763 }
1764 if(f3Prong) {
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();
1772 }
1773 if(f4Prong) {
1774 printf("Reconstruct 4 prong candidates.\n");
1775 printf(" D0->Kpipipi cuts:\n");
1776 if(fCutsD0toKpipipi) fCutsD0toKpipipi->PrintAll();
1777 }
1778 if(fCascades) {
1779 printf("Reconstruct cascades candidates formed with v0s.\n");
1780 printf(" Lc -> k0s P & Lc -> L Pi cuts:\n");
1781 if(fCutsLctoV0) fCutsLctoV0->PrintAll();
1782 }
1783
1784 return;
1785}
1786//-----------------------------------------------------------------------------
1787AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkArray,
1788 Double_t &dispersion,Bool_t useTRefArray) const
1789{
1790 // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
1791
1792 AliESDVertex *vertexESD = 0;
1793 AliAODVertex *vertexAOD = 0;
1794
1795 if(!fSecVtxWithKF) { // AliVertexerTracks
1796
1797 AliVertexerTracks *vertexer = new AliVertexerTracks(fBzkG);
1798 vertexer->SetVtxStart(fV1);
1799 vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1800 delete vertexer; vertexer=NULL;
1801
1802 if(!vertexESD) return vertexAOD;
1803
1804 if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
1805 AliDebug(2,"vertexing failed");
1806 delete vertexESD; vertexESD=NULL;
1807 return vertexAOD;
1808 }
1809
1810 } else { // Kalman Filter vertexer (AliKFParticle)
1811
1812 AliKFParticle::SetField(fBzkG);
1813
1814 AliKFVertex vertexKF;
1815
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);
1821 }
1822 vertexESD = new AliESDVertex(vertexKF.Parameters(),
1823 vertexKF.CovarianceMatrix(),
1824 vertexKF.GetChi2(),
1825 vertexKF.GetNContributors());
1826
1827 }
1828
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;
1836
1837 Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
1838 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
1839
1840 return vertexAOD;
1841}
1842//-----------------------------------------------------------------------------
1843Bool_t AliAnalysisVertexingHF::SelectInvMass(TObjArray *trkArray) const {
1844 // Invariant mass cut on tracks
1845
1846 Int_t nProngs=trkArray->GetEntriesFast();
1847 Int_t retval=kFALSE;
1848
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;
1856
1857 switch(nProngs) {
1858 case 3:
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);
1870 break;
1871 case 4:
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);
1886 break;
1887 default:
1888 printf("SelectInvMass(): wrong decay selection\n");
1889 break;
1890 }
1891
1892 return retval;
1893}
1894//-----------------------------------------------------------------------------
1895Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
1896 Int_t nprongs,
1897 Double_t *px,
1898 Double_t *py,
1899 Double_t *pz) const {
1900 // Check invariant mass cut
1901
1902 UInt_t pdg2[2],pdg3[3],pdg4[4];
1903 Double_t mPDG,minv;
1904
1905 Bool_t retval=kFALSE;
1906 switch (decay)
1907 {
1908 case 0: // D0->Kpi
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;
1917 break;
1918 case 1: // JPSI->ee
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;
1924 break;
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;
1931 // Ds+->KKpi
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;
1939 // Lc->pKpi
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;
1947 break;
1948 case 3: // D*->D0pi
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;
1954 break;
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;
1973 break;
1974 default:
1975 printf("SelectInvMass(): wrong decay selection\n");
1976 break;
1977 }
1978
1979 return retval;
1980}
1981//-----------------------------------------------------------------------------
1982void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
1983 TObjArray &seleTrksArray,Int_t &nSeleTrks,
1984 UChar_t *seleFlags,Int_t *evtNumber)
1985{
1986 // Apply single-track preselection.
1987 // Fill a TObjArray with selected tracks (for displaced vertices or
1988 // soft pion from D*). Selection flag stored in seleFlags.
1989 // Create the AliESDVertex object (convert from AliAODVertex if necessary)
1990 // In case of AOD input, also fill fAODMap for track index<->ID
1991
1992 const AliVVertex *vprimary = event->GetPrimaryVertex();
1993
1994 if(fV1) { delete fV1; fV1=NULL; }
1995 if(fAODMap) { delete fAODMap; fAODMap=NULL; }
1996
1997 Int_t nindices=0;
1998 UShort_t *indices = 0;
1999 Double_t pos[3],cov[6];
2000
2001 if(!fInputAOD) { // ESD
2002 fV1 = new AliESDVertex(*((AliESDVertex*)vprimary));
2003 } else { // AOD
2004 vprimary->GetXYZ(pos);
2005 vprimary->GetCovarianceMatrix(cov);
2006 fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
2007 indices = new UShort_t[event->GetNumberOfTracks()];
2008 fAODMapSize = 100000;
2009 fAODMap = new Int_t[fAODMapSize];
2010 }
2011
2012
2013 Int_t entries = (Int_t)event->GetNumberOfTracks();
2014 Bool_t okDisplaced=kFALSE,okSoftPi=kFALSE;
2015 nSeleTrks=0;
2016
2017 // transfer ITS tracks from event to arrays
2018 for(Int_t i=0; i<entries; i++) {
2019 AliVTrack *track;
2020 track = (AliVTrack*)event->GetTrack(i);
2021
2022 // skip pure ITS SA tracks
2023 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
2024
2025 // TEMPORARY: check that the cov matrix is there
2026 Double_t covtest[21];
2027 if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
2028 //
2029
2030 if(fInputAOD) {
2031 AliAODTrack *aodt = (AliAODTrack*)track;
2032 if(aodt->GetUsedForPrimVtxFit()) {
2033 indices[nindices]=aodt->GetID(); nindices++;
2034 }
2035 fAODMap[(Int_t)aodt->GetID()] = i;
2036 }
2037
2038 AliESDtrack *esdt = 0;
2039
2040 if(!fInputAOD) {
2041 esdt = (AliESDtrack*)track;
2042 } else {
2043 esdt = new AliESDtrack(track);
2044 }
2045
2046 // single track selection
2047 okDisplaced=kFALSE; okSoftPi=kFALSE;
2048 if(fMixEvent){
2049 evtNumber[i]=((AliMixedEvent*)event)->EventIndex(i);
2050 const AliVVertex* eventVtx=((AliMixedEvent*)event)->GetEventVertex(i);
2051 Double_t vtxPos[3],primPos[3],primCov[6],trasl[3];
2052 eventVtx->GetXYZ(vtxPos);
2053 vprimary->GetXYZ(primPos);
2054 eventVtx->GetCovarianceMatrix(primCov);
2055 for(Int_t ind=0;ind<3;ind++){
2056 trasl[ind]=vtxPos[ind]-primPos[ind];
2057 }
2058
2059 Bool_t isTransl=esdt->Translate(trasl,primCov);
2060 if(!isTransl) {
2061 delete esdt;
2062 esdt = NULL;
2063 continue;
2064 }
2065 }
2066
2067 if(SingleTrkCuts(esdt,okDisplaced,okSoftPi)) {
2068 seleTrksArray.AddLast(esdt);
2069 seleFlags[nSeleTrks]=0;
2070 if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
2071 if(okSoftPi) SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
2072 nSeleTrks++;
2073 } else {
2074 if(fInputAOD) delete esdt;
2075 esdt = NULL;
2076 continue;
2077 }
2078
2079 } // end loop on tracks
2080
2081 // primary vertex from AOD
2082 if(fInputAOD) {
2083 delete fV1; fV1=NULL;
2084 vprimary->GetXYZ(pos);
2085 vprimary->GetCovarianceMatrix(cov);
2086 Double_t chi2toNDF = vprimary->GetChi2perNDF();
2087 Int_t ncontr=nindices;
2088 if(!strcmp(vprimary->GetTitle(),"VertexerTracksWithContraint")) ncontr += 1;
2089 Double_t chi2=chi2toNDF*(2.*(Double_t)ncontr-3.);
2090 fV1 = new AliESDVertex(pos,cov,chi2,ncontr,vprimary->GetName());
2091 fV1->SetTitle(vprimary->GetTitle());
2092 fV1->SetIndices(nindices,indices);
2093 }
2094 if(indices) { delete [] indices; indices=NULL; }
2095
2096
2097 return;
2098}
2099//-----------------------------------------------------------------------------
2100Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
2101 Bool_t &okDisplaced,Bool_t &okSoftPi) const
2102{
2103 // Check if track passes some kinematical cuts
2104
2105 // this is needed to store the impact parameters
2106 trk->RelateToVertex(fV1,fBzkG,kVeryBig);
2107
2108 UInt_t selectInfo;
2109 //
2110 // Track selection, displaced tracks
2111 selectInfo = 0;
2112 if(fTrackFilter) {
2113 selectInfo = fTrackFilter->IsSelected(trk);
2114 }
2115 if(selectInfo) okDisplaced=kTRUE;
2116 // Track selection, soft pions
2117 selectInfo = 0;
2118 if(fDstar && fTrackFilterSoftPi) {
2119 selectInfo = fTrackFilterSoftPi->IsSelected(trk);
2120 }
2121 if(selectInfo) okSoftPi=kTRUE;
2122
2123 if(okDisplaced || okSoftPi) return kTRUE;
2124
2125 return kFALSE;
2126}
2127
2128
2129//-----------------------------------------------------------------------------
2130AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArray *twoTrackArrayV0){
2131 //
2132 // Transform ESDv0 to AODv0
2133 //
2134 // this function takes the ESDv0 vertex, computes the DCA variables from the ESDv0
2135 // and creates an AODv0 out of them
2136 //
2137 Double_t vertex[3]; esdV0->GetXYZ(vertex[0],vertex[1],vertex[2]);
2138 AliAODVertex *vertexV0 = new AliAODVertex(vertex,esdV0->GetChi2V0(),AliAODVertex::kV0,2);
2139
2140 // create the v0 neutral track to compute the DCA to the primary vertex
2141 Double_t xyz[3], pxpypz[3];
2142 esdV0->XvYvZv(xyz);
2143 esdV0->PxPyPz(pxpypz);
2144 Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
2145 AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
2146 if(!trackesdV0) return 0;
2147 Double_t d0z0[2],covd0z0[3];
2148 AliAODVertex *primVertexAOD = PrimaryVertex();
2149 trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2150 Double_t dcaV0ToPrimVertex = TMath::Sqrt(covd0z0[0]);
2151 // get the v0 daughters to compute their DCA to the v0 vertex and get their momentum
2152 Double_t dcaV0DaughterToPrimVertex[2];
2153 AliExternalTrackParam *posV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(0);
2154 AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
2155 if( !posV0track || !negV0track) {
2156 if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2157 return 0;
2158 }
2159 posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2160 // if ( covd0z0[0]<=0.) dcaV0DaughterToPrimVertex[0] = 0;
2161 // else
2162 dcaV0DaughterToPrimVertex[0] = TMath::Sqrt(covd0z0[0]);
2163 negV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
2164 // if ( covd0z0[0]<=0.)dcaV0DaughterToPrimVertex[1] = 0;
2165 // else
2166 dcaV0DaughterToPrimVertex[1] = TMath::Sqrt(covd0z0[0]);
2167 Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
2168 Double_t pmom[3],nmom[3];
2169 esdV0->GetNPxPyPz(nmom[0],nmom[1],nmom[2]);
2170 esdV0->GetPPxPyPz(pmom[0],pmom[1],pmom[2]);
2171
2172 AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
2173 aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
2174
2175 if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
2176 if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
2177
2178 return aodV0;
2179}
2180//-----------------------------------------------------------------------------