AliPhysics  master (3d17d9d)
AliAnalysisTaskCombinHF.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2018, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 /* $Id: $ */
17 
18 //*************************************************************************
19 // Class AliAnalysisTaskCombinHF
20 // AliAnalysisTaskSE to build D meson candidates by combining tracks
21 // background is computed LS and track rotations is
22 // Authors: F. Prino, A. Rossi
24 
25 #include <TList.h>
26 #include <TH1F.h>
27 #include <TDatabasePDG.h>
28 #include <TH2F.h>
29 #include <TH3F.h>
30 #include <THnSparse.h>
31 
32 #include "AliAnalysisManager.h"
33 #include "AliInputEventHandler.h"
34 #include "AliPIDResponse.h"
35 #include "AliAODHandler.h"
36 #include "AliAODEvent.h"
37 #include "AliAODMCParticle.h"
38 #include "AliAODMCHeader.h"
39 #include "AliAODVertex.h"
40 #include "AliAODTrack.h"
42 #include "AliVertexingHFUtils.h"
43 #include "AliGenEventHeader.h"
44 #include "AliGenCocktailEventHeader.h"
45 #include "AliGenPythiaEventHeader.h"
47 
49 ClassImp(AliAnalysisTaskCombinHF);
51 
52 //________________________________________________________________________
55  fOutput(0x0),
56  fListCuts(0x0),
57  fHistNEvents(0x0),
58  fHistEventMultCent(0x0),
59  fHistEventMultCentEvSel(0x0),
60  fHistEventMultZv(0x0),
61  fHistEventMultZvEvSel(0x0),
62  fHistXsecVsPtHard(0x0),
63  fHistTrackStatus(0x0),
64  fHistTrackEtaMultZv(0x0),
65  fHistSelTrackPhiPt(0x0),
66  fHistCheckOrigin(0x0),
67  fHistCheckOriginRecoD(0x0),
68  fHistCheckOriginRecoVsGen(0x0),
69  fHistCheckDecChan(0x0),
70  fHistCheckDecChanAcc(0x0),
71  fPtVsYVsMultGenPrompt(0x0),
72  fPtVsYVsMultGenLargeAccPrompt(0x0),
73  fPtVsYVsMultGenLimAccPrompt(0x0),
74  fPtVsYVsMultGenAccPrompt(0x0),
75  fPtVsYVsMultGenAccEvSelPrompt(0x0),
76  fPtVsYVsMultRecoPrompt(0x0),
77  fPtVsYVsMultGenFeeddw(0x0),
78  fPtVsYVsMultGenLargeAccFeeddw(0x0),
79  fPtVsYVsMultGenLimAccFeeddw(0x0),
80  fPtVsYVsMultGenAccFeeddw(0x0),
81  fPtVsYVsMultGenAccEvSelFeeddw(0x0),
82  fPtVsYVsMultRecoFeeddw(0x0),
83  fPtVsYVsPtBGenFeeddw(0x0),
84  fPtVsYVsPtBGenLargeAccFeeddw(0x0),
85  fPtVsYVsPtBGenLimAccFeeddw(0x0),
86  fPtVsYVsPtBGenAccFeeddw(0x0),
87  fPtVsYVsPtBGenAccEvSelFeeddw(0x0),
88  fPtVsYVsPtBRecoFeeddw(0x0),
89  fMassVsPtVsY(0x0),
90  fMassVsPtVsYRot(0x0),
91  fMassVsPtVsYLSpp(0x0),
92  fMassVsPtVsYLSmm(0x0),
93  fMassVsPtVsYSig(0x0),
94  fMassVsPtVsYRefl(0x0),
95  fMassVsPtVsYBkg(0x0),
96  fBMohterPtGen(0x0),
97  fNSelected(0x0),
98  fNormRotated(0x0),
99  fDeltaMass(0x0),
100  fDeltaMassFullAnalysis(0x0),
101  fMassVsPtVsYME(0x0),
102  fMassVsPtVsYMELSpp(0x0),
103  fMassVsPtVsYMELSmm(0x0),
104  fEventsPerPool(0x0),
105  fMixingsPerPool(0x0),
106  fMassVsPtVsCosthSt(0x0),
107  fMassVsPtVsCosthStRot(0x0),
108  fMassVsPtVsCosthStLSpp(0x0),
109  fMassVsPtVsCosthStLSmm(0x0),
110  fMassVsPtVsCosthStSig(0x0),
111  fMassVsPtVsCosthStRefl(0x0),
112  fMassVsPtVsCosthStBkg(0x0),
113  fMassVsPtVsCosthStME(0x0),
114  fMassVsPtVsCosthStMELSpp(0x0),
115  fMassVsPtVsCosthStMELSmm(0x0),
116  fHistonSigmaTPCPion(0x0),
117  fHistonSigmaTPCPionGoodTOF(0x0),
118  fHistonSigmaTOFPion(0x0),
119  fHistonSigmaTPCKaon(0x0),
120  fHistonSigmaTPCKaonGoodTOF(0x0),
121  fHistonSigmaTOFKaon(0x0),
122  fHistoPtKPtPiPtD(0x0),
123  fHistoPtKPtPiPtDSig(0x0),
124  fFilterMask(BIT(4)),
125  fTrackCutsAll(0x0),
126  fTrackCutsPion(0x0),
127  fTrackCutsKaon(0x0),
128  fCutTPCSignalN(0),
129  fFillHistosVsCosThetaStar(kFALSE),
130  fApplyCutCosThetaStar(kFALSE),
131  fCutCosThetaStar(999.),
132  fPhiMassCut(99999.),
133  fCutCos3PiKPhiRFrame(-1.1),
134  fCutCosPiDsLabFrame(1.1),
135  fPidHF(0x0),
136  fAnalysisCuts(0x0),
137  fMinMass(1.720),
138  fMaxMass(2.150),
139  fMaxPt(10.),
140  fPtBinWidth(0.5),
141  fEtaAccCut(0.9),
142  fPtAccCut(0.1),
143  fNRotations(9),
144  fMinAngleForRot(5*TMath::Pi()/6),
145  fMaxAngleForRot(7*TMath::Pi()/6),
146  fNRotations3(9),
147  fMinAngleForRot3(2*TMath::Pi()/6),
148  fMaxAngleForRot3(4*TMath::Pi()/6),
149  fCounter(0x0),
150  fMeson(kDzero),
151  fReadMC(kFALSE),
152  fEnforceMBTrigMaskInMC(kTRUE),
153  fGoUpToQuark(kTRUE),
154  fFullAnalysis(0),
155  fSignalOnlyMC(kFALSE),
156  fSelectPtHardRange(kFALSE),
157  fMinPtHard(0.),
158  fMaxPtHard(999999.),
159  fPIDstrategy(knSigma),
160  fmaxPforIDPion(0.8),
161  fmaxPforIDKaon(2.),
162  fKeepNegID(kFALSE),
163  fPIDselCaseZero(0),
164  fBayesThresKaon(0.4),
165  fBayesThresPion(0.4),
166  fDoEventMixing(1),
167  fNumberOfEventsForMixing(20),
168  fMaxzVertDistForMix(5.),
169  fMaxMultDiffForMix(5.),
170  fNzVertPools(1),
171  fNzVertPoolsLimSize(2),
172  fzVertPoolLims(0x0),
173  fNMultPools(1),
174  fNMultPoolsLimSize(2),
175  fMultPoolLims(0x0),
176  fNOfPools(1),
177  fEventBuffer(0x0),
178  fEventInfo(new TObjString("")),
179  fVtxZ(0),
180  fMultiplicity(0),
181  fNumOfMultBins(200),
182  fMinMultiplicity(-0.5),
183  fMaxMultiplicity(199.5),
184  fKaonTracks(0x0),
185  fPionTracks(0x0)
186 {
188 }
189 
190 //________________________________________________________________________
192  AliAnalysisTaskSE("DmesonCombin"),
193  fOutput(0x0),
194  fListCuts(0x0),
195  fHistNEvents(0x0),
196  fHistEventMultCent(0x0),
198  fHistEventMultZv(0x0),
200  fHistXsecVsPtHard(0x0),
201  fHistTrackStatus(0x0),
202  fHistTrackEtaMultZv(0x0),
203  fHistSelTrackPhiPt(0x0),
204  fHistCheckOrigin(0x0),
207  fHistCheckDecChan(0x0),
227  fMassVsPtVsY(0x0),
228  fMassVsPtVsYRot(0x0),
229  fMassVsPtVsYLSpp(0x0),
230  fMassVsPtVsYLSmm(0x0),
231  fMassVsPtVsYSig(0x0),
232  fMassVsPtVsYRefl(0x0),
233  fMassVsPtVsYBkg(0x0),
234  fBMohterPtGen(0x0),
235  fNSelected(0x0),
236  fNormRotated(0x0),
237  fDeltaMass(0x0),
239  fMassVsPtVsYME(0x0),
240  fMassVsPtVsYMELSpp(0x0),
241  fMassVsPtVsYMELSmm(0x0),
242  fEventsPerPool(0x0),
243  fMixingsPerPool(0x0),
244  fMassVsPtVsCosthSt(0x0),
254  fHistonSigmaTPCPion(0x0),
256  fHistonSigmaTOFPion(0x0),
257  fHistonSigmaTPCKaon(0x0),
259  fHistonSigmaTOFKaon(0x0),
260  fHistoPtKPtPiPtD(0x0),
261  fHistoPtKPtPiPtDSig(0x0),
262  fFilterMask(BIT(4)),
263  fTrackCutsAll(0x0),
264  fTrackCutsPion(0x0),
265  fTrackCutsKaon(0x0),
266  fCutTPCSignalN(0),
268  fApplyCutCosThetaStar(kFALSE),
269  fCutCosThetaStar(999.),
270  fPhiMassCut(99999.),
272  fCutCosPiDsLabFrame(1.1),
273  fPidHF(0x0),
274  fAnalysisCuts(analysiscuts),
275  fMinMass(1.720),
276  fMaxMass(2.150),
277  fMaxPt(10.),
278  fPtBinWidth(0.5),
279  fEtaAccCut(0.9),
280  fPtAccCut(0.1),
281  fNRotations(9),
282  fMinAngleForRot(5*TMath::Pi()/6),
283  fMaxAngleForRot(7*TMath::Pi()/6),
284  fNRotations3(9),
285  fMinAngleForRot3(2*TMath::Pi()/6),
286  fMaxAngleForRot3(4*TMath::Pi()/6),
287  fCounter(0x0),
288  fMeson(meson),
289  fReadMC(kFALSE),
290  fEnforceMBTrigMaskInMC(kTRUE),
291  fGoUpToQuark(kTRUE),
292  fFullAnalysis(0),
293  fSignalOnlyMC(kFALSE),
294  fSelectPtHardRange(kFALSE),
295  fMinPtHard(0.),
296  fMaxPtHard(999999.),
298  fmaxPforIDPion(0.8),
299  fmaxPforIDKaon(2.),
300  fKeepNegID(kFALSE),
301  fPIDselCaseZero(0),
302  fBayesThresKaon(0.4),
303  fBayesThresPion(0.4),
304  fDoEventMixing(1),
307  fMaxMultDiffForMix(5.),
308  fNzVertPools(1),
310  fzVertPoolLims(0x0),
311  fNMultPools(1),
313  fMultPoolLims(0x0),
314  fNOfPools(1),
315  fEventBuffer(0x0),
316  fEventInfo(new TObjString("")),
317  fVtxZ(0),
318  fMultiplicity(0),
319  fNumOfMultBins(200),
320  fMinMultiplicity(-0.5),
321  fMaxMultiplicity(199.5),
322  fKaonTracks(0x0),
323  fPionTracks(0x0)
324 {
326  DefineOutput(1,TList::Class()); //My private output
327  DefineOutput(2,AliNormalizationCounter::Class());
328  DefineOutput(3,TList::Class());
329 }
330 
331 //________________________________________________________________________
333 {
334  //
336  //
337  if(fOutput && !fOutput->IsOwner()){
338  delete fHistNEvents;
339  delete fHistEventMultCent;
341  delete fHistEventMultZv;
342  delete fHistEventMultZvEvSel;
343  delete fHistXsecVsPtHard;
344  delete fHistTrackStatus;
345  delete fHistTrackEtaMultZv;
346  delete fHistSelTrackPhiPt;
347  delete fHistCheckOrigin;
348  delete fHistCheckOriginRecoD;
350  delete fHistCheckDecChan;
351  delete fHistCheckDecChanAcc;
352  delete fPtVsYVsMultGenPrompt;
357  delete fPtVsYVsMultRecoPrompt;
358  delete fPtVsYVsMultGenFeeddw;
363  delete fPtVsYVsMultRecoFeeddw;
364  delete fPtVsYVsPtBGenFeeddw;
369  delete fPtVsYVsPtBRecoFeeddw;
370  delete fMassVsPtVsY;
371  delete fMassVsPtVsYLSpp;
372  delete fMassVsPtVsYLSmm;
373  delete fMassVsPtVsYRot;
374  delete fMassVsPtVsYSig;
375  delete fMassVsPtVsYRefl;
376  delete fMassVsPtVsYBkg;
377  delete fBMohterPtGen;
378  delete fNSelected;
379  delete fNormRotated;
380  delete fDeltaMass;
381  delete fDeltaMassFullAnalysis;
382  delete fMassVsPtVsYME;
383  delete fMassVsPtVsYMELSpp;
384  delete fMassVsPtVsYMELSmm;
385  delete fMassVsPtVsCosthSt;
386  delete fMassVsPtVsCosthStRot;
387  delete fMassVsPtVsCosthStLSpp;
388  delete fMassVsPtVsCosthStLSmm;
389  delete fMassVsPtVsCosthStSig;
390  delete fMassVsPtVsCosthStRefl;
391  delete fMassVsPtVsCosthStBkg;
392  delete fMassVsPtVsCosthStME;
395  delete fHistonSigmaTPCPion;
397  delete fHistonSigmaTOFPion;
398  delete fHistonSigmaTPCKaon;
400  delete fHistonSigmaTOFKaon;
401  delete fHistoPtKPtPiPtD;
402  delete fHistoPtKPtPiPtDSig;
403  }
404 
405  delete fOutput;
406  if (fListCuts) delete fListCuts;
407  delete fCounter;
408  delete fTrackCutsAll;
409  delete fTrackCutsPion;
410  delete fTrackCutsKaon;
411  delete fAnalysisCuts;
412  if(fKaonTracks) fKaonTracks->Delete();
413  if(fPionTracks) fPionTracks->Delete();
414  delete fKaonTracks;
415  delete fPionTracks;
416 
417  if(fEventBuffer){
418  for(Int_t i=0; i<fNOfPools; i++) delete fEventBuffer[i];
419  delete fEventBuffer;
420  }
421  delete fEventInfo;
422  delete [] fzVertPoolLims;
423  delete [] fMultPoolLims;
424 }
425 
426 //________________________________________________________________________
428 {
430  if(fzVertPoolLims) delete [] fzVertPoolLims;
431  fNzVertPools=nPools;
432  fNzVertPoolsLimSize=nPools+1;
434  for(Int_t ib=0; ib<fNzVertPoolsLimSize; ib++) fzVertPoolLims[ib]=zVertLimits[ib];
435  return;
436 }
437 //________________________________________________________________________
439 {
440  // sets the pools for event mizing in zvertex
441  if(fMultPoolLims) delete [] fMultPoolLims;
442  fNMultPools=nPools;
443  fNMultPoolsLimSize=nPools+1;
445  for(Int_t ib=0; ib<nPools+1; ib++) fMultPoolLims[ib]=multLimits[ib];
446  return;
447 }
448 //________________________________________________________________________
450 {
452  //
453  if(fDebug > 1) printf("AnalysisTaskCombinHF::UserCreateOutputObjects() \n");
454 
455  fOutput = new TList();
456  fOutput->SetOwner();
457  fOutput->SetName("OutputHistos");
458 
459  fHistNEvents = new TH1F("hNEvents", "number of events ",10,-0.5,9.5);
460  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
461  fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
462  fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
463  fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to phys sel");
464  fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected due to not reco vertex");
465  fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for contr vertex");
466  fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for zSPD-zTrack");
467  fHistNEvents->GetXaxis()->SetBinLabel(8,"n. rejected for vertex out of accept");
468  fHistNEvents->GetXaxis()->SetBinLabel(9,"n. rejected for pileup");
469  fHistNEvents->GetXaxis()->SetBinLabel(10,"no. of out centrality events");
470 
471  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
472  fHistNEvents->SetMinimum(0);
473  fOutput->Add(fHistNEvents);
474 
475  fHistEventMultCent = new TH2F("hEventMultCent"," ; Centrality (V0M) ; N_{tracklets} (|#eta|<1)",100,0.,100.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
477 
478  fHistEventMultCentEvSel = new TH2F("hEventMultCentEvSel"," ; Centrality (V0M) ; N_{tracklets} (|#eta|<1)",100,0.,100.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
480 
481  fHistEventMultZv = new TH2F("hEventMultZv"," ; z_{vertex} (cm) ; N_{tracklets} (|#eta|<1)",30,-15.,15.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
483 
484  fHistEventMultZvEvSel = new TH2F("hEventMultZvEvSel"," ; z_{vertex} (cm) ; N_{tracklets} (|#eta|<1)",30,-15.,15.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
486 
487  fHistXsecVsPtHard = new TH1F("hXsecVsPtHard", " ; pthard (GeV/c) ; Xsec", 200,0.,100.);
489 
490  fHistTrackStatus = new TH1F("hTrackStatus", "",8,-0.5,7.5);
491  fHistTrackStatus->GetXaxis()->SetBinLabel(1,"Not OK");
492  fHistTrackStatus->GetXaxis()->SetBinLabel(2,"Track OK");
493  fHistTrackStatus->GetXaxis()->SetBinLabel(3,"Kaon, Not OK");
494  fHistTrackStatus->GetXaxis()->SetBinLabel(4,"Kaon OK");
495  fHistTrackStatus->GetXaxis()->SetBinLabel(5,"Pion, Not OK");
496  fHistTrackStatus->GetXaxis()->SetBinLabel(6,"Pion OK");
497  fHistTrackStatus->GetXaxis()->SetBinLabel(7,"Kaon||Pion, Not OK");
498  fHistTrackStatus->GetXaxis()->SetBinLabel(8,"Kaon||Pion OK");
499  fHistTrackStatus->GetXaxis()->SetNdivisions(1,kFALSE);
500  fHistTrackStatus->SetMinimum(0);
502 
503  fHistTrackEtaMultZv = new TH3F("hTrackEtaMultZv","",40,-1.,1.,30,-15.,15.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
505  fHistSelTrackPhiPt = new TH2F("hSelTrackPhiPt"," ; #varphi ; p_{T} (GeV/c)",180,0.,2.*TMath::Pi(),20,0.,10.);
507 
508  Int_t nPtBins = (Int_t)(fMaxPt/fPtBinWidth+0.001);
510 
511  if(fReadMC){
512 
513  fHistCheckOrigin=new TH2F("hCheckOrigin"," ; origin ; generator",7,-1.5,5.5,2,-0.5,1.5);
514  fHistCheckOrigin->GetYaxis()->SetBinLabel(1,"Hijing");
515  fHistCheckOrigin->GetYaxis()->SetBinLabel(2,"Injected");
516  fHistCheckOriginRecoD=new TH2F("hCheckOriginRecoD"," ; origin ; generator",7,-1.5,5.5,2,-0.5,1.5);
517  fHistCheckOriginRecoD->GetYaxis()->SetBinLabel(1,"Hijing");
518  fHistCheckOriginRecoD->GetYaxis()->SetBinLabel(2,"Injected");
519  fHistCheckOriginRecoVsGen=new TH2F("hCheckOriginRecoVsGen"," ; Origin (reco D) ; Origin (gen D)",7,-1.5,5.5,7,-1.5,5.5);
523 
524  fHistCheckDecChan=new TH1F("hCheckDecChan","",7,-2.5,4.5);
525  fHistCheckDecChan->SetMinimum(0);
527 
528  fHistCheckDecChanAcc=new TH1F("hCheckDecChanAcc","",7,-2.5,4.5);
529  fHistCheckDecChanAcc->SetMinimum(0);
531 
532  fPtVsYVsMultGenPrompt = new TH3F("hPtVsYVsMultGenPrompt","",nPtBins,0.,maxPt,20,-1.,1.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
534 
535  fPtVsYVsMultGenLargeAccPrompt = new TH3F("hPtVsYVsMultGenLargeAccPrompt","",nPtBins,0.,maxPt,20,-1.,1.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
537 
538  fPtVsYVsMultGenLimAccPrompt = new TH3F("hPtVsYVsMultGenLimAccPrompt","",nPtBins,0.,maxPt,20,-1.,1.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
540 
541  fPtVsYVsMultGenAccPrompt = new TH3F("hPtVsYVsMultGenAccPrompt","",nPtBins,0.,maxPt,20,-1.,1.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
543 
544  fPtVsYVsMultGenAccEvSelPrompt = new TH3F("hPtVsYVsMultGenAccEvSelPrompt","",nPtBins,0.,maxPt,20,-1.,1.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
546 
547  fPtVsYVsMultRecoPrompt = new TH3F("hPtVsYVsMultRecoPrompt","",nPtBins,0.,maxPt,20,-1.,1.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
549 
550  fPtVsYVsMultGenFeeddw = new TH3F("hPtVsYVsMultGenFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
552 
553  fPtVsYVsMultGenLargeAccFeeddw = new TH3F("hPtVsYVsMultGenLargeAccFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
555 
556  fPtVsYVsMultGenLimAccFeeddw = new TH3F("hPtVsYVsMultGenLimAccFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
558 
559  fPtVsYVsMultGenAccFeeddw = new TH3F("hPtVsYVsMultGenAccFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
561 
562  fPtVsYVsMultGenAccEvSelFeeddw = new TH3F("hPtVsYVsMultGenAccEvSelFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
564 
565  fPtVsYVsMultRecoFeeddw = new TH3F("hPtVsYVsMultRecoFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,fNumOfMultBins,fMinMultiplicity,fMaxMultiplicity);
567 
568  fPtVsYVsPtBGenFeeddw = new TH3F("hPtVsYVsPtBGenFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,100,0.,50.);
570 
571  fPtVsYVsPtBGenLargeAccFeeddw = new TH3F("hPtVsYVsPtBGenLargeAccFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,100,0.,50.);
573 
574  fPtVsYVsPtBGenLimAccFeeddw = new TH3F("hPtVsYVsPtBGenLimAccFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,100,0.,50.);
576 
577  fPtVsYVsPtBGenAccFeeddw = new TH3F("hPtVsYVsPtBGenAccFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,100,0.,50.);
579 
580  fPtVsYVsPtBGenAccEvSelFeeddw = new TH3F("hPtVsYVsPtBGenAccEvSelFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,100,0.,50.);
582 
583  fPtVsYVsPtBRecoFeeddw = new TH3F("hPtVsYVsPtBRecoFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,100,0.,50.);
585  }
586 
587 
588  Int_t nMassBins=static_cast<Int_t>(fMaxMass*1000.-fMinMass*1000.);
589  Double_t maxm=fMinMass+nMassBins*0.001;
590  fMassVsPtVsY=new TH3F("hMassVsPtVsY","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
591  fOutput->Add(fMassVsPtVsY);
592 
593  fMassVsPtVsYRot=new TH3F("hMassVsPtVsYRot","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
594  fOutput->Add(fMassVsPtVsYRot);
595 
596  fMassVsPtVsYLSpp=new TH3F("hMassVsPtVsYLSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
598  fMassVsPtVsYLSmm=new TH3F("hMassVsPtVsYLSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
600 
601  fMassVsPtVsYSig=new TH3F("hMassVsPtVsYSig","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
602  fOutput->Add(fMassVsPtVsYSig);
603 
604  fMassVsPtVsYRefl=new TH3F("hMassVsPtVsYRefl","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
606 
607  fMassVsPtVsYBkg=new TH3F("hMassVsPtVsYBkg","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
608  fOutput->Add(fMassVsPtVsYBkg);
609 
610  fBMohterPtGen=new TH1F("hBMohterPtGen","",100,0.,50.);
611  fOutput->Add(fBMohterPtGen);
612 
613  fNSelected=new TH1F("hNSelected","",100,-0.5,99.5);
614  fOutput->Add(fNSelected);
615 
616  fNormRotated=new TH1F("hNormRotated","",11,-0.5,10.5);
617  fOutput->Add(fNormRotated);
618 
619  fDeltaMass=new TH1F("hDeltaMass","",100,-0.4,0.4);
620  fOutput->Add(fDeltaMass);
621 
622  Int_t binSparseDMassRot[5]={nMassBins,100,24,40,20};
623  Double_t edgeLowSparseDMassRot[5]={fMinMass,-0.4,0.,-4.,0};
624  Double_t edgeHighSparseDMassRot[5]={maxm,0.4,12.,4.,3.14};
625  fDeltaMassFullAnalysis=new THnSparseF("fDeltaMassFullAnalysis","fDeltaMassFullAnalysis;inv mass (GeV/c);#Delta inv mass (GeV/c) ; p_{T}^{D} (GeV/c); #Delta p_{T} (GeV/c); daughter angle (2prongs) (rad);",5,binSparseDMassRot,edgeLowSparseDMassRot,edgeHighSparseDMassRot);
627 
628  fMassVsPtVsYME=new TH3F("hMassVsPtVsYME","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
629  fOutput->Add(fMassVsPtVsYME);
630 
631  fMassVsPtVsYMELSpp=new TH3F("hMassVsPtVsYMELSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
633 
634  fMassVsPtVsYMELSmm=new TH3F("hMassVsPtVsYMELSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
636 
639  if(fDoEventMixing==2) fNOfPools=1;
641  fEventsPerPool=new TH2F("hEventsPerPool","hEventsPerPool",fNzVertPools,fzVertPoolLims,fNMultPools,fMultPoolLims);
642  fMixingsPerPool=new TH2F("hMixingsPerPool","hMixingsPerPool",fNzVertPools,fzVertPoolLims,fNMultPools,fMultPoolLims);
643  }else{
644  fEventsPerPool=new TH2F("hEventsPerPool","hEventsPerPool",1,-10.,10.,1,-0.5,2000.5);
645  fMixingsPerPool=new TH2F("hMixingsPerPool","hMixingsPerPool",1,-10.,10.,1,-0.5,2000.5);
646  }
647  fOutput->Add(fEventsPerPool);
648  fOutput->Add(fMixingsPerPool);
649 
650  fMassVsPtVsCosthSt=new TH3F("hMassVsPtVsCosthSt","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,6,0.4,1.);
651  fMassVsPtVsCosthStRot=new TH3F("hMassVsPtVsCosthStRot","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,6,0.4,1.);
652  fMassVsPtVsCosthStLSpp=new TH3F("hMassVsPtVsCosthStLSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,6,0.4,1.);
653  fMassVsPtVsCosthStLSmm=new TH3F("hMassVsPtVsCosthStLSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,6,0.4,1.);
654  fMassVsPtVsCosthStSig=new TH3F("hMassVsPtVsCosthStSig","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,6,0.4,1.);
655  fMassVsPtVsCosthStRefl=new TH3F("hMassVsPtVsCosthStRefl","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,6,0.4,1.);
656  fMassVsPtVsCosthStBkg=new TH3F("hMassVsPtVsCosthStBkg","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,6,0.4,1.);
657  fMassVsPtVsCosthStME=new TH3F("hMassVsPtVsCosthStME","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,6,0.4,1.);
658  fMassVsPtVsCosthStMELSpp=new TH3F("hMassVsPtVsCosthStMELSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,6,0.4,1.);
659  fMassVsPtVsCosthStMELSmm=new TH3F("hMassVsPtVsCosthStMELSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,6,0.4,1.);
670 
671  fHistonSigmaTPCPion=new TH2F("hnSigmaTPCPion"," ; p (GeV/c) ; n#sigma^{#pi}_{TPC}",20,0.,10.,100,-5.,5.);
672  fHistonSigmaTPCPionGoodTOF=new TH2F("hnSigmaTPCPionGoodTOF"," ; p (GeV/c) ; n#sigma^{#pi}_{TPC}",20,0.,10.,100,-5.,5.);
673  fHistonSigmaTOFPion=new TH2F("hnSigmaTOFPion"," ; p (GeV/c) ; n#sigma^{#pi}_{TOF}",20,0.,10.,100,-5.,5.);
674  fHistonSigmaTPCKaon=new TH2F("hnSigmaTPCKaon"," ; p (GeV/c) ; n#sigma^{K}_{TPC}",20,0.,10.,100,-5.,5.);
675  fHistonSigmaTPCKaonGoodTOF=new TH2F("hnSigmaTPCKaonGoodTOF"," ; p (GeV/c) ; n#sigma^{K}_{TPC}",20,0.,10.,100,-5.,5.);
676  fHistonSigmaTOFKaon=new TH2F("hnSigmaTOFKaon"," ; p (GeV/c) ; n#sigma^{K}_{TOF}",20,0.,10.,100,-5.,5.);
683 
684  fHistoPtKPtPiPtD = new TH3F("hPtKPtPiPtD"," ; p_{T}(D) ; p_{T}(K) ; p_{T}(#pi)",32,0.,16.,60,0.,15.,60,0.,15.);
685  fHistoPtKPtPiPtDSig = new TH3F("hPtKPtPiPtDSig"," ; p_{T}(D) ; p_{T}(K) ; p_{T}(#pi)",32,0.,16.,60,0.,15.,60,0.,15.);
688 
689  //Counter for Normalization
690  fCounter = new AliNormalizationCounter("NormalizationCounter");
691  fCounter->Init();
692 
693  fListCuts = new TList();
694  fListCuts->SetOwner();
695  if(fTrackCutsAll){
696  AliESDtrackCuts* tatosave=new AliESDtrackCuts(*fTrackCutsAll);
697  fListCuts->Add(tatosave);
698  }
699  if(fTrackCutsPion){
700  AliESDtrackCuts* tptosave=new AliESDtrackCuts(*fTrackCutsPion);
701  tptosave->SetName(Form("%sForPions",fTrackCutsPion->GetName()));
702  fListCuts->Add(tptosave);
703  }
704  if(fTrackCutsKaon){
705  AliESDtrackCuts* tktosave=new AliESDtrackCuts(*fTrackCutsKaon);
706  tktosave->SetName(Form("%sForKaons",fTrackCutsKaon->GetName()));
707  fListCuts->Add(tktosave);
708  }
709 
710  if(fAnalysisCuts->GetPidHF()){
711  AliAODPidHF* pidtosave=new AliAODPidHF(*(fAnalysisCuts->GetPidHF()));
712  fListCuts->Add(pidtosave);
713  }
714  TH1F* hCutValues = new TH1F("hCutValues","",8,0.5,8.5);
715  hCutValues->SetBinContent(1,fFilterMask);
716  hCutValues->GetXaxis()->SetBinLabel(1,"Filter bit");
717  hCutValues->SetBinContent(2,fCutTPCSignalN);
718  hCutValues->GetXaxis()->SetBinLabel(2,"n TPC clu for PID");
719  hCutValues->SetBinContent(3,(Float_t)fApplyCutCosThetaStar);
720  hCutValues->GetXaxis()->SetBinLabel(3,"Use costhetastar (D0)");
721  hCutValues->SetBinContent(4,fCutCosThetaStar);
722  hCutValues->GetXaxis()->SetBinLabel(4,"costhetastar (D0)");
723  hCutValues->SetBinContent(5,fPhiMassCut);
724  hCutValues->GetXaxis()->SetBinLabel(5,"phi mass (Ds)");
725  hCutValues->SetBinContent(6,fCutCos3PiKPhiRFrame);
726  hCutValues->GetXaxis()->SetBinLabel(6,"cos3piK (Ds)");
727  hCutValues->SetBinContent(7,fCutCosPiDsLabFrame);
728  hCutValues->GetXaxis()->SetBinLabel(7,"cospiDs (Ds)");
729  hCutValues->SetBinContent(8,fAnalysisCuts->GetUseTimeRangeCutForPbPb2018());
730  hCutValues->GetXaxis()->SetBinLabel(8,"TimeRangeCut");
731  fListCuts->Add(hCutValues);
732  PostData(3, fListCuts);
733 
734 
735  fKaonTracks = new TObjArray();
736  fPionTracks=new TObjArray();
737  fKaonTracks->SetOwner();
738  fPionTracks->SetOwner();
739 
740  fEventBuffer = new TTree*[fNOfPools];
741  for(Int_t i=0; i<fNOfPools; i++){
742  fEventBuffer[i]=new TTree(Form("EventBuffer_%d",i), "Temporary buffer for event mixing");
743  fEventBuffer[i]->Branch("zVertex", &fVtxZ);
744  fEventBuffer[i]->Branch("multiplicity", &fMultiplicity);
745  fEventBuffer[i]->Branch("eventInfo", "TObjString",&fEventInfo);
746  fEventBuffer[i]->Branch("karray", "TObjArray", &fKaonTracks);
747  fEventBuffer[i]->Branch("parray", "TObjArray", &fPionTracks);
748  }
749 
750  PostData(1,fOutput);
751  PostData(2,fCounter);
752 }
753 
754 //________________________________________________________________________
757 
758  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
759  if(!aod && AODEvent() && IsStandardAOD()) {
760  // In case there is an AOD handler writing a standard AOD, use the AOD
761  // event in memory rather than the input (ESD) event.
762  aod = dynamic_cast<AliAODEvent*> (AODEvent());
763  }
764  if(!aod){
765  printf("AliAnalysisTaskCombinHF::UserExec: AOD not found!\n");
766  return;
767  }
768 
769  // fix for temporary bug in ESDfilter
770  // the AODs with null vertex pointer didn't pass the PhysSel
771  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
772 
773  if(fReadMC){
774  // Reject events with trigger mask 0 of the LHC13d3 production
775  // For these events the ITS layers are skipped in the trakcing
776  // and the vertex reconstruction efficiency from tracks is biased
777  Int_t runnumber = aod->GetRunNumber();
778  if(aod->GetTriggerMask()==0 &&
779  (runnumber>=195344 && runnumber<=195677)){
780  return;
781  }
782  // Set the trigger mask for physics selection to kMB in the MC
784  // printf("Enforce trigger mask to kMB, previous mask = %d\n",fAnalysisCuts->GetTriggerMask());
785  fAnalysisCuts->SetTriggerMask(AliVEvent::kMB);
786  }
787  for(Int_t j=0; j<200000; j++) fOrigContainer[j]=-1;
788  }
789 
790  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
791  AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
792 
793 
794  fHistNEvents->Fill(0); // count event
795  // Post the data already here
796  PostData(1,fOutput);
797 
799 
800  Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
804  }else{
806  fHistNEvents->Fill(9);
807  }else{
812  }else{
815  }
816  }
817  }
818 
819  // PID object should be taken AFTER IsEventSelected to have proper call to SetupPid !!
820  AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
822  fPidHF->SetPidResponse(pidResp);
823  //
824 
825  Int_t ntracks=aod->GetNumberOfTracks();
826  fVtxZ = aod->GetPrimaryVertex()->GetZ();
828  Float_t evCentr=fAnalysisCuts->GetCentrality(aod);
831  fHistEventMultCent->Fill(evCentr,fMultiplicity);
832  }
833 
835  // events not passing the centrality selection can be removed immediately. For the others we must count the generated D mesons
836 
837 
838  TClonesArray *arrayMC=0;
839  AliAODMCHeader *mcHeader=0;
840  if(fReadMC){
841  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
842  if(!arrayMC) {
843  printf("AliAnalysisTaskCombinHF::UserExec: MC particles branch not found!\n");
844  return;
845  }
846 
847  // load MC header
848  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
849  if(!mcHeader) {
850  printf("AliAnalysisTaskCombinHF::UserExec: MC header branch not found!\n");
851  return;
852  }
853  // selection on pt hard bins in Pb-Pb
854  if(fSelectPtHardRange){
855  TList *lh=mcHeader->GetCocktailHeaders();
856  if(lh){
857  Int_t nh=lh->GetEntries();
858  for(Int_t i=0;i<nh;i++){
859  AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
860  TString genname=gh->GetName();
861  if(genname.Contains("ythia") || genname.Contains("YTHIA")){
862  AliGenPythiaEventHeader* pyth=(AliGenPythiaEventHeader*)lh->At(i);
863  Double_t ptha=pyth->GetPtHard();
864  Double_t xsec=pyth->GetXsection();
865  if(ptha<fMinPtHard || ptha>fMaxPtHard) return;
866  fHistXsecVsPtHard->SetBinContent(fHistXsecVsPtHard->GetXaxis()->FindBin(ptha),xsec);
867  }
868  }
869  }
870  }
871  Double_t zMCVertex = mcHeader->GetVtxZ();
872  if (TMath::Abs(zMCVertex) < fAnalysisCuts->GetMaxVtxZ()){ // only cut on zVertex applied to count the signal
873  FillGenHistos(arrayMC,mcHeader,isEvSel);
874  }
875  fHistEventMultZv->Fill(zMCVertex,fMultiplicity);
876  if(isEvSel) fHistEventMultZvEvSel->Fill(zMCVertex,fMultiplicity);
877  // switch off event mixing in case of signal only MC
879  }else{
881  if(isEvSel) fHistEventMultZvEvSel->Fill(fVtxZ,fMultiplicity);
882  }
883 
884 
885  if(!isEvSel)return;
886 
887  fHistNEvents->Fill(1);
888  fHistEventMultCentEvSel->Fill(evCentr,fMultiplicity);
889 
890 
891  Int_t pdgOfD=421;
892  if(fMeson==kDplus) pdgOfD=411;
893  else if(fMeson==kDs) pdgOfD=431;
894 
895  // select and flag tracks
896  UChar_t* status = new UChar_t[ntracks];
897  for(Int_t iTr=0; iTr<ntracks; iTr++){
898  status[iTr]=0;
899  AliAODTrack* track=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr));
900  if(!track){
901  AliWarning("Error in casting track to AOD track. Not a standard AOD?");
902  continue;
903  }
904  if(fReadMC && fSignalOnlyMC && arrayMC){
905  // for fast MC analysis we skip tracks not coming from charm hadrons
906  Bool_t isCharm=AliVertexingHFUtils::IsTrackFromHadronDecay(pdgOfD,track,arrayMC);
907  if(!isCharm) continue;
908  }
909  if(IsTrackSelected(track)) status[iTr]+=1;
910 
911  // PID
912  if (fPIDstrategy == knSigma) {
913  // nsigma PID
914  Double_t trmom=track->P();
915  Bool_t okTOF=fPidHF->CheckTOFPIDStatus(track);
916  if(IsKaon(track)){
917  Double_t nstpc,nstof;
918  fPidHF->GetnSigmaTPC(track,AliPID::kKaon,nstpc);
919  fPidHF->GetnSigmaTOF(track,AliPID::kKaon,nstof);
920  fHistonSigmaTPCKaon->Fill(trmom,nstpc);
921  if(okTOF) fHistonSigmaTPCKaonGoodTOF->Fill(trmom,nstpc);
922  fHistonSigmaTOFKaon->Fill(trmom,nstof);
923  status[iTr]+=2;
924  }
925  if(IsPion(track)){
926  Double_t nstpc,nstof;
927  fPidHF->GetnSigmaTPC(track,AliPID::kPion,nstpc);
928  fPidHF->GetnSigmaTOF(track,AliPID::kPion,nstof);
929  fHistonSigmaTPCPion->Fill(trmom,nstpc);
930  if(okTOF) fHistonSigmaTPCPionGoodTOF->Fill(trmom,nstpc);
931  fHistonSigmaTOFPion->Fill(trmom,nstof);
932  status[iTr]+=4;
933  }
934  }
936  // Bayesian PID
937  Double_t *weights = new Double_t[AliPID::kSPECIES];
938  fPidHF->GetPidCombined()->ComputeProbabilities(track, fPidHF->GetPidResponse(), weights);
940  if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kKaon]) status[iTr] += 2;
941  if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kPion]) status[iTr] += 4;
942  }
943  if (fPIDstrategy == kBayesianThres) {
944  if (weights[AliPID::kKaon] > fBayesThresKaon) status[iTr] += 2;
945  if (weights[AliPID::kPion] > fBayesThresPion) status[iTr] += 4;
946  }
947  delete[] weights;
948  }
949 
950  fHistTrackStatus->Fill(status[iTr]);
951  fHistTrackEtaMultZv->Fill(track->Eta(),fVtxZ,fMultiplicity);
952  if(status[iTr]>0) fHistSelTrackPhiPt->Fill(track->Phi(),track->Pt());
953  }
954 
955  // build the combinatorics
956  Int_t nSelected=0;
957  Int_t nFiltered=0;
958  Double_t dummypos[3]={0.,0.,0.};
959  AliAODVertex* v2=new AliAODVertex(dummypos,999.,-1,2);
960  AliAODVertex* v3=new AliAODVertex(dummypos,999.,-1,3);
961  // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
962  Double_t d02[2]={0.,0.};
963  Double_t d03[3]={0.,0.,0.};
964  AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
965  AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
966  UInt_t pdg0[2]={321,211};
967  UInt_t pdgp[3]={321,211,211};
968  UInt_t pdgs[3]={321,211,321};
969  Double_t tmpp[3];
970  Double_t px[3],py[3],pz[3];
971  Int_t dgLabels[3];
972  fKaonTracks->Delete();
973  fPionTracks->Delete();
974 
975  for(Int_t iTr1=0; iTr1<ntracks; iTr1++){
976  AliAODTrack* trK=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr1));
977  if(!trK){
978  AliWarning("Error in casting track to AOD track. Not a standard AOD?");
979  continue;
980  }
981  if((status[iTr1] & 1)==0) continue;
982  if(fDoEventMixing>0){
983  if(status[iTr1] & 2) fKaonTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
984  if(status[iTr1] & 4) fPionTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
985  }
986  if((status[iTr1] & 2)==0) continue;
987  Int_t chargeK=trK->Charge();
988  trK->GetPxPyPz(tmpp);
989  px[0] = tmpp[0];
990  py[0] = tmpp[1];
991  pz[0] = tmpp[2];
992  dgLabels[0]=trK->GetLabel();
993  for(Int_t iTr2=0; iTr2<ntracks; iTr2++){
994  if((status[iTr2] & 1)==0) continue;
995  if((status[iTr2] & 4)==0) continue;
996  if(iTr1==iTr2) continue;
997  AliAODTrack* trPi1=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr2));
998  if(!trPi1){
999  AliWarning("Error in casting track to AOD track. Not a standard AOD?");
1000  continue;
1001  }
1002  Int_t chargePi1=trPi1->Charge();
1003  trPi1->GetPxPyPz(tmpp);
1004  px[1] = tmpp[0];
1005  py[1] = tmpp[1];
1006  pz[1] = tmpp[2];
1007  dgLabels[1]=trPi1->GetLabel();
1008  if(fMeson==kDzero){
1009  if(chargePi1==chargeK){
1010  // LS candidate
1011  FillLSHistos(421,2,tmpRD2,px,py,pz,pdg0,chargePi1);
1012  }else{
1013  // OS candidate
1014  nFiltered++;
1015  v2->AddDaughter(trK);
1016  v2->AddDaughter(trPi1);
1017  tmpRD2->SetSecondaryVtx(v2);
1018  Bool_t ok=FillHistos(421,2,tmpRD2,px,py,pz,pdg0,arrayMC,mcHeader,dgLabels);
1019  v2->RemoveDaughters();
1020  if(ok) nSelected++;
1021  }
1022  }else{
1023  for(Int_t iTr3=iTr2+1; iTr3<ntracks; iTr3++){
1024  if((status[iTr3] & 1)==0) continue;
1025  if(fMeson==kDplus && (status[iTr3] & 4)==0) continue;
1026  if(fMeson==kDs && (status[iTr3] & 2)==0) continue;
1027  if(iTr1==iTr3) continue;
1028  AliAODTrack* trPi2=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr3));
1029  if(!trPi2){
1030  AliWarning("Error in casting track to AOD track. Not a standard AOD?");
1031  continue;
1032  }
1033  Int_t chargePi2=trPi2->Charge();
1034  trPi2->GetPxPyPz(tmpp);
1035  px[2] = tmpp[0];
1036  py[2] = tmpp[1];
1037  pz[2] = tmpp[2];
1038  dgLabels[2]=trPi2->GetLabel();
1039  if(fMeson==kDs){
1040  Double_t massKK=ComputeInvMassKK(trK,trPi2);
1041  Double_t deltaMass=massKK-TDatabasePDG::Instance()->GetParticle(333)->Mass();
1042  if(TMath::Abs(deltaMass)>fPhiMassCut) continue;
1043  tmpRD3->SetPxPyPzProngs(3,px,py,pz);
1044  Double_t cos1=((AliAODRecoDecayHF3Prong*)tmpRD3)->CosPiKPhiRFrameKpiK();
1045  Double_t kincutPiKPhi=TMath::Abs(cos1*cos1*cos1);
1046  if(kincutPiKPhi<fCutCos3PiKPhiRFrame) continue;
1047  Double_t cosPiDsLabFrame=((AliAODRecoDecayHF3Prong*)tmpRD3)->CosPiDsLabFrameKpiK();
1048  if(cosPiDsLabFrame>fCutCosPiDsLabFrame) continue;
1049  }
1050  Bool_t isThreeLS=kFALSE;
1051  if(chargePi1==chargeK && chargePi2==chargeK){
1052  isThreeLS=kTRUE;
1053  if(fMeson==kDplus)FillLSHistos(411,3,tmpRD3,px,py,pz,pdgp,chargePi1);
1054  else if(fMeson==kDs)FillLSHistos(431,3,tmpRD3,px,py,pz,pdgs,chargePi1);
1055  }
1056  Bool_t acceptOS=kFALSE;
1057  if(fMeson==kDplus){
1058  if(chargePi1!=chargeK && chargePi2!=chargeK)acceptOS=kTRUE;
1059  }else if(fMeson==kDs){
1060  if(chargePi2!=chargeK && !isThreeLS) acceptOS=kTRUE;
1061  }
1062  if(acceptOS){
1063  nFiltered++;
1064  v3->AddDaughter(trK);
1065  v3->AddDaughter(trPi1);
1066  v3->AddDaughter(trPi2);
1067  tmpRD3->SetSecondaryVtx(v3);
1068  Bool_t ok=kFALSE;
1069  if(fMeson==kDplus) ok=FillHistos(411,3,tmpRD3,px,py,pz,pdgp,arrayMC,mcHeader,dgLabels);
1070  else if(fMeson==kDs) ok=FillHistos(431,3,tmpRD3,px,py,pz,pdgs,arrayMC,mcHeader,dgLabels);
1071  v3->RemoveDaughters();
1072  if(ok) nSelected++;
1073  }
1074  }
1075  }
1076  }
1077  }
1078 
1079  delete [] status;
1080  delete v2;
1081  delete v3;
1082  delete tmpRD2;
1083  delete tmpRD3;
1084 
1085  fNSelected->Fill(nSelected);
1086 
1087  fCounter->StoreCandidates(aod,nFiltered,kTRUE);
1088  fCounter->StoreCandidates(aod,nSelected,kFALSE);
1089  fEventInfo->SetString(Form("Ev%d_esd%d_Pi%d_K%d",mgr->GetNcalls(),((AliAODHeader*)aod->GetHeader())->GetEventNumberESDFile(),fPionTracks->GetEntries(),fKaonTracks->GetEntries()));
1090  if(fDoEventMixing==1){
1092  if(ind>=0 && ind<fNOfPools){
1094  fEventBuffer[ind]->Fill();
1095  if(fEventBuffer[ind]->GetEntries() >= fNumberOfEventsForMixing){
1097  DoMixingWithPools(ind);
1098  ResetPool(ind);
1099  }
1100  }
1101  }else if(fDoEventMixing==2){ // mix with cuts, no pools
1102  fEventBuffer[0]->Fill();
1103  }
1104  PostData(1,fOutput);
1105  PostData(2,fCounter);
1106 
1107  return;
1108 }
1109 
1110 //________________________________________________________________________
1111 void AliAnalysisTaskCombinHF::FillLSHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, Int_t charge){
1113 
1114  if(fReadMC && fSignalOnlyMC) return;
1115  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
1116  Double_t pt = tmpRD->Pt();
1117  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
1118  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
1119  Double_t rapid = tmpRD->Y(pdgD);
1120  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
1121  Bool_t fillLS=kTRUE;
1122  Double_t costhst=0;
1123  Double_t absCosThSt=0;
1124  if(TMath::Abs(pdgD)==421 && (fApplyCutCosThetaStar || fFillHistosVsCosThetaStar)){
1125  costhst=tmpRD->CosThetaStar(0,421,321,211); // kaon is the first daughter
1126  absCosThSt=TMath::Abs(costhst);
1127  if(fApplyCutCosThetaStar && absCosThSt>fCutCosThetaStar) fillLS=kFALSE;
1128  }
1129  if(fillLS){
1130  Double_t invMass=TMath::Sqrt(minv2);
1131  if(charge>0) fMassVsPtVsYLSpp->Fill(invMass,pt,rapid);
1132  else fMassVsPtVsYLSmm->Fill(invMass,pt,rapid);
1134  if(charge>0) fMassVsPtVsCosthStLSpp->Fill(invMass,pt,absCosThSt);
1135  else fMassVsPtVsCosthStLSmm->Fill(invMass,pt,absCosThSt);
1136  }
1137  }
1138  }
1139  }
1140  return;
1141 }
1142 
1143 //________________________________________________________________________
1144 void AliAnalysisTaskCombinHF::FillGenHistos(TClonesArray* arrayMC, AliAODMCHeader *mcHeader, Bool_t isEvSel){
1146  Int_t totPart=arrayMC->GetEntriesFast();
1147  Int_t thePDG=411;
1148  Int_t nProng=3;
1149  if(fMeson==kDzero){
1150  thePDG=421;
1151  nProng=2;
1152  }else if(fMeson==kDs){
1153  thePDG=431;
1154  nProng=3;
1155  }
1156  for(Int_t ip=0; ip<totPart; ip++){
1157  AliAODMCParticle *part = (AliAODMCParticle*)arrayMC->At(ip);
1158  if(TMath::Abs(part->GetPdgCode())==thePDG){
1160  if(ip<200000) fOrigContainer[ip]=orig;
1161  Bool_t isInj=AliVertexingHFUtils::IsTrackInjected(ip,mcHeader,arrayMC);
1162  Int_t deca=0;
1163  Bool_t isGoodDecay=kFALSE;
1164  Int_t labDau[4]={-1,-1,-1,-1};
1165  if(fMeson==kDzero){
1166  deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,part,labDau);
1167  if(part->GetNDaughters()!=2) continue;
1168  if(deca==1) isGoodDecay=kTRUE;
1169  }else if(fMeson==kDplus){
1170  deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,part,labDau);
1171  if(deca>0) isGoodDecay=kTRUE;
1172  }else if(fMeson==kDs){
1173  deca=AliVertexingHFUtils::CheckDsDecay(arrayMC,part,labDau);
1174  if(deca==1) isGoodDecay=kTRUE;
1175  }
1176  if(labDau[0]==-1){
1177  // printf(Form("Meson %d Label of daughters not filled correctly -- %d\n",fMeson,isGoodDecay));
1178  continue; //protection against unfilled array of labels
1179  }
1180  fHistCheckDecChan->Fill(deca);
1181  Bool_t isInAcc=CheckAcceptance(arrayMC,nProng,labDau);
1182  if(isInAcc) fHistCheckDecChanAcc->Fill(deca);
1183  if(isGoodDecay){
1184  Double_t ptgen=part->Pt();
1185  Double_t ygen=part->Y();
1186  Double_t ptbmoth=0.;
1187  if(orig==5) ptbmoth=AliVertexingHFUtils::GetBeautyMotherPt(arrayMC,part);
1188  if(fAnalysisCuts->IsInFiducialAcceptance(ptgen,ygen)){
1189  fHistCheckOrigin->Fill(orig,isInj);
1190  if(orig==4){
1191  fPtVsYVsMultGenPrompt->Fill(ptgen,ygen,fMultiplicity);
1192  if(TMath::Abs(ygen)<0.5) fPtVsYVsMultGenLimAccPrompt->Fill(ptgen,ygen,fMultiplicity);
1193  if(isInAcc) fPtVsYVsMultGenAccPrompt->Fill(ptgen,ygen,fMultiplicity);
1194  if(isEvSel && isInAcc) fPtVsYVsMultGenAccEvSelPrompt->Fill(ptgen,ygen,fMultiplicity);
1195  }else if(orig==5){
1196  fPtVsYVsMultGenFeeddw->Fill(ptgen,ygen,fMultiplicity);
1197  fPtVsYVsPtBGenFeeddw->Fill(ptgen,ygen,ptbmoth);
1198  if(TMath::Abs(ygen)<0.5){
1199  fPtVsYVsMultGenLimAccFeeddw->Fill(ptgen,ygen,fMultiplicity);
1200  fPtVsYVsPtBGenLimAccFeeddw->Fill(ptgen,ygen,ptbmoth);
1201  fBMohterPtGen->Fill(ptbmoth);
1202  }
1203  if(isInAcc){
1204  fPtVsYVsMultGenAccFeeddw->Fill(ptgen,ygen,fMultiplicity);
1205  fPtVsYVsPtBGenAccFeeddw->Fill(ptgen,ygen,ptbmoth);
1206  }
1207  if(isEvSel && isInAcc){
1208  fPtVsYVsMultGenAccEvSelFeeddw->Fill(ptgen,ygen,fMultiplicity);
1209  fPtVsYVsMultGenAccEvSelFeeddw->Fill(ptgen,ygen,ptbmoth);
1210  }
1211  }
1212  }
1213  if(TMath::Abs(ygen)<0.9){
1214  if(orig==4) fPtVsYVsMultGenLargeAccPrompt->Fill(ptgen,ygen,fMultiplicity);
1215  else if(orig==5){
1216  fPtVsYVsMultGenLargeAccFeeddw->Fill(ptgen,ygen,fMultiplicity);
1217  fPtVsYVsPtBGenLargeAccFeeddw->Fill(ptgen,ygen,ptbmoth);
1218  }
1219  }
1220  }
1221  }
1222  }
1223 }
1224 
1225 //________________________________________________________________________
1226 Bool_t AliAnalysisTaskCombinHF::FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, TClonesArray *arrayMC, AliAODMCHeader *mcHeader, Int_t* dgLabels){
1228 
1229  Bool_t accept=kFALSE;
1230 
1231  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
1232  Double_t pt = tmpRD->Pt();
1233  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
1234  Double_t mass=TMath::Sqrt(minv2);
1235 
1236  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
1237  Double_t rapid = tmpRD->Y(pdgD);
1238  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
1239  accept=kTRUE;
1240  Double_t costhst=0;
1241  Double_t absCosThSt=0;
1242  Double_t ptK=0;
1243  Double_t ptPi=0;
1244  if(TMath::Abs(pdgD)==421 && (fApplyCutCosThetaStar || fFillHistosVsCosThetaStar)){
1245  costhst=tmpRD->CosThetaStar(0,421,321,211); // kaon is the first daughter
1246  absCosThSt=TMath::Abs(costhst);
1247  if(fApplyCutCosThetaStar && absCosThSt>fCutCosThetaStar) accept=kFALSE;
1248  }
1249  if(accept){
1250  fMassVsPtVsY->Fill(mass,pt,rapid);
1251  if(fFillHistosVsCosThetaStar) fMassVsPtVsCosthSt->Fill(mass,pt,absCosThSt);
1252  if(pdgD==421){
1253  ptK=TMath::Sqrt(px[0]*px[0]+py[0]*py[0]);
1254  ptPi=TMath::Sqrt(px[1]*px[1]+py[1]*py[1]);
1255  if(TMath::Abs(mass-1.865)<0.025) fHistoPtKPtPiPtD->Fill(pt,ptK,ptPi);
1256  }
1257  if(fReadMC){
1258  Int_t signPdg[3]={0,0,0};
1259  for(Int_t iii=0; iii<nProngs; iii++) signPdg[iii]=pdgdau[iii];
1260  Int_t labD = tmpRD->MatchToMC(pdgD,arrayMC,nProngs,signPdg);
1261  if(labD>=0){
1262  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(TMath::Abs(dgLabels[0])));
1263  if(part){
1264  Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
1265  if(pdgCode==321){ // if the first daughter is a Kaon, this is signal with correct mass assignment
1266  fMassVsPtVsYSig->Fill(mass,pt,rapid);
1267  if(fFillHistosVsCosThetaStar) fMassVsPtVsCosthStSig->Fill(mass,pt,absCosThSt);
1268  if(pdgD==421) fHistoPtKPtPiPtDSig->Fill(pt,ptK,ptPi);
1269  AliAODMCParticle* dmes = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labD));
1270  if(dmes){
1272  Bool_t isInj=AliVertexingHFUtils::IsTrackInjected(labD,mcHeader,arrayMC);
1273  fHistCheckOriginRecoD->Fill(orig,isInj);
1274  if(labD<200000) fHistCheckOriginRecoVsGen->Fill(fOrigContainer[labD],orig);
1275  if(orig==4) fPtVsYVsMultRecoPrompt->Fill(dmes->Pt(),dmes->Y(),fMultiplicity);
1276  else if(orig==5){
1277  Double_t ptbmoth=AliVertexingHFUtils::GetBeautyMotherPt(arrayMC,dmes);
1278  fPtVsYVsMultRecoFeeddw->Fill(dmes->Pt(),dmes->Y(),fMultiplicity);
1279  fPtVsYVsPtBRecoFeeddw->Fill(dmes->Pt(),dmes->Y(),ptbmoth);
1280  }
1281  }
1282  }else{ // if the first daughter is not a kaon, it is a reflection
1283  fMassVsPtVsYRefl->Fill(mass,pt,rapid);
1284  if(fFillHistosVsCosThetaStar) fMassVsPtVsCosthStRefl->Fill(mass,pt,absCosThSt);
1285  }
1286  }
1287  }else{
1288  if(fSignalOnlyMC) accept=kFALSE;
1289  else fMassVsPtVsYBkg->Fill(mass,pt,rapid);
1290  if(fFillHistosVsCosThetaStar) fMassVsPtVsCosthStBkg->Fill(mass,pt,absCosThSt);
1291  }
1292  }
1293  }
1294  }
1295  }
1296  // skip track rotations in case of signal only MC
1297  if(fReadMC && fSignalOnlyMC) return accept;
1298 
1299  // Track rotations to estimate the background
1300  Int_t nRotated=0;
1301  Double_t massRot=0;// calculated later only if candidate is acceptable
1302  Double_t angleProngXY;
1303  if(TMath::Abs(pdgD)==421)angleProngXY=TMath::ACos((px[0]*px[1]+py[0]*py[1])/TMath::Sqrt((px[0]*px[0]+py[0]*py[0])*(px[1]*px[1]+py[1]*py[1])));
1304  else if(TMath::Abs(pdgD)==431) {
1305  Double_t px_phi = px[0]+px[2];
1306  Double_t py_phi = py[0]+py[2];
1307  Double_t pz_phi = pz[0]+pz[2];
1308  angleProngXY=TMath::ACos((pz_phi*px[1]+py_phi*py[1])/TMath::Sqrt((px_phi*px_phi+py_phi*py_phi)*(px[1]*px[1]+py[1]*py[1])));
1309  }
1310  else {//angle between pion and phi meson
1311  angleProngXY=TMath::ACos(((px[0]+px[1])*px[2]+(py[0]+py[1])*py[2])/TMath::Sqrt(((px[0]+px[1])*(px[0]+px[1])+(py[0]+py[1])*(py[0]+py[1]))*(px[2]*px[2]+py[2]*py[2])));
1312  }
1313  Double_t ptOrig=pt;
1314 
1315 
1316  Double_t rotStep=0.;
1317  if(fNRotations>1) rotStep=(fMaxAngleForRot-fMinAngleForRot)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
1318  if(TMath::Abs(pdgD)==421 || TMath::Abs(pdgD)==431) fNRotations3=1;
1319  Double_t rotStep3=0.;
1320  if(fNRotations3>1) rotStep3=(fMaxAngleForRot3-fMinAngleForRot3)/(fNRotations3-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
1321 
1322  for(Int_t irot=0; irot<fNRotations; irot++){
1323  Double_t phirot=fMinAngleForRot+rotStep*irot;
1324  Double_t tmpx=px[0];
1325  Double_t tmpy=py[0];
1326  Double_t tmpx2=px[2];
1327  Double_t tmpy2=py[2];
1328  if(pdgD==431) {
1329  //rotate pion w.r.t. phi meson
1330  tmpx=px[1];
1331  tmpy=py[1];
1332  px[1]=tmpx*TMath::Cos(phirot)-tmpy*TMath::Sin(phirot);
1333  py[1]=tmpx*TMath::Sin(phirot)+tmpy*TMath::Cos(phirot);
1334  }
1335  else {
1336  px[0]=tmpx*TMath::Cos(phirot)-tmpy*TMath::Sin(phirot);
1337  py[0]=tmpx*TMath::Sin(phirot)+tmpy*TMath::Cos(phirot);
1338  }
1339  for(Int_t irot3=0; irot3<fNRotations3; irot3++){
1340  if(pdgD==411){
1341  Double_t phirot2=fMaxAngleForRot3-rotStep3*irot;
1342  px[2]=tmpx*TMath::Cos(phirot2)-tmpy*TMath::Sin(phirot2);
1343  py[2]=tmpx*TMath::Sin(phirot2)+tmpy*TMath::Cos(phirot2);
1344  }
1345  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
1346  pt = tmpRD->Pt();
1347  minv2 = tmpRD->InvMass2(nProngs,pdgdau);
1348  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
1349  Double_t rapid = tmpRD->Y(pdgD);
1350  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
1351  Bool_t fillRotCase=kTRUE;
1352  Double_t costhst=0;
1353  Double_t absCosThSt=0;
1354  if(TMath::Abs(pdgD)==421 && (fApplyCutCosThetaStar || fFillHistosVsCosThetaStar)){
1355  costhst=tmpRD->CosThetaStar(0,421,321,211); // kaon is the first daughter
1356  absCosThSt=TMath::Abs(costhst);
1357  if(fApplyCutCosThetaStar && absCosThSt>fCutCosThetaStar) fillRotCase=kFALSE;
1358  }
1359  if(fillRotCase){
1360  massRot=TMath::Sqrt(minv2);
1361  fMassVsPtVsYRot->Fill(massRot,pt,rapid);
1362  if(fFillHistosVsCosThetaStar) fMassVsPtVsCosthStRot->Fill(massRot,pt,absCosThSt);
1363  nRotated++;
1364  fDeltaMass->Fill(massRot-mass);
1365  if(fFullAnalysis){
1366  Double_t pointRot[5]={mass,massRot-mass,ptOrig,pt-ptOrig,angleProngXY};
1367  fDeltaMassFullAnalysis->Fill(pointRot);
1368  }
1369  }
1370  }
1371  }
1372  }
1373  if(pdgD==431) {
1374  px[1]=tmpx;
1375  py[1]=tmpy;
1376  }
1377  else {
1378  px[0]=tmpx;
1379  py[0]=tmpy;
1380  if(pdgD==411){
1381  px[2]=tmpx2;
1382  py[2]=tmpy2;
1383  }
1384  }
1385  }
1386  fNormRotated->Fill(nRotated);
1387 
1388  return accept;
1389 
1390 }
1391 //________________________________________________________________________
1392 void AliAnalysisTaskCombinHF::FillMEHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau){
1394 
1395  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
1396  Double_t pt = tmpRD->Pt();
1397  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
1398  Double_t mass=TMath::Sqrt(minv2);
1399 
1400  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
1401  Double_t rapid = tmpRD->Y(pdgD);
1402  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
1403  Bool_t fillME=kTRUE;
1404  Double_t costhst=0;
1405  Double_t absCosThSt=0;
1406  if(TMath::Abs(pdgD)==421 && (fApplyCutCosThetaStar || fFillHistosVsCosThetaStar)){
1407  costhst=tmpRD->CosThetaStar(0,421,321,211); // kaon is the first daughter
1408  absCosThSt=TMath::Abs(costhst);
1409  if(fApplyCutCosThetaStar && absCosThSt>fCutCosThetaStar) fillME=kFALSE;
1410  }
1411  if(fillME){
1412  fMassVsPtVsYME->Fill(mass,pt,rapid);
1413  if(fFillHistosVsCosThetaStar) fMassVsPtVsCosthStME->Fill(mass,pt,absCosThSt);
1414  }
1415  }
1416  }
1417  return;
1418 }
1419 //________________________________________________________________________
1420 void AliAnalysisTaskCombinHF::FillMEHistosLS(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, Int_t charge){
1422 
1423  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
1424  Double_t pt = tmpRD->Pt();
1425  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
1426  Double_t mass=TMath::Sqrt(minv2);
1427 
1428  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
1429  Double_t rapid = tmpRD->Y(pdgD);
1430  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
1431  Bool_t fillME=kTRUE;
1432  Double_t costhst=0;
1433  Double_t absCosThSt=0;
1434  if(TMath::Abs(pdgD)==421 && (fApplyCutCosThetaStar || fFillHistosVsCosThetaStar)){
1435  costhst=tmpRD->CosThetaStar(0,421,321,211); // kaon is the first daughter
1436  absCosThSt=TMath::Abs(costhst);
1437  if(fApplyCutCosThetaStar && absCosThSt>fCutCosThetaStar) fillME=kFALSE;
1438  }
1439  if(fillME){
1440  if(charge>0) fMassVsPtVsYMELSpp->Fill(mass,pt,rapid);
1441  else if(charge<0) fMassVsPtVsYMELSmm->Fill(mass,pt,rapid);
1443  if(charge>0) fMassVsPtVsCosthStMELSpp->Fill(mass,pt,absCosThSt);
1444  else if(charge<0) fMassVsPtVsCosthStMELSmm->Fill(mass,pt,absCosThSt);
1445  }
1446  }
1447  }
1448  }
1449  return;
1450 }
1451 //________________________________________________________________________
1454 
1455  if(track->Charge()==0) return kFALSE;
1456  if(track->GetID()<0&&!fKeepNegID) return kFALSE;
1457  if(fFilterMask>=0){
1458  if(!(track->TestFilterMask(fFilterMask))) return kFALSE;
1459  }
1460  if(fCutTPCSignalN>0 && track->GetTPCsignalN()<fCutTPCSignalN) return kFALSE;
1461  if(!SelectAODTrack(track,fTrackCutsAll)) return kFALSE;
1462  return kTRUE;
1463 }
1464 
1465 //________________________________________________________________________
1468 
1469  if(!fPidHF) return kTRUE;
1470  Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
1471  Double_t mom=track->P();
1472  if(SelectAODTrack(track,fTrackCutsKaon)) {
1473  if(isKaon>=1) return kTRUE;
1474  if(isKaon<=-1) return kFALSE;
1475  switch(fPIDselCaseZero){// isKaon=0
1476  case 0:
1477  {
1478  return kTRUE;// always accept
1479  }
1480  break;
1481  case 1:
1482  {
1483  if(isKaon>=0 && track->P()>fmaxPforIDKaon)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDKaon
1484  }
1485  break;
1486  case 2:
1487  {
1488  if(track->P()>fmaxPforIDKaon)return kTRUE;
1489  AliPIDResponse *pidResp=fPidHF->GetPidResponse();// more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
1490  Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kKaon);
1491  if(nsigma>-2.&& nsigma<3. && mom<0.6)isKaon=1;
1492  else if(nsigma>-1.&& nsigma<3.&& mom<0.8)isKaon=1;
1493  if(isKaon==1)return kTRUE;
1494  }
1495  break;
1496  default:
1497  {
1498  AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
1499  return kFALSE;// actually case 0 could be set as the default and return kTRUE
1500  }
1501  }
1502  }
1503 
1504  return kFALSE;
1505 }
1506 //_______________________________________________________________________
1509 
1510  if(!fPidHF) return kTRUE;
1511  Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
1512  Double_t mom=track->P();
1513  if(SelectAODTrack(track,fTrackCutsPion)) {
1514  if(isPion>=1) return kTRUE;
1515  if(isPion<=-1) return kFALSE;
1516  switch(fPIDselCaseZero){// isPion=0
1517  case 0:
1518  {
1519  return kTRUE;// always accept
1520  }
1521  break;
1522  case 1:
1523  {
1524  if(track->P()>fmaxPforIDPion)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDPion
1525  }
1526  break;
1527  case 2:
1528  {
1529  // more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
1530  if(track->P()>fmaxPforIDPion)return kTRUE;
1531  AliPIDResponse *pidResp=fPidHF->GetPidResponse();
1532  Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kPion);
1533  if(nsigma<2.&& nsigma>-3. && mom<0.6)isPion=1;
1534  else if(nsigma<1. && nsigma> -3. && mom<0.8)isPion=1;
1535  if(isPion==1)return kTRUE;
1536  }
1537  break;
1538  default:
1539  {
1540  AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
1541  return kFALSE;// actually case 0 could be set as the default and return kTRUE
1542  }
1543  }
1544  }
1545 
1546  return kFALSE;
1547 }
1548 
1549 //________________________________________________________________________
1550 Bool_t AliAnalysisTaskCombinHF::SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts){
1552 
1553  if(!cuts) return kTRUE;
1554  // conversion to ESD track no longer needed after updates in AliESDtrackCuts to deal with AOD tracks
1555  // AliESDtrack esdTrack(track);
1556  // // set the TPC cluster info
1557  // esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
1558  // esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
1559  // esdTrack.SetTPCPointsF(track->GetTPCNclsF());
1560  // if(!cuts->IsSelected(&esdTrack)) return kFALSE;
1561  if(!cuts->IsSelected(track)) return kFALSE;
1562  return kTRUE;
1563 }
1564 
1565 //_________________________________________________________________
1566 Bool_t AliAnalysisTaskCombinHF::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1568  for (Int_t iProng = 0; iProng<nProng; iProng++){
1569  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1570  if(!mcPartDaughter) return kFALSE;
1571  Double_t eta = mcPartDaughter->Eta();
1572  Double_t pt = mcPartDaughter->Pt();
1573  if (TMath::Abs(eta) > fEtaAccCut || pt < fPtAccCut) return kFALSE;
1574  }
1575  return kTRUE;
1576 }
1577 //_________________________________________________________________
1580  if(!fzVertPoolLims || !fMultPoolLims) return 0;
1581  Int_t theBinZ=TMath::BinarySearch(fNzVertPoolsLimSize,fzVertPoolLims,zvert);
1582  if(theBinZ<0 || theBinZ>=fNzVertPoolsLimSize) return -1;
1583  Int_t theBinM=TMath::BinarySearch(fNMultPoolsLimSize,fMultPoolLims,mult);
1584  if(theBinM<0 || theBinM>=fNMultPoolsLimSize) return -1;
1585  return fNMultPools*theBinZ+theBinM;
1586 }
1587 //_________________________________________________________________
1590  if(poolIndex<0 || poolIndex>=fNOfPools) return;
1591  delete fEventBuffer[poolIndex];
1592  fEventBuffer[poolIndex]=new TTree(Form("EventBuffer_%d",poolIndex), "Temporary buffer for event mixing");
1593  fEventBuffer[poolIndex]->Branch("zVertex", &fVtxZ);
1594  fEventBuffer[poolIndex]->Branch("multiplicity", &fMultiplicity);
1595  fEventBuffer[poolIndex]->Branch("eventInfo", "TObjString",&fEventInfo);
1596  fEventBuffer[poolIndex]->Branch("karray", "TObjArray", &fKaonTracks);
1597  fEventBuffer[poolIndex]->Branch("parray", "TObjArray", &fPionTracks);
1598  return;
1599 }
1600 //_________________________________________________________________
1603  if(TMath::Abs(zv2-zv1)>fMaxzVertDistForMix) return kFALSE;
1604  if(TMath::Abs(mult2-mult1)>fMaxMultDiffForMix) return kFALSE;
1605  return kTRUE;
1606 }
1607 //_________________________________________________________________
1610 
1611  if(fDoEventMixing==0) return;
1612  Int_t nEvents=fEventBuffer[0]->GetEntries();
1613  if(fDebug > 1) printf("AnalysisTaskCombinHF::DoMixingWithCuts Start Event Mixing of %d events\n",nEvents);
1614 
1615  TObjArray* karray=0x0;
1616  TObjArray* parray=0x0;
1617  Double_t zVertex,mult;
1618  TObjString* eventInfo=0x0;
1619  fEventBuffer[0]->SetBranchAddress("karray", &karray);
1620  fEventBuffer[0]->SetBranchAddress("parray", &parray);
1621  fEventBuffer[0]->SetBranchAddress("eventInfo",&eventInfo);
1622  fEventBuffer[0]->SetBranchAddress("zVertex", &zVertex);
1623  fEventBuffer[0]->SetBranchAddress("multiplicity", &mult);
1624  Double_t d02[2]={0.,0.};
1625  Double_t d03[3]={0.,0.,0.};
1626  AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
1627  AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
1628  UInt_t pdg0[2]={321,211};
1629  // UInt_t pdgp[3]={321,211,211};
1630  Double_t px[3],py[3],pz[3];
1631  Int_t evId1,esdId1,nk1,np1;
1632  Int_t evId2,esdId2,nk2,np2;
1633 
1634  for(Int_t iEv1=0; iEv1<nEvents; iEv1++){
1635  fEventBuffer[0]->GetEvent(iEv1);
1636  TObjArray* karray1=(TObjArray*)karray->Clone();
1637  Double_t zVertex1=zVertex;
1638  Double_t mult1=mult;
1639  Int_t nKaons=karray1->GetEntries();
1640  Int_t nPionsForCheck=parray->GetEntries();
1641  sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId1,&esdId1,&np1,&nk1);
1642  if(nk1!=nKaons || np1!=nPionsForCheck){
1643  printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: read event does not match to the stored one\n");
1644  delete karray1;
1645  continue;
1646  }
1647  for(Int_t iEv2=0; iEv2<fNumberOfEventsForMixing; iEv2++){
1648  Int_t iToMix=iEv1+iEv2+1;
1649  if(iEv1>=(nEvents-fNumberOfEventsForMixing)) iToMix=iEv1-iEv2-1;
1650  if(iToMix<0) continue;
1651  if(iToMix==iEv1) continue;
1652  if(iToMix<iEv1) continue;
1653  fEventBuffer[0]->GetEvent(iToMix);
1654  Double_t zVertex2=zVertex;
1655  Double_t mult2=mult;
1656  if(TMath::Abs(zVertex2-zVertex1)<0.0001 && TMath::Abs(mult2-mult1)<0.001){
1657  printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: same event in mixing??? %d %d %f %f %f %f\n",iEv1,iEv2,zVertex1,zVertex2,mult1,mult2);
1658  continue;
1659  }
1660  TObjArray* parray2=(TObjArray*)parray->Clone();
1661  Int_t nPions=parray2->GetEntries();
1662  Int_t nKaonsForCheck=karray->GetEntries();
1663  sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId2,&esdId2,&np2,&nk2);
1664  if(nk2!=nKaonsForCheck || np2!=nPions){
1665  printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: read event does not match to the stored one\n");
1666  delete parray2;
1667  continue;
1668  }
1669  if(evId2==evId1 && esdId2==esdId1){
1670  printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: same event in mixing??? %d %d nK=%d %d nPi=%d %d\n",evId1,evId2,nKaons,nKaonsForCheck,nPionsForCheck,nPions);
1671  delete parray2;
1672  continue;
1673  }
1674  if(CanBeMixed(zVertex1,zVertex2,mult1,mult2)){
1675  for(Int_t iTr1=0; iTr1<nKaons; iTr1++){
1676  TLorentzVector* trK=(TLorentzVector*)karray1->At(iTr1);
1677  Double_t chargeK=trK->T();
1678  px[0] = trK->Px();
1679  py[0] = trK->Py();
1680  pz[0] = trK->Pz();
1681  for(Int_t iTr2=0; iTr2<nPions; iTr2++){
1682  TLorentzVector* trPi1=(TLorentzVector*)parray2->At(iTr2);
1683  Double_t chargePi1=trPi1->T();
1684  px[1] = trPi1->Px();
1685  py[1] = trPi1->Py();
1686  pz[1] = trPi1->Pz();
1687  if(fMeson==kDzero && chargePi1*chargeK<0){
1688  FillMEHistos(421,2,tmpRD2,px,py,pz,pdg0);
1689  }
1690  if(fMeson==kDzero && chargePi1*chargeK>0){
1691  FillMEHistosLS(421,2,tmpRD2,px,py,pz,pdg0,(Int_t)chargePi1);
1692  }
1693  }
1694  }
1695  }
1696  delete parray2;
1697  }
1698  delete karray1;
1699  }
1700  delete tmpRD2;
1701  delete tmpRD3;
1702 }
1703 //_________________________________________________________________
1706 
1707  if(fDoEventMixing==0) return;
1708  if(poolIndex<0 || poolIndex>fNzVertPools*fNMultPools) return;
1709 
1710  Int_t nEvents=fEventBuffer[poolIndex]->GetEntries();
1711  if(fDebug > 1) printf("AliAnalysisTaskCombinHF::DoMixingWithPools Start Event Mixing of %d events\n",nEvents);
1712  TObjArray* karray=0x0;
1713  TObjArray* parray=0x0;
1714  Double_t zVertex,mult;
1715  TObjString* eventInfo=0x0;
1716  fEventBuffer[poolIndex]->SetBranchAddress("karray", &karray);
1717  fEventBuffer[poolIndex]->SetBranchAddress("parray", &parray);
1718  fEventBuffer[poolIndex]->SetBranchAddress("eventInfo",&eventInfo);
1719  fEventBuffer[poolIndex]->SetBranchAddress("zVertex", &zVertex);
1720  fEventBuffer[poolIndex]->SetBranchAddress("multiplicity", &mult);
1721 
1722  // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
1723  Double_t d02[2]={0.,0.};
1724  Double_t d03[3]={0.,0.,0.};
1725  AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
1726  AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
1727  UInt_t pdg0[2]={321,211};
1728  UInt_t pdgp[3]={321,211,211};
1729  UInt_t pdgs[3]={321,211,321};
1730  Double_t px[3],py[3],pz[3];
1731  Int_t evId1,esdId1,nk1,np1;
1732  Int_t evId2,esdId2,nk2,np2;
1733 
1734  for(Int_t iEv1=0; iEv1<nEvents; iEv1++){
1735  fEventBuffer[poolIndex]->GetEvent(iEv1);
1736  TObjArray* karray1=(TObjArray*)karray->Clone();
1737  Double_t zVertex1=zVertex;
1738  Double_t mult1=mult;
1739  Int_t nKaons=karray1->GetEntries();
1740  Int_t nPionsForCheck=parray->GetEntries();
1741  sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId1,&esdId1,&np1,&nk1);
1742  if(nk1!=nKaons || np1!=nPionsForCheck){
1743  printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: read event does not match to the stored one\n");
1744  delete karray1;
1745  continue;
1746  }
1747  for(Int_t iEv2=0; iEv2<nEvents; iEv2++){
1748  if(iEv2==iEv1) continue;
1749  fEventBuffer[poolIndex]->GetEvent(iEv2);
1750  Double_t zVertex2=zVertex;
1751  Double_t mult2=mult;
1752  if(TMath::Abs(zVertex2-zVertex1)<0.0001 && TMath::Abs(mult2-mult1)<0.001){
1753  printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: same event in mixing??? %d %d %f %f %f %f\n",iEv1,iEv2,zVertex1,zVertex2,mult1,mult2);
1754  continue;
1755  }
1756  TObjArray* parray2=(TObjArray*)parray->Clone();
1757  Int_t nPions=parray2->GetEntries();
1758  Int_t nKaonsForCheck=karray->GetEntries();
1759  sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId2,&esdId2,&np2,&nk2);
1760  if(nk2!=nKaonsForCheck || np2!=nPions){
1761  printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: read event does not match to the stored one\n");
1762  delete parray2;
1763  continue;
1764  }
1765  if(evId2==evId1 && esdId2==esdId1){
1766  printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: same event in mixing??? %d %d nK=%d %d nPi=%d %d\n",evId1,evId2,nKaons,nKaonsForCheck,nPionsForCheck,nPions);
1767  delete parray2;
1768  continue;
1769  }
1770  TObjArray* parray3=0x0;
1771  Int_t nPions3=0;
1772  if(fMeson!=kDzero && fMeson!=kDs){
1773  Int_t iEv3=iEv2+1;
1774  if(iEv3==iEv1) iEv3=iEv2+2;
1775  if(iEv3>=nEvents) iEv3=iEv2-3;
1776  if(nEvents==2) iEv3=iEv1;
1777  if(iEv3<0) iEv3=iEv2-1;
1778  fEventBuffer[poolIndex]->GetEvent(iEv3);
1779  parray3=(TObjArray*)parray->Clone();
1780  nPions3=parray3->GetEntries();
1781  }
1782  for(Int_t iTr1=0; iTr1<nKaons; iTr1++){
1783  TLorentzVector* trK=(TLorentzVector*)karray1->At(iTr1);
1784  Double_t chargeK=trK->T();
1785  px[0] = trK->Px();
1786  py[0] = trK->Py();
1787  pz[0] = trK->Pz();
1788  for(Int_t iTr2=0; iTr2<nPions; iTr2++){
1789  TLorentzVector* trPi1=(TLorentzVector*)parray2->At(iTr2);
1790  Double_t chargePi1=trPi1->T();
1791  px[1] = trPi1->Px();
1792  py[1] = trPi1->Py();
1793  pz[1] = trPi1->Pz();
1794  if(chargePi1*chargeK<0){
1795  if(fMeson==kDzero){
1796  FillMEHistos(421,2,tmpRD2,px,py,pz,pdg0);
1797  }else if(fMeson==kDs) {
1798  for(Int_t iTr3=iTr1+1; iTr3<nKaons; iTr3++){
1799  TLorentzVector* trK2=(TLorentzVector*)karray1->At(iTr3);
1800  Double_t chargeK2=trK2->T();
1801  px[2] = trK2->Px();
1802  py[2] = trK2->Py();
1803  pz[2] = trK2->Pz();
1804  Double_t massKK=ComputeInvMassKK(trK,trK2);
1805  Double_t deltaMass=massKK-TDatabasePDG::Instance()->GetParticle(333)->Mass();
1806  Double_t cos1=CosPiKPhiRFrame(trK,trK2,trPi1);
1807  Double_t kincutPiKPhi=TMath::Abs(cos1*cos1*cos1);
1808  Double_t cosPiDsLabFrame=CosPiDsLabFrame(trK,trK2,trPi1);
1809  if(chargeK2*chargeK<0 && TMath::Abs(deltaMass)<fPhiMassCut && kincutPiKPhi>fCutCos3PiKPhiRFrame && cosPiDsLabFrame<fCutCosPiDsLabFrame){
1810  FillMEHistos(431,3,tmpRD3,px,py,pz,pdgs);
1811  }
1812  }
1813  } else{
1814  if(parray3){
1815  for(Int_t iTr3=iTr2+1; iTr3<nPions3; iTr3++){
1816  TLorentzVector* trPi2=(TLorentzVector*)parray3->At(iTr3);
1817  Double_t chargePi2=trPi2->T();
1818  px[2] = trPi2->Px();
1819  py[2] = trPi2->Py();
1820  pz[2] = trPi2->Pz();
1821  if(chargePi2*chargeK<0){
1822  if(fMeson==kDplus) FillMEHistos(411,3,tmpRD3,px,py,pz,pdgp);
1823  }
1824  }
1825  }
1826  }
1827  }else if(chargePi1*chargeK>0){
1828  if(fMeson==kDzero) FillMEHistosLS(421,2,tmpRD2,px,py,pz,pdg0,(Int_t)chargePi1);
1829  }
1830  }
1831  }
1832  delete parray3;
1833  delete parray2;
1834  }
1835  delete karray1;
1836  }
1837  delete tmpRD2;
1838  delete tmpRD3;
1839 }
1840 //_________________________________________________________________
1842 {
1844  if(fDoEventMixing==0) return;
1845  printf("AliAnalysisTaskCombinHF: FinishTaskOutput\n");
1846 
1847  if(fDoEventMixing==1){
1848  for(Int_t i=0; i<fNOfPools; i++){
1849  Int_t nEvents=fEventBuffer[i]->GetEntries();
1850  if(nEvents>1) DoMixingWithPools(i);
1851  }
1852  }else if(fDoEventMixing==2){
1853  DoMixingWithCuts();
1854  }
1855 }
1856 //_________________________________________________________________
1857 Double_t AliAnalysisTaskCombinHF::ComputeInvMassKK(AliAODTrack* tr1, AliAODTrack* tr2) const{
1859  Double_t massK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
1860  Double_t p1=tr1->P();
1861  Double_t p2=tr2->P();
1862  Double_t pxtot=tr1->Px()+tr2->Px();
1863  Double_t pytot=tr1->Py()+tr2->Py();
1864  Double_t pztot=tr1->Pz()+tr2->Pz();
1865  Double_t e1=TMath::Sqrt(massK*massK+p1*p1);
1866  Double_t e2=TMath::Sqrt(massK*massK+p2*p2);
1867  Double_t etot=e1+e2;
1868  Double_t m2=etot*etot-(pxtot*pxtot+pytot*pytot+pztot*pztot);
1869  return TMath::Sqrt(m2);
1870 }
1871 //_________________________________________________________________
1872 Double_t AliAnalysisTaskCombinHF::ComputeInvMassKK(TLorentzVector* tr1, TLorentzVector* tr2) const{
1874  Double_t massK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
1875  Double_t p1=tr1->P();
1876  Double_t p2=tr2->P();
1877  Double_t pxtot=tr1->Px()+tr2->Px();
1878  Double_t pytot=tr1->Py()+tr2->Py();
1879  Double_t pztot=tr1->Pz()+tr2->Pz();
1880  Double_t e1=TMath::Sqrt(massK*massK+p1*p1);
1881  Double_t e2=TMath::Sqrt(massK*massK+p2*p2);
1882  Double_t etot=e1+e2;
1883  Double_t m2=etot*etot-(pxtot*pxtot+pytot*pytot+pztot*pztot);
1884  return TMath::Sqrt(m2);
1885 }
1886 
1887 //----------------------------------------------------------------------
1888 Double_t AliAnalysisTaskCombinHF::CosPiKPhiRFrame(TLorentzVector* dauK1, TLorentzVector* dauK2, TLorentzVector* daupi) const {
1890 
1891  Double_t massK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
1892  Double_t masspi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
1893 
1894  Double_t eK1 = TMath::Sqrt(massK*massK+dauK1->P()*dauK1->P());
1895  Double_t eK2 = TMath::Sqrt(massK*massK+dauK2->P()*dauK2->P());
1896  Double_t epi = TMath::Sqrt(masspi*masspi+daupi->P()*daupi->P());
1897 
1898  Double_t ePhi=eK1+eK2;
1899  Double_t pxPhi=dauK1->Px()+dauK2->Px();
1900  Double_t pyPhi=dauK1->Py()+dauK2->Py();
1901  Double_t pzPhi=dauK1->Pz()+dauK2->Pz();
1902  Double_t bxPhi=pxPhi/ePhi;
1903  Double_t byPhi=pyPhi/ePhi;
1904  Double_t bzPhi=pzPhi/ePhi;
1905 
1906  TVector3 vecK1Phiframe;
1907  TLorentzVector vecK1(dauK1->Px(),dauK1->Py(),dauK1->Pz(),eK1);
1908  vecK1.Boost(-bxPhi,-byPhi,-bzPhi);
1909  vecK1.Boost(vecK1Phiframe);
1910  vecK1Phiframe=vecK1.BoostVector();
1911 
1912  TVector3 vecPiPhiframe;
1913  TLorentzVector vecPi(daupi->Px(),daupi->Py(),daupi->Pz(),epi);
1914  vecPi.Boost(-bxPhi,-byPhi,-bzPhi);
1915  vecPi.Boost(vecPiPhiframe);
1916  vecPiPhiframe=vecPi.BoostVector();
1917 
1918  Double_t innera=vecPiPhiframe.Dot(vecK1Phiframe);
1919  Double_t norm1a=TMath::Sqrt(vecPiPhiframe.Dot(vecPiPhiframe));
1920  Double_t norm2a=TMath::Sqrt(vecK1Phiframe.Dot(vecK1Phiframe));
1921  Double_t cosK1PhiFrame=innera/(norm1a*norm2a);
1922 
1923  return cosK1PhiFrame;
1924 }
1925 
1926 //----------------------------------------------------------------------
1927 Double_t AliAnalysisTaskCombinHF::CosPiDsLabFrame(TLorentzVector* dauK1, TLorentzVector* dauK2, TLorentzVector* daupi) const {
1929 
1930  Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1931  Double_t masspi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
1932 
1933  Double_t pxD=dauK1->Px()+dauK2->Px()+daupi->Px();
1934  Double_t pyD=dauK1->Px()+dauK2->Px()+daupi->Px();
1935  Double_t pzD=dauK1->Px()+dauK2->Px()+daupi->Px();
1936  Double_t pD=TMath::Sqrt(pxD*pxD+pyD*pyD+pzD*pzD);
1937  Double_t eD=TMath::Sqrt(massDs*massDs+pD*pD);
1938  Double_t epi = TMath::Sqrt(masspi*masspi+daupi->P()*daupi->P());
1939  Double_t bxD=pxD/eD;
1940  Double_t byD=pyD/eD;
1941  Double_t bzD=pzD/eD;
1942 
1943  TVector3 piDsframe;
1944  TLorentzVector vecPi(daupi->Px(),daupi->Py(),daupi->Pz(),epi);
1945  vecPi.Boost(-bxD,-byD,-bzD);
1946  vecPi.Boost(piDsframe);
1947  piDsframe=vecPi.BoostVector();
1948 
1949  TVector3 vecDs(pxD,pyD,pzD);
1950 
1951  Double_t inner=vecDs.Dot(piDsframe);
1952  Double_t norm1=TMath::Sqrt(vecDs.Dot(vecDs));
1953  Double_t norm2=TMath::Sqrt(piDsframe.Dot(piDsframe));
1954  Double_t cosPiDsFrame=inner/(norm1*norm2);
1955 
1956  return cosPiDsFrame;
1957 }
1958 
1959 //_________________________________________________________________
1961 {
1963  //
1964  if(fDebug > 1) printf("AliAnalysisTaskCombinHF: Terminate() \n");
1965  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1966  if (!fOutput) {
1967  printf("ERROR: fOutput not available\n");
1968  return;
1969  }
1970  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
1971  if(fHistNEvents){
1972  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1973  }else{
1974  printf("ERROR: fHistNEvents not available\n");
1975  return;
1976  }
1977  return;
1978 }
1979 
Int_t charge
TH3F * fHistoPtKPtPiPtD
! hist. for propagation of tracking unc
Bool_t IsEventRejectedDueToCentrality() const
Definition: AliRDHFCuts.h:361
Double_t fEtaAccCut
width of pt bin (GeV/c)
TH3F * fPtVsYVsPtBGenLargeAccFeeddw
! hist. of Y vs. Pt vs. PtB generated (|y|<0.9)
Bool_t IsEventRejectedDueToZVertexOutsideFiducialRegion() const
Definition: AliRDHFCuts.h:355
Bool_t IsTrackSelected(AliAODTrack *track)
static Int_t CheckD0Decay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
Bool_t IsEventRejectedDueToNotRecoVertex() const
Definition: AliRDHFCuts.h:346
double Double_t
Definition: External.C:58
Definition: External.C:260
Double_t fMaxzVertDistForMix
maximum number of events to be used in event mixing
void FillMEHistosLS(Int_t pdgD, Int_t nProngs, AliAODRecoDecay *tmpRD, Double_t *px, Double_t *py, Double_t *pz, UInt_t *pdgdau, Int_t charge)
Int_t fNzVertPoolsLimSize
number of pools in z vertex for event mixing
void StoreCandidates(AliVEvent *, Int_t nCand=0, Bool_t flagFilter=kTRUE)
TH3F * fMassVsPtVsYMELSpp
! hist. of Y vs. Pt vs. Mass (mixedevents)
Bool_t IsKaon(AliAODTrack *track)
Definition: External.C:236
TH2F * fHistCheckOriginRecoD
!hist. of origin (c/b) of D meson (reco)
Int_t GetnSigmaTOF(AliAODTrack *track, Int_t species, Double_t &sigma) const
Bool_t CheckAcceptance(TClonesArray *arrayMC, Int_t nProng, Int_t *labDau)
static Int_t CheckDplusDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
TH1F * fHistTrackStatus
!hist. of status of tracks
AliNormalizationCounter * fCounter
maximum angle for track rotation (3rd prong)
Bool_t IsPion(AliAODTrack *track)
Int_t IsEventSelectedInCentrality(AliVEvent *event)
Int_t fNMultPoolsLimSize
number of pools in multiplicity for event mixing
Double_t fCutCos3PiKPhiRFrame
cut on the KK inv mass for phi selection
void FillMEHistos(Int_t pdgD, Int_t nProngs, AliAODRecoDecay *tmpRD, Double_t *px, Double_t *py, Double_t *pz, UInt_t *pdgdau)
TH3F * fPtVsYVsMultGenLargeAccPrompt
! hist. of Y vs. Pt vs. Mult generated (|y|<0.9)
static Int_t CheckDsDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
void ConfigureMultiplicityPools(Int_t nPools, Double_t *multLimits)
Double_t fMaxMultiplicity
lower limit for multiplcities in MC histos
TH1F * fNSelected
! hist. of n. of selected D+
TList * fOutput
! list with output histograms
Bool_t fReadMC
mesonSpecies (see enum)
TH3F * fMassVsPtVsY
! hist. of Y vs. Pt vs. Mass (all cand)
Bool_t FillHistos(Int_t pdgD, Int_t nProngs, AliAODRecoDecay *tmpRD, Double_t *px, Double_t *py, Double_t *pz, UInt_t *pdgdau, TClonesArray *arrayMC, AliAODMCHeader *mcHeader, Int_t *dgLabels)
Int_t GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &sigma) const
TH2F * fHistEventMultZvEvSel
!hist. of evnt Mult vs. Zv for selected ev
Double_t fMinAngleForRot
number of rotations
TH3F * fMassVsPtVsCosthStLSmm
! hist. of Pt vs. Mass vs. cos(th*) (like sign –)
Double_t mass
Bool_t CheckTOFPIDStatus(AliAODTrack *track) const
TH3F * fPtVsYVsMultGenPrompt
! hist. of Y vs. Pt vs. Mult generated (all D)
TH3F * fPtVsYVsMultGenLimAccPrompt
! hist. of Y vs. Pt vs. Mult generated (|y|<0.5)
TH3F * fPtVsYVsMultRecoFeeddw
! hist. of Y vs. Pt vs. Mult generated (Reco D)
Double_t fPtAccCut
eta limits for acceptance step
Double_t fMinPtHard
flag to select specific phard range in MC
TH3F * fPtVsYVsMultRecoPrompt
! hist. of Y vs. Pt vs. Mult generated (Reco D)
TH3F * fPtVsYVsPtBGenLimAccFeeddw
! hist. of Y vs. Pt vs. PtB generated (|y|<0.5)
Bool_t IsEventRejectedDueToVertexContributors() const
Definition: AliRDHFCuts.h:349
TH3F * fPtVsYVsMultGenLargeAccFeeddw
! hist. of Y vs. Pt vs. Mult generated (|y|<0.9)
Double_t fPtBinWidth
maximum pT value for inv. mass histograms
TH3F * fPtVsYVsPtBGenAccFeeddw
! hist. of Y vs. Pt vs. PtB generated (D in acc)
void ConfigureZVertPools(Int_t nPools, Double_t *zVertLimits)
Bool_t fKeepNegID
flag for upper p limit for id band for kaon
TTree ** fEventBuffer
number of pools
TH3F * fMassVsPtVsCosthStMELSmm
! hist. of Pt vs. Mass vs. cos(th*) (mixedevents)
Double_t fBayesThresPion
threshold for kaon identification via Bayesian PID
Double_t fMaxPt
maximum value of invariant mass
Int_t fNRotations3
maximum angle for track rotation
Bool_t fEnforceMBTrigMaskInMC
flag for access to MC
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Functions to check the decay tree.
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:274
Double_t CosPiDsLabFrame(TLorentzVector *dauK1, TLorentzVector *dauK2, TLorentzVector *daupi) const
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:264
TH2F * fHistEventMultZv
!hist. of evnt Mult vs. Zv for all events
TH3F * fMassVsPtVsCosthStSig
! hist. of Pt vs. Mass vs. cos(th*) (signal)
Bool_t GetUseTimeRangeCutForPbPb2018() const
Definition: AliRDHFCuts.h:310
Bool_t CanBeMixed(Double_t zv1, Double_t zv2, Double_t mult1, Double_t mult2)
Double_t fPhiMassCut
cos(theta*) cut
Bool_t fApplyCutCosThetaStar
flag to control cos(theta*) cut
TH2F * fHistonSigmaTPCKaon
! hist. of nSigmaTPC kaon
TH2F * fHistonSigmaTPCKaonGoodTOF
! hist. of nSigmaTPC kaon
static Int_t GetNumberOfTrackletsInEtaRange(AliAODEvent *ev, Double_t mineta, Double_t maxeta)
Int_t fPIDselCaseZero
flag to keep also track with negative ID (default kFALSE, change it only if you know what you are doi...
TH1F * fHistNEvents
!hist. for No. of events
TH1F * fNormRotated
! hist. rotated/selected D+
TH3F * fMassVsPtVsCosthStBkg
! hist. of Pt vs. Mass vs. cos(th*) (background)
TH3F * fMassVsPtVsYLSmm
! hist. of Y vs. Pt vs. Mass (like sign –)
TH3F * fMassVsPtVsCosthStMELSpp
! hist. of Pt vs. Mass vs. cos(th*) (mixedevents)
int Int_t
Definition: External.C:63
TH1F * fHistXsecVsPtHard
!hist. of xsec vs pthard (MC)
TList * fListCuts
! list with cut values
Bool_t fFillHistosVsCosThetaStar
min. value of number of TPC clusters for PID, cut if !=0
Double_t CosPiKPhiRFrame(TLorentzVector *dauK1, TLorentzVector *dauK2, TLorentzVector *daupi) const
unsigned int UInt_t
Definition: External.C:33
AliPIDCombined * GetPidCombined() const
Definition: AliAODPidHF.h:174
TH2F * fHistonSigmaTOFPion
! hist. of nSigmaTOF pion
float Float_t
Definition: External.C:68
TH1F * fBMohterPtGen
! hist. of beauty mother pt
Double_t fMaxAngleForRot
minimum angle for track rotation
Int_t GetPoolIndex(Double_t zvert, Double_t mult)
Double_t fMaxPtHard
minimum pthard
Double_t fMinAngleForRot3
number of rotations (3rd prong)
AliESDtrackCuts * fTrackCutsAll
FilterMask.
TH2F * fHistCheckOrigin
!hist. of origin (c/b) of D meson (gen)
TH3F * fPtVsYVsMultGenAccEvSelFeeddw
! hist. of Y vs. Pt vs. Mult generated (D in acc, sel ev.)
TH3F * fMassVsPtVsCosthSt
! hist. of Pt vs. Mass vs. cos(th*) (all cand)
AliRDHFCuts * fAnalysisCuts
PID configuration.
static Bool_t IsTrackFromHadronDecay(Int_t pdgMoth, AliAODTrack *tr, TClonesArray *arrayMC)
Bool_t IsEventRejectedDueToBadTrackVertex() const
Definition: AliRDHFCuts.h:373
THnSparse * fDeltaMassFullAnalysis
! hist. mass difference after rotations with more details
Int_t fNzVertPools
cut on multiplicity difference for event mixing with cuts
TH3F * fMassVsPtVsYSig
! hist. of Y vs. Pt vs. Mass (signal)
Int_t fFullAnalysis
flag for definition of c,b origin
void FillGenHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader, Bool_t isEvSel)
Double_t * fzVertPoolLims
number of pools in z vertex for event mixing +1
TH3F * fPtVsYVsMultGenAccEvSelPrompt
! hist. of Y vs. Pt vs. Mult generated (D in acc, sel ev.)
TH3F * fMassVsPtVsCosthStLSpp
! hist. of Pt vs. Mass vs. cos(th*) (like sign ++)
TH3F * fMassVsPtVsCosthStME
! hist. of Pt vs. Mass vs. cos(th*) (mixedevents)
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:275
Bool_t IsEventRejectedDueToPileup() const
Definition: AliRDHFCuts.h:358
Double_t nsigma
Double_t nEvents
plot quality messages
TH3F * fMassVsPtVsCosthStRot
! hist. of Pt vs. Mass vs. cos(th*) (rotations)
void FillLSHistos(Int_t pdgD, Int_t nProngs, AliAODRecoDecay *tmpRD, Double_t *px, Double_t *py, Double_t *pz, UInt_t *pdgdau, Int_t charge)
TH3F * fPtVsYVsMultGenAccFeeddw
! hist. of Y vs. Pt vs. Mult generated (D in acc)
virtual void UserExec(Option_t *option)
Int_t MakeRawPid(AliAODTrack *track, Int_t specie)
AliPIDResponse * GetPidResponse() const
Definition: AliAODPidHF.h:173
virtual void Terminate(Option_t *option)
TH3F * fPtVsYVsMultGenFeeddw
! hist. of Y vs. Pt vs. Mult generated (all D)
Double_t fCutCosThetaStar
flag to control cos(theta*) cut
TH3F * fPtVsYVsPtBRecoFeeddw
! hist. of Y vs. Pt vs. PtB generated (Reco D)
Double_t ComputeInvMassKK(AliAODTrack *tr1, AliAODTrack *tr2) const
Double_t fMinMultiplicity
number of bins for multiplcities in MC histos
TH3F * fMassVsPtVsCosthStRefl
! hist. of Pt vs. Mass vs. cos(th*) (reflections)
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Bool_t IsEventRejectedDueToBadPrimaryVertex() const
Definition: AliRDHFCuts.h:382
Bool_t fSelectPtHardRange
flag to speed up the MC
Bool_t IsEventRejectedDuePhysicsSelection() const
Definition: AliRDHFCuts.h:379
TH3F * fHistoPtKPtPiPtDSig
! hist. for propagation of tracking unc
TH2F * fHistEventMultCentEvSel
!hist. for evnt Mult vs. centrality (sel)
Double_t fMinMass
Cuts for candidates.
TH2F * fHistCheckOriginRecoVsGen
!hist. of origin (c/b) of D meson
Int_t fPIDstrategy
maximum pthard
Double_t fVtxZ
unique event Id for event mixing checks
AliESDtrackCuts * fTrackCutsKaon
pion track selection
static Double_t GetBeautyMotherPt(TClonesArray *arrayMC, AliAODMCParticle *mcPart)
TH2F * fHistonSigmaTPCPionGoodTOF
! hist. of nSigmaTPC pion
Bool_t fGoUpToQuark
if true force the MC to use
Bool_t SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts)
TH3F * fPtVsYVsMultGenLimAccFeeddw
! hist. of Y vs. Pt vs. Mult generated (|y|<0.5)
Bool_t IsEventSelected(AliVEvent *event)
Bool_t fSignalOnlyMC
flag to set analysis level (0 is the fastest)
TH1F * fHistCheckDecChanAcc
!hist. of decay channel of D meson in acc.
Int_t fDoEventMixing
container for checks
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
TH3F * fMassVsPtVsYME
! hist. of Y vs. Pt vs. Mass (mixedevents)
TH3F * fPtVsYVsPtBGenFeeddw
! hist. of Y vs. Pt vs. PtB generated (all D)
Double_t fmaxPforIDKaon
flag for upper p limit for id band for pion
static Bool_t IsTrackInjected(Int_t label, AliAODMCHeader *header, TClonesArray *arrayMC)
TH3F * fPtVsYVsMultGenAccPrompt
! hist. of Y vs. Pt vs. Mult generated (D in acc)
TH2F * fEventsPerPool
! hist with number of events per pool
TH3F * fMassVsPtVsYBkg
! hist. of Y vs. Pt vs. Mass (background)
TH2F * fMixingsPerPool
! hist with number of mixings per pool
TH3F * fMassVsPtVsYLSpp
! hist. of Y vs. Pt vs. Mass (like sign ++)
Int_t fNumberOfEventsForMixing
flag for event mixing
TH2F * fHistonSigmaTOFKaon
! hist. of nSigmaTOF kaon
TObjArray * fPionTracks
array of kaon-compatible tracks (TLorentzVectors)
void DoMixingWithPools(Int_t poolIndex)
Int_t fOrigContainer[200000]
threshold for pion identification via Bayesian PID
Int_t fNRotations
pt limits for acceptance step
const char Option_t
Definition: External.C:48
TH2F * fHistEventMultCent
!hist. for evnt Mult vs. centrality (all)
Int_t fCutTPCSignalN
kaon track selection
Double_t fBayesThresKaon
flag to change PID strategy
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:343
TH1F * fDeltaMass
! hist. mass difference after rotations
TH2F * fHistonSigmaTPCPion
! hist. of nSigmaTPC pion
TObjArray * fKaonTracks
upper limit for multiplcities in MC histos
bool Bool_t
Definition: External.C:53
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:338
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:297
TH3F * fPtVsYVsPtBGenAccEvSelFeeddw
! hist. of Y vs. Pt vs. PtB generated (D in acc, sel ev.)
Double_t fMaxMultDiffForMix
cut on zvertex distance for event mixing with cuts
Double_t fMaxMass
minimum value of invariant mass
TString meson
TH3F * fMassVsPtVsYMELSmm
! hist. of Y vs. Pt vs. Mass (mixedevents)
TH1F * fHistCheckDecChan
!hist. of decay channel of D meson
Double_t fmaxPforIDPion
knSigma, kBayesianMaxProb, kBayesianThres
TH3F * fMassVsPtVsYRefl
! hist. of Y vs. Pt vs. Mass (reflections)
Double_t * fMultPoolLims
number of pools in multiplicity for event mixing +1
void SetTriggerMask(ULong64_t mask=0)
Definition: AliRDHFCuts.h:71
Double_t fMaxAngleForRot3
minimum angle for track rotation (3rd prong)
void SetPidResponse(AliPIDResponse *pidResp)
Definition: AliAODPidHF.h:124
TH3F * fMassVsPtVsYRot
! hist. of Y vs. Pt vs. Mass (rotations)