AliPhysics  5a28df1 (5a28df1)
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"
44 
46 ClassImp(AliAnalysisTaskCombinHF);
48 
49 //________________________________________________________________________
52  fOutput(0x0),
53  fHistNEvents(0x0),
54  fHistEventMultZv(0x0),
55  fHistEventMultZvEvSel(0x0),
56  fHistTrackStatus(0x0),
57  fHistTrackEtaMultZv(0x0),
58  fHistCheckOrigin(0x0),
59  fHistCheckDecChan(0x0),
60  fHistCheckDecChanAcc(0x0),
61  fPtVsYVsMultGenPrompt(0x0),
62  fPtVsYVsMultGenLargeAccPrompt(0x0),
63  fPtVsYVsMultGenLimAccPrompt(0x0),
64  fPtVsYVsMultGenAccPrompt(0x0),
65  fPtVsYVsMultGenAccEvSelPrompt(0x0),
66  fPtVsYVsMultRecoPrompt(0x0),
67  fPtVsYVsMultGenFeeddw(0x0),
68  fPtVsYVsMultGenLargeAccFeeddw(0x0),
69  fPtVsYVsMultGenLimAccFeeddw(0x0),
70  fPtVsYVsMultGenAccFeeddw(0x0),
71  fPtVsYVsMultGenAccEvSelFeeddw(0x0),
72  fPtVsYVsMultRecoFeeddw(0x0),
73  fMassVsPtVsY(0x0),
74  fMassVsPtVsYRot(0x0),
75  fMassVsPtVsYLSpp(0x0),
76  fMassVsPtVsYLSmm(0x0),
77  fMassVsPtVsYSig(0x0),
78  fMassVsPtVsYRefl(0x0),
79  fMassVsPtVsYBkg(0x0),
80  fNSelected(0x0),
81  fNormRotated(0x0),
82  fDeltaMass(0x0),
83  fDeltaMassFullAnalysis(0x0),
84  fMassVsPtVsYME(0x0),
85  fMassVsPtVsYMELSpp(0x0),
86  fMassVsPtVsYMELSmm(0x0),
87  fEventsPerPool(0x0),
88  fMixingsPerPool(0x0),
89  fFilterMask(BIT(4)),
90  fTrackCutsAll(0x0),
91  fTrackCutsPion(0x0),
92  fTrackCutsKaon(0x0),
93  fPhiMassCut(99999.),
94  fCutCos3PiKPhiRFrame(-1.1),
95  fCutCosPiDsLabFrame(1.1),
96  fPidHF(new AliAODPidHF()),
97  fAnalysisCuts(0x0),
98  fMinMass(1.720),
99  fMaxMass(2.150),
100  fMaxPt(10.),
101  fPtBinWidth(0.5),
102  fEtaAccCut(0.9),
103  fPtAccCut(0.1),
104  fNRotations(9),
105  fMinAngleForRot(5*TMath::Pi()/6),
106  fMaxAngleForRot(7*TMath::Pi()/6),
107  fNRotations3(9),
108  fMinAngleForRot3(2*TMath::Pi()/6),
109  fMaxAngleForRot3(4*TMath::Pi()/6),
110  fCounter(0x0),
111  fMeson(kDzero),
112  fReadMC(kFALSE),
113  fGoUpToQuark(kTRUE),
114  fFullAnalysis(0),
115  fPIDstrategy(knSigma),
116  fmaxPforIDPion(0.8),
117  fmaxPforIDKaon(2.),
118  fKeepNegID(kFALSE),
119  fPIDselCaseZero(0),
120  fBayesThresKaon(0.4),
121  fBayesThresPion(0.4),
122  fDoEventMixing(1),
123  fNumberOfEventsForMixing(20),
124  fMaxzVertDistForMix(5.),
125  fMaxMultDiffForMix(5.),
126  fNzVertPools(1),
127  fNzVertPoolsLimSize(2),
128  fzVertPoolLims(0x0),
129  fNMultPools(1),
130  fNMultPoolsLimSize(2),
131  fMultPoolLims(0x0),
132  fNOfPools(1),
133  fEventBuffer(0x0),
134  fEventInfo(new TObjString("")),
135  fVtxZ(0),
136  fMultiplicity(0),
137  fMinMultiplicity(-0.5),
138  fMaxMultiplicity(199.5),
139  fKaonTracks(0x0),
140  fPionTracks(0x0)
141 {
143 }
144 
145 //________________________________________________________________________
147  AliAnalysisTaskSE("DmesonCombin"),
148  fOutput(0x0),
149  fHistNEvents(0x0),
150  fHistEventMultZv(0x0),
152  fHistTrackStatus(0x0),
153  fHistTrackEtaMultZv(0x0),
154  fHistCheckOrigin(0x0),
155  fHistCheckDecChan(0x0),
169  fMassVsPtVsY(0x0),
170  fMassVsPtVsYRot(0x0),
171  fMassVsPtVsYLSpp(0x0),
172  fMassVsPtVsYLSmm(0x0),
173  fMassVsPtVsYSig(0x0),
174  fMassVsPtVsYRefl(0x0),
175  fMassVsPtVsYBkg(0x0),
176  fNSelected(0x0),
177  fNormRotated(0x0),
178  fDeltaMass(0x0),
180  fMassVsPtVsYME(0x0),
181  fMassVsPtVsYMELSpp(0x0),
182  fMassVsPtVsYMELSmm(0x0),
183  fEventsPerPool(0x0),
184  fMixingsPerPool(0x0),
185  fFilterMask(BIT(4)),
186  fTrackCutsAll(0x0),
187  fTrackCutsPion(0x0),
188  fTrackCutsKaon(0x0),
189  fPhiMassCut(99999.),
191  fCutCosPiDsLabFrame(1.1),
192  fPidHF(new AliAODPidHF()),
193  fAnalysisCuts(analysiscuts),
194  fMinMass(1.720),
195  fMaxMass(2.150),
196  fMaxPt(10.),
197  fPtBinWidth(0.5),
198  fEtaAccCut(0.9),
199  fPtAccCut(0.1),
200  fNRotations(9),
201  fMinAngleForRot(5*TMath::Pi()/6),
202  fMaxAngleForRot(7*TMath::Pi()/6),
203  fNRotations3(9),
204  fMinAngleForRot3(2*TMath::Pi()/6),
205  fMaxAngleForRot3(4*TMath::Pi()/6),
206  fCounter(0x0),
207  fMeson(meson),
208  fReadMC(kFALSE),
209  fGoUpToQuark(kTRUE),
210  fFullAnalysis(0),
212  fmaxPforIDPion(0.8),
213  fmaxPforIDKaon(2.),
214  fKeepNegID(kFALSE),
215  fPIDselCaseZero(0),
216  fBayesThresKaon(0.4),
217  fBayesThresPion(0.4),
218  fDoEventMixing(1),
221  fMaxMultDiffForMix(5.),
222  fNzVertPools(1),
224  fzVertPoolLims(0x0),
225  fNMultPools(1),
227  fMultPoolLims(0x0),
228  fNOfPools(1),
229  fEventBuffer(0x0),
230  fEventInfo(new TObjString("")),
231  fVtxZ(0),
232  fMultiplicity(0),
233  fMinMultiplicity(-0.5),
234  fMaxMultiplicity(199.5),
235  fKaonTracks(0x0),
236  fPionTracks(0x0)
237 {
239 
240  DefineOutput(1,TList::Class()); //My private output
241  DefineOutput(2,AliNormalizationCounter::Class());
242 }
243 
244 //________________________________________________________________________
246 {
247  //
249  //
250  if(fOutput && !fOutput->IsOwner()){
251  delete fHistNEvents;
252  delete fHistEventMultZv;
253  delete fHistEventMultZvEvSel;
254  delete fHistTrackStatus;
255  delete fHistTrackEtaMultZv;
256  delete fHistCheckOrigin;
257  delete fHistCheckDecChan;
258  delete fHistCheckDecChanAcc;
259  delete fPtVsYVsMultGenPrompt;
264  delete fPtVsYVsMultRecoPrompt;
265  delete fPtVsYVsMultGenFeeddw;
270  delete fPtVsYVsMultRecoFeeddw;
271  delete fMassVsPtVsY;
272  delete fMassVsPtVsYLSpp;
273  delete fMassVsPtVsYLSmm;
274  delete fMassVsPtVsYRot;
275  delete fMassVsPtVsYSig;
276  delete fMassVsPtVsYRefl;
277  delete fMassVsPtVsYBkg;
278  delete fNSelected;
279  delete fNormRotated;
280  delete fDeltaMass;
281  delete fDeltaMassFullAnalysis;
282  delete fMassVsPtVsYME;
283  delete fMassVsPtVsYMELSpp;
284  delete fMassVsPtVsYMELSmm;
285  }
286 
287  delete fOutput;
288  delete fCounter;
289  delete fTrackCutsAll;
290  delete fTrackCutsPion;
291  delete fTrackCutsKaon;
292  delete fPidHF;
293  delete fAnalysisCuts;
294  if(fKaonTracks) fKaonTracks->Delete();
295  if(fPionTracks) fPionTracks->Delete();
296  delete fKaonTracks;
297  delete fPionTracks;
298 
299  if(fEventBuffer){
300  for(Int_t i=0; i<fNOfPools; i++) delete fEventBuffer[i];
301  delete fEventBuffer;
302  }
303  delete fEventInfo;
304  delete [] fzVertPoolLims;
305  delete [] fMultPoolLims;
306 }
307 
308 //________________________________________________________________________
310 {
312  if(fzVertPoolLims) delete [] fzVertPoolLims;
313  fNzVertPools=nPools;
314  fNzVertPoolsLimSize=nPools+1;
316  for(Int_t ib=0; ib<fNzVertPoolsLimSize; ib++) fzVertPoolLims[ib]=zVertLimits[ib];
317  return;
318 }
319 //________________________________________________________________________
321 {
322  // sets the pools for event mizing in zvertex
323  if(fMultPoolLims) delete [] fMultPoolLims;
324  fNMultPools=nPools;
325  fNMultPoolsLimSize=nPools+1;
327  for(Int_t ib=0; ib<nPools+1; ib++) fMultPoolLims[ib]=multLimits[ib];
328  return;
329 }
330 //________________________________________________________________________
332 {
334  //
335  if(fDebug > 1) printf("AnalysisTaskCombinHF::UserCreateOutputObjects() \n");
336 
337  fOutput = new TList();
338  fOutput->SetOwner();
339  fOutput->SetName("OutputHistos");
340 
341  fHistNEvents = new TH1F("hNEvents", "number of events ",9,-0.5,8.5);
342  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
343  fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
344  fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
345  fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to phys sel");
346  fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected due to not reco vertex");
347  fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for contr vertex");
348  fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for vertex out of accept");
349  fHistNEvents->GetXaxis()->SetBinLabel(8,"n. rejected for pileup events");
350  fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of out centrality events");
351 
352  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
353  fHistNEvents->SetMinimum(0);
354  fOutput->Add(fHistNEvents);
355 
356  fHistEventMultZv = new TH2F("hEventMultZv","",30,-15.,15.,200,fMinMultiplicity,fMaxMultiplicity);
358 
359  fHistEventMultZvEvSel = new TH2F("hEventMultZvEvSel","",30,-15.,15.,200,fMinMultiplicity,fMaxMultiplicity);
361 
362  fHistTrackStatus = new TH1F("hTrackStatus", "",8,-0.5,7.5);
363  fHistTrackStatus->GetXaxis()->SetBinLabel(1,"Not OK");
364  fHistTrackStatus->GetXaxis()->SetBinLabel(2,"Track OK");
365  fHistTrackStatus->GetXaxis()->SetBinLabel(3,"Kaon, Not OK");
366  fHistTrackStatus->GetXaxis()->SetBinLabel(4,"Kaon OK");
367  fHistTrackStatus->GetXaxis()->SetBinLabel(5,"Pion, Not OK");
368  fHistTrackStatus->GetXaxis()->SetBinLabel(6,"Pion OK");
369  fHistTrackStatus->GetXaxis()->SetBinLabel(7,"Kaon||Pion, Not OK");
370  fHistTrackStatus->GetXaxis()->SetBinLabel(8,"Kaon||Pion OK");
371  fHistTrackStatus->GetXaxis()->SetNdivisions(1,kFALSE);
372  fHistTrackStatus->SetMinimum(0);
374 
375  fHistTrackEtaMultZv = new TH3F("hTrackEtaMultZv","",40,-1.,1.,30,-15.,15.,200,fMinMultiplicity,fMaxMultiplicity);
377 
378  Int_t nPtBins = (Int_t)(fMaxPt/fPtBinWidth+0.001);
380 
381  if(fReadMC){
382 
383  fHistCheckOrigin=new TH1F("hCheckOrigin","",7,-1.5,5.5);
384  fHistCheckOrigin->SetMinimum(0);
386 
387  fHistCheckDecChan=new TH1F("hCheckDecChan","",7,-2.5,4.5);
388  fHistCheckDecChan->SetMinimum(0);
390 
391  fHistCheckDecChanAcc=new TH1F("hCheckDecChanAcc","",7,-2.5,4.5);
392  fHistCheckDecChanAcc->SetMinimum(0);
394 
395  fPtVsYVsMultGenPrompt = new TH3F("hPtVsYVsMultGenPrompt","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
396  fPtVsYVsMultGenPrompt->Sumw2();
397  fPtVsYVsMultGenPrompt->SetMinimum(0);
399 
400  fPtVsYVsMultGenLargeAccPrompt = new TH3F("hPtVsYVsMultGenLargeAccPrompt","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
402  fPtVsYVsMultGenLargeAccPrompt->SetMinimum(0);
404 
405  fPtVsYVsMultGenLimAccPrompt = new TH3F("hPtVsYVsMultGenLimAccPrompt","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
407  fPtVsYVsMultGenLimAccPrompt->SetMinimum(0);
409 
410  fPtVsYVsMultGenAccPrompt = new TH3F("hPtVsYVsMultGenAccPrompt","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
411  fPtVsYVsMultGenAccPrompt->Sumw2();
412  fPtVsYVsMultGenAccPrompt->SetMinimum(0);
414 
415  fPtVsYVsMultGenAccEvSelPrompt = new TH3F("hPtVsYVsMultGenAccEvSelPrompt","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
417  fPtVsYVsMultGenAccEvSelPrompt->SetMinimum(0);
419 
420  fPtVsYVsMultRecoPrompt = new TH3F("hPtVsYVsMultRecoPrompt","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
421  fPtVsYVsMultRecoPrompt->Sumw2();
422  fPtVsYVsMultRecoPrompt->SetMinimum(0);
424 
425  fPtVsYVsMultGenFeeddw = new TH3F("hPtVsYVsMultGenFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
426  fPtVsYVsMultGenFeeddw->Sumw2();
427  fPtVsYVsMultGenFeeddw->SetMinimum(0);
429 
430  fPtVsYVsMultGenLargeAccFeeddw = new TH3F("hPtVsYVsMultGenLargeAccFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
432  fPtVsYVsMultGenLargeAccFeeddw->SetMinimum(0);
434 
435  fPtVsYVsMultGenLimAccFeeddw = new TH3F("hPtVsYVsMultGenLimAccFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
437  fPtVsYVsMultGenLimAccFeeddw->SetMinimum(0);
439 
440  fPtVsYVsMultGenAccFeeddw = new TH3F("hPtVsYVsMultGenAccFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
441  fPtVsYVsMultGenAccFeeddw->Sumw2();
442  fPtVsYVsMultGenAccFeeddw->SetMinimum(0);
444 
445  fPtVsYVsMultGenAccEvSelFeeddw = new TH3F("hPtVsYVsMultGenAccEvSelFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
447  fPtVsYVsMultGenAccEvSelFeeddw->SetMinimum(0);
449 
450  fPtVsYVsMultRecoFeeddw = new TH3F("hPtVsYVsMultRecoFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
451  fPtVsYVsMultRecoFeeddw->Sumw2();
452  fPtVsYVsMultRecoFeeddw->SetMinimum(0);
454 
455  }
456 
457 
458  Int_t nMassBins=static_cast<Int_t>(fMaxMass*1000.-fMinMass*1000.);
459  Double_t maxm=fMinMass+nMassBins*0.001;
460  fMassVsPtVsY=new TH3F("hMassVsPtVsY","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
461  fMassVsPtVsY->Sumw2();
462  fMassVsPtVsY->SetMinimum(0);
463  fOutput->Add(fMassVsPtVsY);
464 
465  fMassVsPtVsYRot=new TH3F("hMassVsPtVsYRot","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
466  fMassVsPtVsYRot->Sumw2();
467  fMassVsPtVsYRot->SetMinimum(0);
468  fOutput->Add(fMassVsPtVsYRot);
469 
470  fMassVsPtVsYLSpp=new TH3F("hMassVsPtVsYLSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
471  fMassVsPtVsYLSpp->Sumw2();
472  fMassVsPtVsYLSpp->SetMinimum(0);
474  fMassVsPtVsYLSmm=new TH3F("hMassVsPtVsYLSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
475  fMassVsPtVsYLSmm->Sumw2();
476  fMassVsPtVsYLSmm->SetMinimum(0);
478 
479  fMassVsPtVsYSig=new TH3F("hMassVsPtVsYSig","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
480  fMassVsPtVsYSig->Sumw2();
481  fMassVsPtVsYSig->SetMinimum(0);
482  fOutput->Add(fMassVsPtVsYSig);
483 
484  fMassVsPtVsYRefl=new TH3F("hMassVsPtVsYRefl","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
485  fMassVsPtVsYRefl->Sumw2();
486  fMassVsPtVsYRefl->SetMinimum(0);
488 
489  fMassVsPtVsYBkg=new TH3F("hMassVsPtVsYBkg","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
490  fMassVsPtVsYBkg->Sumw2();
491  fMassVsPtVsYBkg->SetMinimum(0);
492  fOutput->Add(fMassVsPtVsYBkg);
493 
494  fNSelected=new TH1F("hNSelected","",100,-0.5,99.5);
495  fNSelected->Sumw2();
496  fNSelected->SetMinimum(0);
497  fOutput->Add(fNSelected);
498 
499  fNormRotated=new TH1F("hNormRotated","",11,-0.5,10.5);
500  fNormRotated->Sumw2();
501  fNormRotated->SetMinimum(0);
502  fOutput->Add(fNormRotated);
503 
504  fDeltaMass=new TH1F("hDeltaMass","",100,-0.4,0.4);
505  fDeltaMass->Sumw2();
506  fDeltaMass->SetMinimum(0);
507  fOutput->Add(fDeltaMass);
508 
509  Int_t binSparseDMassRot[5]={nMassBins,100,24,40,20};
510  Double_t edgeLowSparseDMassRot[5]={fMinMass,-0.4,0.,-4.,0};
511  Double_t edgeHighSparseDMassRot[5]={maxm,0.4,12.,4.,3.14};
512  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);
514 
515  fMassVsPtVsYME=new TH3F("hMassVsPtVsYME","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
516  fMassVsPtVsYME->Sumw2();
517  fMassVsPtVsYME->SetMinimum(0);
518  fOutput->Add(fMassVsPtVsYME);
519 
520  fMassVsPtVsYMELSpp=new TH3F("hMassVsPtVsYMELSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
521  fMassVsPtVsYMELSpp->Sumw2();
522  fMassVsPtVsYMELSpp->SetMinimum(0);
524 
525  fMassVsPtVsYMELSmm=new TH3F("hMassVsPtVsYMELSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
526  fMassVsPtVsYMELSmm->Sumw2();
527  fMassVsPtVsYMELSmm->SetMinimum(0);
529 
532  if(fDoEventMixing==2) fNOfPools=1;
534  fEventsPerPool=new TH2F("hEventsPerPool","hEventsPerPool",fNzVertPools,fzVertPoolLims,fNMultPools,fMultPoolLims);
535  fMixingsPerPool=new TH2F("hMixingsPerPool","hMixingsPerPool",fNzVertPools,fzVertPoolLims,fNMultPools,fMultPoolLims);
536  }else{
537  fEventsPerPool=new TH2F("hEventsPerPool","hEventsPerPool",1,-10.,10.,1,-0.5,2000.5);
538  fMixingsPerPool=new TH2F("hMixingsPerPool","hMixingsPerPool",1,-10.,10.,1,-0.5,2000.5);
539  }
540  fEventsPerPool->Sumw2();
541  fEventsPerPool->SetMinimum(0);
542  fOutput->Add(fEventsPerPool);
543  fMixingsPerPool->Sumw2();
544  fMixingsPerPool->SetMinimum(0);
545  fOutput->Add(fMixingsPerPool);
546 
547  //Counter for Normalization
548  fCounter = new AliNormalizationCounter("NormalizationCounter");
549  fCounter->Init();
550 
551  fKaonTracks = new TObjArray();
552  fPionTracks=new TObjArray();
553  fKaonTracks->SetOwner();
554  fPionTracks->SetOwner();
555 
556  fEventBuffer = new TTree*[fNOfPools];
557  for(Int_t i=0; i<fNOfPools; i++){
558  fEventBuffer[i]=new TTree(Form("EventBuffer_%d",i), "Temporary buffer for event mixing");
559  fEventBuffer[i]->Branch("zVertex", &fVtxZ);
560  fEventBuffer[i]->Branch("multiplicity", &fMultiplicity);
561  fEventBuffer[i]->Branch("eventInfo", "TObjString",&fEventInfo);
562  fEventBuffer[i]->Branch("karray", "TObjArray", &fKaonTracks);
563  fEventBuffer[i]->Branch("parray", "TObjArray", &fPionTracks);
564  }
565 
566  PostData(1,fOutput);
567  PostData(2,fCounter);
568 }
569 
570 //________________________________________________________________________
573 
574  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
575  if(!aod && AODEvent() && IsStandardAOD()) {
576  // In case there is an AOD handler writing a standard AOD, use the AOD
577  // event in memory rather than the input (ESD) event.
578  aod = dynamic_cast<AliAODEvent*> (AODEvent());
579  }
580  if(!aod){
581  printf("AliAnalysisTaskCombinHF::UserExec: AOD not found!\n");
582  return;
583  }
584 
585  // fix for temporary bug in ESDfilter
586  // the AODs with null vertex pointer didn't pass the PhysSel
587  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
588 
589  // Reject events with trigger mask 0 of the LHC13d3 production
590  // For these events the ITS layers are skipped in the trakcing
591  // and the vertex reconstruction efficiency from tracks is biased
592  if(fReadMC){
593  Int_t runnumber = aod->GetRunNumber();
594  if(aod->GetTriggerMask()==0 &&
595  (runnumber>=195344 && runnumber<=195677)){
596  return;
597  }
598  }
599 
600  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
601  AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
602  AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
603  fPidHF->SetPidResponse(pidResp);
604 
605 
606  fHistNEvents->Fill(0); // count event
607  // Post the data already here
608  PostData(1,fOutput);
609 
611 
612  Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
616  }else{
618  fHistNEvents->Fill(8);
619  }else{
623  }else{
626  }
627  }
628  }
629 
631  // events not passing the centrality selection can be removed immediately. For the others we must count the generated D mesons
632 
633  Int_t ntracks=aod->GetNumberOfTracks();
634  fVtxZ = aod->GetPrimaryVertex()->GetZ();
636 
637  TClonesArray *arrayMC=0;
638  AliAODMCHeader *mcHeader=0;
639  if(fReadMC){
640  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
641  if(!arrayMC) {
642  printf("AliAnalysisTaskCombinHF::UserExec: MC particles branch not found!\n");
643  return;
644  }
645 
646  // load MC header
647  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
648  if(!mcHeader) {
649  printf("AliAnalysisTaskCombinHF::UserExec: MC header branch not found!\n");
650  return;
651  }
652  Double_t zMCVertex = mcHeader->GetVtxZ();
653  if (TMath::Abs(zMCVertex) < fAnalysisCuts->GetMaxVtxZ()){ // only cut on zVertex applied to count the signal
654  FillGenHistos(arrayMC,isEvSel);
655  }
656  fHistEventMultZv->Fill(zMCVertex,fMultiplicity);
657  if(isEvSel) fHistEventMultZvEvSel->Fill(zMCVertex,fMultiplicity);
658  }else{
660  if(isEvSel) fHistEventMultZvEvSel->Fill(fVtxZ,fMultiplicity);
661  }
662 
663 
664  if(!isEvSel)return;
665 
666  fHistNEvents->Fill(1);
667 
668 
669  // select and flag tracks
670  UChar_t* status = new UChar_t[ntracks];
671  for(Int_t iTr=0; iTr<ntracks; iTr++){
672  status[iTr]=0;
673  AliAODTrack* track=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr));
674  if(!track){
675  AliWarning("Error in casting track to AOD track. Not a standard AOD?");
676  continue;
677  }
678  if(IsTrackSelected(track)) status[iTr]+=1;
679 
680  // PID
681  if (fPIDstrategy == knSigma) {
682  // nsigma PID
683  if(IsKaon(track)) status[iTr]+=2;
684  if(IsPion(track)) status[iTr]+=4;
685  }
687  // Bayesian PID
688  Double_t *weights = new Double_t[AliPID::kSPECIES];
689  fPidHF->GetPidCombined()->ComputeProbabilities(track, fPidHF->GetPidResponse(), weights);
691  if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kKaon]) status[iTr] += 2;
692  if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kPion]) status[iTr] += 4;
693  }
694  if (fPIDstrategy == kBayesianThres) {
695  if (weights[AliPID::kKaon] > fBayesThresKaon) status[iTr] += 2;
696  if (weights[AliPID::kPion] > fBayesThresPion) status[iTr] += 4;
697  }
698  delete[] weights;
699  }
700 
701  fHistTrackStatus->Fill(status[iTr]);
702  fHistTrackEtaMultZv->Fill(track->Eta(),fVtxZ,fMultiplicity);
703  }
704 
705  // build the combinatorics
706  Int_t nSelected=0;
707  Int_t nFiltered=0;
708  Double_t dummypos[3]={0.,0.,0.};
709  AliAODVertex* v2=new AliAODVertex(dummypos,999.,-1,2);
710  AliAODVertex* v3=new AliAODVertex(dummypos,999.,-1,3);
711  // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
712  Double_t d02[2]={0.,0.};
713  Double_t d03[3]={0.,0.,0.};
714  AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
715  AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
716  UInt_t pdg0[2]={321,211};
717  UInt_t pdgp[3]={321,211,211};
718  UInt_t pdgs[3]={321,211,321};
719  Double_t tmpp[3];
720  Double_t px[3],py[3],pz[3];
721  Int_t dgLabels[3];
722  fKaonTracks->Delete();
723  fPionTracks->Delete();
724 
725  for(Int_t iTr1=0; iTr1<ntracks; iTr1++){
726  AliAODTrack* trK=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr1));
727  if(!trK){
728  AliWarning("Error in casting track to AOD track. Not a standard AOD?");
729  continue;
730  }
731  if((status[iTr1] & 1)==0) continue;
732  if(fDoEventMixing>0){
733  if(status[iTr1] & 2) fKaonTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
734  if(status[iTr1] & 4) fPionTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
735  }
736  if((status[iTr1] & 2)==0) continue;
737  Int_t chargeK=trK->Charge();
738  trK->GetPxPyPz(tmpp);
739  px[0] = tmpp[0];
740  py[0] = tmpp[1];
741  pz[0] = tmpp[2];
742  dgLabels[0]=trK->GetLabel();
743  for(Int_t iTr2=0; iTr2<ntracks; iTr2++){
744  if((status[iTr2] & 1)==0) continue;
745  if((status[iTr2] & 4)==0) continue;
746  if(iTr1==iTr2) continue;
747  AliAODTrack* trPi1=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr2));
748  if(!trPi1){
749  AliWarning("Error in casting track to AOD track. Not a standard AOD?");
750  continue;
751  }
752  Int_t chargePi1=trPi1->Charge();
753  trPi1->GetPxPyPz(tmpp);
754  px[1] = tmpp[0];
755  py[1] = tmpp[1];
756  pz[1] = tmpp[2];
757  dgLabels[1]=trPi1->GetLabel();
758  if(fMeson==kDzero){
759  if(chargePi1==chargeK){
760  // LS candidate
761  FillLSHistos(421,2,tmpRD2,px,py,pz,pdg0,chargePi1);
762  }else{
763  // OS candidate
764  nFiltered++;
765  v2->AddDaughter(trK);
766  v2->AddDaughter(trPi1);
767  tmpRD2->SetSecondaryVtx(v2);
768  Bool_t ok=FillHistos(421,2,tmpRD2,px,py,pz,pdg0,arrayMC,dgLabels);
769  v2->RemoveDaughters();
770  if(ok) nSelected++;
771  }
772  }else{
773  for(Int_t iTr3=iTr2+1; iTr3<ntracks; iTr3++){
774  if((status[iTr3] & 1)==0) continue;
775  if(fMeson==kDplus && (status[iTr3] & 4)==0) continue;
776  if(fMeson==kDs && (status[iTr3] & 2)==0) continue;
777  if(iTr1==iTr3) continue;
778  AliAODTrack* trPi2=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr3));
779  if(!trPi2){
780  AliWarning("Error in casting track to AOD track. Not a standard AOD?");
781  continue;
782  }
783  Int_t chargePi2=trPi2->Charge();
784  trPi2->GetPxPyPz(tmpp);
785  px[2] = tmpp[0];
786  py[2] = tmpp[1];
787  pz[2] = tmpp[2];
788  dgLabels[2]=trPi2->GetLabel();
789  if(fMeson==kDs){
790  Double_t massKK=ComputeInvMassKK(trK,trPi2);
791  Double_t deltaMass=massKK-TDatabasePDG::Instance()->GetParticle(333)->Mass();
792  if(TMath::Abs(deltaMass)>fPhiMassCut) continue;
793  tmpRD3->SetPxPyPzProngs(3,px,py,pz);
794  Double_t cos1=((AliAODRecoDecayHF3Prong*)tmpRD3)->CosPiKPhiRFrameKpiK();
795  Double_t kincutPiKPhi=TMath::Abs(cos1*cos1*cos1);
796  if(kincutPiKPhi<fCutCos3PiKPhiRFrame) continue;
797  Double_t cosPiDsLabFrame=((AliAODRecoDecayHF3Prong*)tmpRD3)->CosPiDsLabFrameKpiK();
798  if(cosPiDsLabFrame>fCutCosPiDsLabFrame) continue;
799  }
800  Bool_t isThreeLS=kFALSE;
801  if(chargePi1==chargeK && chargePi2==chargeK){
802  isThreeLS=kTRUE;
803  if(fMeson==kDplus)FillLSHistos(411,3,tmpRD3,px,py,pz,pdgp,chargePi1);
804  else if(fMeson==kDs)FillLSHistos(431,3,tmpRD3,px,py,pz,pdgs,chargePi1);
805  }
806  Bool_t acceptOS=kFALSE;
807  if(fMeson==kDplus){
808  if(chargePi1!=chargeK && chargePi2!=chargeK)acceptOS=kTRUE;
809  }else if(fMeson==kDs){
810  if(chargePi2!=chargeK && !isThreeLS) acceptOS=kTRUE;
811  }
812  if(acceptOS){
813  nFiltered++;
814  v3->AddDaughter(trK);
815  v3->AddDaughter(trPi1);
816  v3->AddDaughter(trPi2);
817  tmpRD3->SetSecondaryVtx(v3);
818  Bool_t ok=kFALSE;
819  if(fMeson==kDplus) ok=FillHistos(411,3,tmpRD3,px,py,pz,pdgp,arrayMC,dgLabels);
820  else if(fMeson==kDs) ok=FillHistos(431,3,tmpRD3,px,py,pz,pdgs,arrayMC,dgLabels);
821  v3->RemoveDaughters();
822  if(ok) nSelected++;
823  }
824  }
825  }
826  }
827  }
828 
829  delete [] status;
830  delete v2;
831  delete v3;
832  delete tmpRD2;
833  delete tmpRD3;
834 
835  fNSelected->Fill(nSelected);
836 
837  fCounter->StoreCandidates(aod,nFiltered,kTRUE);
838  fCounter->StoreCandidates(aod,nSelected,kFALSE);
839  fEventInfo->SetString(Form("Ev%d_esd%d_Pi%d_K%d",mgr->GetNcalls(),((AliAODHeader*)aod->GetHeader())->GetEventNumberESDFile(),fPionTracks->GetEntries(),fKaonTracks->GetEntries()));
840  if(fDoEventMixing==1){
842  if(ind>=0 && ind<fNOfPools){
844  fEventBuffer[ind]->Fill();
845  if(fEventBuffer[ind]->GetEntries() >= fNumberOfEventsForMixing){
847  DoMixingWithPools(ind);
848  ResetPool(ind);
849  }
850  }
851  }else if(fDoEventMixing==2){ // mix with cuts, no pools
852  fEventBuffer[0]->Fill();
853  }
854  PostData(1,fOutput);
855  PostData(2,fCounter);
856 
857  return;
858 }
859 
860 //________________________________________________________________________
861 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){
863 
864  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
865  Double_t pt = tmpRD->Pt();
866  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
867  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
868  Double_t rapid = tmpRD->Y(pdgD);
869  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
870  if(charge>0) fMassVsPtVsYLSpp->Fill(TMath::Sqrt(minv2),pt,rapid);
871  else fMassVsPtVsYLSmm->Fill(TMath::Sqrt(minv2),pt,rapid);
872  }
873  }
874  return;
875 }
876 
877 //________________________________________________________________________
878 void AliAnalysisTaskCombinHF::FillGenHistos(TClonesArray* arrayMC, Bool_t isEvSel){
880  Int_t totPart=arrayMC->GetEntriesFast();
881  Int_t thePDG=411;
882  Int_t nProng=3;
883  if(fMeson==kDzero){
884  thePDG=421;
885  nProng=2;
886  }else if(fMeson==kDs){
887  thePDG=431;
888  nProng=3;
889  }
890  for(Int_t ip=0; ip<totPart; ip++){
891  AliAODMCParticle *part = (AliAODMCParticle*)arrayMC->At(ip);
892  if(TMath::Abs(part->GetPdgCode())==thePDG){
894  fHistCheckOrigin->Fill(orig);
895  Int_t deca=0;
896  Bool_t isGoodDecay=kFALSE;
897  Int_t labDau[4]={-1,-1,-1,-1};
898  if(fMeson==kDzero){
899  deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,part,labDau);
900  if(part->GetNDaughters()!=2) continue;
901  if(deca==1) isGoodDecay=kTRUE;
902  }else if(fMeson==kDplus){
903  deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,part,labDau);
904  if(deca>0) isGoodDecay=kTRUE;
905  }else if(fMeson==kDs){
906  deca=AliVertexingHFUtils::CheckDsDecay(arrayMC,part,labDau);
907  if(deca==1) isGoodDecay=kTRUE;
908  }
909  fHistCheckDecChan->Fill(deca);
910  if(labDau[0]==-1){
911  // printf(Form("Meson %d Label of daughters not filled correctly -- %d\n",fMeson,isGoodDecay));
912  continue; //protection against unfilled array of labels
913  }
914  Bool_t isInAcc=CheckAcceptance(arrayMC,nProng,labDau);
915  if(isInAcc) fHistCheckDecChanAcc->Fill(deca);
916  if(isGoodDecay){
917  Double_t ptgen=part->Pt();
918  Double_t ygen=part->Y();
919  if(fAnalysisCuts->IsInFiducialAcceptance(ptgen,ygen)){
920  if(orig==4){
921  fPtVsYVsMultGenPrompt->Fill(ptgen,ygen,fMultiplicity);
922  if(TMath::Abs(ygen)<0.5) fPtVsYVsMultGenLimAccPrompt->Fill(ptgen,ygen,fMultiplicity);
923  if(isInAcc) fPtVsYVsMultGenAccPrompt->Fill(ptgen,ygen,fMultiplicity);
924  if(isEvSel && isInAcc) fPtVsYVsMultGenAccEvSelPrompt->Fill(ptgen,ygen,fMultiplicity);
925  }else if(orig==5){
926  fPtVsYVsMultGenFeeddw->Fill(ptgen,ygen,fMultiplicity);
927  if(TMath::Abs(ygen)<0.5) fPtVsYVsMultGenLimAccFeeddw->Fill(ptgen,ygen,fMultiplicity);
928  if(isInAcc) fPtVsYVsMultGenAccFeeddw->Fill(ptgen,ygen,fMultiplicity);
929  if(isEvSel && isInAcc) fPtVsYVsMultGenAccEvSelFeeddw->Fill(ptgen,ygen,fMultiplicity);
930  }
931  }
932  if(TMath::Abs(ygen)<0.9){
933  if(orig==4) fPtVsYVsMultGenLargeAccPrompt->Fill(ptgen,ygen,fMultiplicity);
934  else if(orig==5) fPtVsYVsMultGenLargeAccFeeddw->Fill(ptgen,ygen,fMultiplicity);
935  }
936  }
937  }
938  }
939 }
940 
941 //________________________________________________________________________
942 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, Int_t* dgLabels){
944 
945  Bool_t accept=kFALSE;
946 
947  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
948  Double_t pt = tmpRD->Pt();
949  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
950  Double_t mass=TMath::Sqrt(minv2);
951 
952  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
953  Double_t rapid = tmpRD->Y(pdgD);
954  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
955  fMassVsPtVsY->Fill(mass,pt,rapid);
956  accept=kTRUE;
957  if(fReadMC){
958  Int_t signPdg[3]={0,0,0};
959  for(Int_t iii=0; iii<nProngs; iii++) signPdg[iii]=pdgdau[iii];
960  Int_t labD = tmpRD->MatchToMC(pdgD,arrayMC,nProngs,signPdg);
961  if(labD>=0){
962  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(TMath::Abs(dgLabels[0])));
963  if(part){
965  Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
966  if(pdgCode==321){
967  fMassVsPtVsYSig->Fill(mass,pt,rapid);
968  AliAODMCParticle* dmes = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labD));
969  if(dmes){
970  if(orig==4) fPtVsYVsMultRecoPrompt->Fill(dmes->Pt(),dmes->Y(),fMultiplicity);
971  else if(orig==5) fPtVsYVsMultRecoFeeddw->Fill(dmes->Pt(),dmes->Y(),fMultiplicity);
972  }
973  }else{
974  fMassVsPtVsYRefl->Fill(mass,pt,rapid);
975  }
976  }
977  }else{
978  fMassVsPtVsYBkg->Fill(mass,pt,rapid);
979  }
980  }
981  }
982  }
983 
984  Int_t nRotated=0;
985  Double_t massRot=0;// calculated later only if candidate is acceptable
986  Double_t angleProngXY;
987  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])));
988  else if(TMath::Abs(pdgD)==431) {
989  Double_t px_phi = px[0]+px[2];
990  Double_t py_phi = py[0]+py[2];
991  Double_t pz_phi = pz[0]+pz[2];
992  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])));
993  }
994  else {//angle between pion and phi meson
995  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])));
996  }
997  Double_t ptOrig=pt;
998 
999 
1000  Double_t rotStep=0.;
1001  if(fNRotations>1) rotStep=(fMaxAngleForRot-fMinAngleForRot)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
1002  if(TMath::Abs(pdgD)==421 || TMath::Abs(pdgD)==431) fNRotations3=1;
1003  Double_t rotStep3=0.;
1004  if(fNRotations3>1) rotStep3=(fMaxAngleForRot3-fMinAngleForRot3)/(fNRotations3-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
1005 
1006  for(Int_t irot=0; irot<fNRotations; irot++){
1007  Double_t phirot=fMinAngleForRot+rotStep*irot;
1008  Double_t tmpx=px[0];
1009  Double_t tmpy=py[0];
1010  Double_t tmpx2=px[2];
1011  Double_t tmpy2=py[2];
1012  if(pdgD==431) {
1013  //rotate pion w.r.t. phi meson
1014  tmpx=px[1];
1015  tmpy=py[1];
1016  px[1]=tmpx*TMath::Cos(phirot)-tmpy*TMath::Sin(phirot);
1017  py[1]=tmpx*TMath::Sin(phirot)+tmpy*TMath::Cos(phirot);
1018  }
1019  else {
1020  px[0]=tmpx*TMath::Cos(phirot)-tmpy*TMath::Sin(phirot);
1021  py[0]=tmpx*TMath::Sin(phirot)+tmpy*TMath::Cos(phirot);
1022  }
1023  for(Int_t irot3=0; irot3<fNRotations3; irot3++){
1024  if(pdgD==411){
1025  Double_t phirot2=fMaxAngleForRot3-rotStep3*irot;
1026  px[2]=tmpx*TMath::Cos(phirot2)-tmpy*TMath::Sin(phirot2);
1027  py[2]=tmpx*TMath::Sin(phirot2)+tmpy*TMath::Cos(phirot2);
1028  }
1029  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
1030  pt = tmpRD->Pt();
1031  minv2 = tmpRD->InvMass2(nProngs,pdgdau);
1032  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
1033  Double_t rapid = tmpRD->Y(pdgD);
1034  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
1035  massRot=TMath::Sqrt(minv2);
1036  fMassVsPtVsYRot->Fill(massRot,pt,rapid);
1037  nRotated++;
1038  fDeltaMass->Fill(massRot-mass);
1039  if(fFullAnalysis){
1040  Double_t pointRot[5]={mass,massRot-mass,ptOrig,pt-ptOrig,angleProngXY};
1041  fDeltaMassFullAnalysis->Fill(pointRot);
1042  }
1043  }
1044  }
1045  }
1046  if(pdgD==431) {
1047  px[1]=tmpx;
1048  py[1]=tmpy;
1049  }
1050  else {
1051  px[0]=tmpx;
1052  py[0]=tmpy;
1053  if(pdgD==411){
1054  px[2]=tmpx2;
1055  py[2]=tmpy2;
1056  }
1057  }
1058  }
1059  fNormRotated->Fill(nRotated);
1060 
1061  return accept;
1062 
1063 }
1064 //________________________________________________________________________
1065 void AliAnalysisTaskCombinHF::FillMEHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau){
1067 
1068  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
1069  Double_t pt = tmpRD->Pt();
1070  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
1071  Double_t mass=TMath::Sqrt(minv2);
1072 
1073  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
1074  Double_t rapid = tmpRD->Y(pdgD);
1075  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
1076  fMassVsPtVsYME->Fill(mass,pt,rapid);
1077  }
1078  }
1079  return;
1080 }
1081 //________________________________________________________________________
1082 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){
1084 
1085  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
1086  Double_t pt = tmpRD->Pt();
1087  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
1088  Double_t mass=TMath::Sqrt(minv2);
1089 
1090  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
1091  Double_t rapid = tmpRD->Y(pdgD);
1092  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
1093  if(charge>0) fMassVsPtVsYMELSpp->Fill(mass,pt,rapid);
1094  if(charge<0) fMassVsPtVsYMELSmm->Fill(mass,pt,rapid);
1095  }
1096  }
1097  return;
1098 }
1099 //________________________________________________________________________
1102 
1103  if(track->Charge()==0) return kFALSE;
1104  if(track->GetID()<0&&!fKeepNegID)return kFALSE;
1105  if(fFilterMask>=0){
1106  if(!(track->TestFilterMask(fFilterMask))) return kFALSE;
1107  }
1108  if(!SelectAODTrack(track,fTrackCutsAll)) return kFALSE;
1109  return kTRUE;
1110 }
1111 
1112 //________________________________________________________________________
1115 
1116  if(!fPidHF) return kTRUE;
1117  Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
1118  Double_t mom=track->P();
1119  if(SelectAODTrack(track,fTrackCutsKaon)) {
1120  if(isKaon>=1) return kTRUE;
1121  if(isKaon<=-1) return kFALSE;
1122  switch(fPIDselCaseZero){// isKaon=0
1123  case 0:
1124  {
1125  return kTRUE;// always accept
1126  }
1127  break;
1128  case 1:
1129  {
1130  if(isKaon>=0 && track->P()>fmaxPforIDKaon)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDKaon
1131  }
1132  break;
1133  case 2:
1134  {
1135  if(track->P()>fmaxPforIDKaon)return kTRUE;
1136  AliPIDResponse *pidResp=fPidHF->GetPidResponse();// more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
1137  Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kKaon);
1138  if(nsigma>-2.&& nsigma<3. && mom<0.6)isKaon=1;
1139  else if(nsigma>-1.&& nsigma<3.&& mom<0.8)isKaon=1;
1140  if(isKaon==1)return kTRUE;
1141  }
1142  break;
1143  default:
1144  {
1145  AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
1146  return kFALSE;// actually case 0 could be set as the default and return kTRUE
1147  }
1148  }
1149  }
1150 
1151  return kFALSE;
1152 }
1153 //_______________________________________________________________________
1156 
1157  if(!fPidHF) return kTRUE;
1158  Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
1159  Double_t mom=track->P();
1160  if(SelectAODTrack(track,fTrackCutsPion)) {
1161  if(isPion>=1) return kTRUE;
1162  if(isPion<=-1) return kFALSE;
1163  switch(fPIDselCaseZero){// isPion=0
1164  case 0:
1165  {
1166  return kTRUE;// always accept
1167  }
1168  break;
1169  case 1:
1170  {
1171  if(track->P()>fmaxPforIDPion)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDPion
1172  }
1173  break;
1174  case 2:
1175  {
1176  // more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
1177  if(track->P()>fmaxPforIDPion)return kTRUE;
1178  AliPIDResponse *pidResp=fPidHF->GetPidResponse();
1179  Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kPion);
1180  if(nsigma<2.&& nsigma>-3. && mom<0.6)isPion=1;
1181  else if(nsigma<1. && nsigma> -3. && mom<0.8)isPion=1;
1182  if(isPion==1)return kTRUE;
1183  }
1184  break;
1185  default:
1186  {
1187  AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
1188  return kFALSE;// actually case 0 could be set as the default and return kTRUE
1189  }
1190  }
1191  }
1192 
1193  return kFALSE;
1194 }
1195 
1196 //________________________________________________________________________
1197 Bool_t AliAnalysisTaskCombinHF::SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts){
1199 
1200  if(!cuts) return kTRUE;
1201 
1202  AliESDtrack esdTrack(track);
1203  // set the TPC cluster info
1204  esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
1205  esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
1206  esdTrack.SetTPCPointsF(track->GetTPCNclsF());
1207  if(!cuts->IsSelected(&esdTrack)) return kFALSE;
1208 
1209  return kTRUE;
1210 }
1211 
1212 //_________________________________________________________________
1213 Bool_t AliAnalysisTaskCombinHF::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1215  for (Int_t iProng = 0; iProng<nProng; iProng++){
1216  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1217  if(!mcPartDaughter) return kFALSE;
1218  Double_t eta = mcPartDaughter->Eta();
1219  Double_t pt = mcPartDaughter->Pt();
1220  if (TMath::Abs(eta) > fEtaAccCut || pt < fPtAccCut) return kFALSE;
1221  }
1222  return kTRUE;
1223 }
1224 //_________________________________________________________________
1227  if(!fzVertPoolLims || !fMultPoolLims) return 0;
1228  Int_t theBinZ=TMath::BinarySearch(fNzVertPoolsLimSize,fzVertPoolLims,zvert);
1229  if(theBinZ<0 || theBinZ>=fNzVertPoolsLimSize) return -1;
1230  Int_t theBinM=TMath::BinarySearch(fNMultPoolsLimSize,fMultPoolLims,mult);
1231  if(theBinM<0 || theBinM>=fNMultPoolsLimSize) return -1;
1232  return fNMultPools*theBinZ+theBinM;
1233 }
1234 //_________________________________________________________________
1237  if(poolIndex<0 || poolIndex>=fNOfPools) return;
1238  delete fEventBuffer[poolIndex];
1239  fEventBuffer[poolIndex]=new TTree(Form("EventBuffer_%d",poolIndex), "Temporary buffer for event mixing");
1240  fEventBuffer[poolIndex]->Branch("zVertex", &fVtxZ);
1241  fEventBuffer[poolIndex]->Branch("multiplicity", &fMultiplicity);
1242  fEventBuffer[poolIndex]->Branch("eventInfo", "TObjString",&fEventInfo);
1243  fEventBuffer[poolIndex]->Branch("karray", "TObjArray", &fKaonTracks);
1244  fEventBuffer[poolIndex]->Branch("parray", "TObjArray", &fPionTracks);
1245  return;
1246 }
1247 //_________________________________________________________________
1250  if(TMath::Abs(zv2-zv1)>fMaxzVertDistForMix) return kFALSE;
1251  if(TMath::Abs(mult2-mult1)>fMaxMultDiffForMix) return kFALSE;
1252  return kTRUE;
1253 }
1254 //_________________________________________________________________
1257 
1258  if(fDoEventMixing==0) return;
1259  Int_t nEvents=fEventBuffer[0]->GetEntries();
1260  if(fDebug > 1) printf("AnalysisTaskCombinHF::DoMixingWithCuts Start Event Mixing of %d events\n",nEvents);
1261 
1262  TObjArray* karray=0x0;
1263  TObjArray* parray=0x0;
1264  Double_t zVertex,mult;
1265  TObjString* eventInfo=0x0;
1266  fEventBuffer[0]->SetBranchAddress("karray", &karray);
1267  fEventBuffer[0]->SetBranchAddress("parray", &parray);
1268  fEventBuffer[0]->SetBranchAddress("eventInfo",&eventInfo);
1269  fEventBuffer[0]->SetBranchAddress("zVertex", &zVertex);
1270  fEventBuffer[0]->SetBranchAddress("multiplicity", &mult);
1271  Double_t d02[2]={0.,0.};
1272  Double_t d03[3]={0.,0.,0.};
1273  AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
1274  AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
1275  UInt_t pdg0[2]={321,211};
1276  // UInt_t pdgp[3]={321,211,211};
1277  Double_t px[3],py[3],pz[3];
1278  Int_t evId1,esdId1,nk1,np1;
1279  Int_t evId2,esdId2,nk2,np2;
1280 
1281  for(Int_t iEv1=0; iEv1<nEvents; iEv1++){
1282  fEventBuffer[0]->GetEvent(iEv1);
1283  TObjArray* karray1=(TObjArray*)karray->Clone();
1284  Double_t zVertex1=zVertex;
1285  Double_t mult1=mult;
1286  Int_t nKaons=karray1->GetEntries();
1287  Int_t nPionsForCheck=parray->GetEntries();
1288  sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId1,&esdId1,&np1,&nk1);
1289  if(nk1!=nKaons || np1!=nPionsForCheck){
1290  printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: read event does not match to the stored one\n");
1291  delete karray1;
1292  continue;
1293  }
1294  for(Int_t iEv2=0; iEv2<fNumberOfEventsForMixing; iEv2++){
1295  Int_t iToMix=iEv1+iEv2+1;
1296  if(iEv1>=(nEvents-fNumberOfEventsForMixing)) iToMix=iEv1-iEv2-1;
1297  if(iToMix<0) continue;
1298  if(iToMix==iEv1) continue;
1299  if(iToMix<iEv1) continue;
1300  fEventBuffer[0]->GetEvent(iToMix);
1301  Double_t zVertex2=zVertex;
1302  Double_t mult2=mult;
1303  if(TMath::Abs(zVertex2-zVertex1)<0.0001 && TMath::Abs(mult2-mult1)<0.001){
1304  printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: same event in mixing??? %d %d %f %f %f %f\n",iEv1,iEv2,zVertex1,zVertex2,mult1,mult2);
1305  continue;
1306  }
1307  TObjArray* parray2=(TObjArray*)parray->Clone();
1308  Int_t nPions=parray2->GetEntries();
1309  Int_t nKaonsForCheck=karray->GetEntries();
1310  sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId2,&esdId2,&np2,&nk2);
1311  if(nk2!=nKaonsForCheck || np2!=nPions){
1312  printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: read event does not match to the stored one\n");
1313  delete parray2;
1314  continue;
1315  }
1316  if(evId2==evId1 && esdId2==esdId1){
1317  printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: same event in mixing??? %d %d nK=%d %d nPi=%d %d\n",evId1,evId2,nKaons,nKaonsForCheck,nPionsForCheck,nPions);
1318  delete parray2;
1319  continue;
1320  }
1321  if(CanBeMixed(zVertex1,zVertex2,mult1,mult2)){
1322  for(Int_t iTr1=0; iTr1<nKaons; iTr1++){
1323  TLorentzVector* trK=(TLorentzVector*)karray1->At(iTr1);
1324  Double_t chargeK=trK->T();
1325  px[0] = trK->Px();
1326  py[0] = trK->Py();
1327  pz[0] = trK->Pz();
1328  for(Int_t iTr2=0; iTr2<nPions; iTr2++){
1329  TLorentzVector* trPi1=(TLorentzVector*)parray2->At(iTr2);
1330  Double_t chargePi1=trPi1->T();
1331  px[1] = trPi1->Px();
1332  py[1] = trPi1->Py();
1333  pz[1] = trPi1->Pz();
1334  if(fMeson==kDzero && chargePi1*chargeK<0){
1335  FillMEHistos(421,2,tmpRD2,px,py,pz,pdg0);
1336  }
1337  if(fMeson==kDzero && chargePi1*chargeK>0){
1338  FillMEHistosLS(421,2,tmpRD2,px,py,pz,pdg0,(Int_t)chargePi1);
1339  }
1340  }
1341  }
1342  }
1343  delete parray2;
1344  }
1345  delete karray1;
1346  }
1347  delete tmpRD2;
1348  delete tmpRD3;
1349 }
1350 //_________________________________________________________________
1353 
1354  if(fDoEventMixing==0) return;
1355  if(poolIndex<0 || poolIndex>fNzVertPools*fNMultPools) return;
1356 
1357  Int_t nEvents=fEventBuffer[poolIndex]->GetEntries();
1358  if(fDebug > 1) printf("AliAnalysisTaskCombinHF::DoMixingWithPools Start Event Mixing of %d events\n",nEvents);
1359  TObjArray* karray=0x0;
1360  TObjArray* parray=0x0;
1361  Double_t zVertex,mult;
1362  TObjString* eventInfo=0x0;
1363  fEventBuffer[poolIndex]->SetBranchAddress("karray", &karray);
1364  fEventBuffer[poolIndex]->SetBranchAddress("parray", &parray);
1365  fEventBuffer[poolIndex]->SetBranchAddress("eventInfo",&eventInfo);
1366  fEventBuffer[poolIndex]->SetBranchAddress("zVertex", &zVertex);
1367  fEventBuffer[poolIndex]->SetBranchAddress("multiplicity", &mult);
1368 
1369  // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
1370  Double_t d02[2]={0.,0.};
1371  Double_t d03[3]={0.,0.,0.};
1372  AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
1373  AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
1374  UInt_t pdg0[2]={321,211};
1375  UInt_t pdgp[3]={321,211,211};
1376  UInt_t pdgs[3]={321,211,321};
1377  Double_t px[3],py[3],pz[3];
1378  Int_t evId1,esdId1,nk1,np1;
1379  Int_t evId2,esdId2,nk2,np2;
1380 
1381  for(Int_t iEv1=0; iEv1<nEvents; iEv1++){
1382  fEventBuffer[poolIndex]->GetEvent(iEv1);
1383  TObjArray* karray1=(TObjArray*)karray->Clone();
1384  Double_t zVertex1=zVertex;
1385  Double_t mult1=mult;
1386  Int_t nKaons=karray1->GetEntries();
1387  Int_t nPionsForCheck=parray->GetEntries();
1388  sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId1,&esdId1,&np1,&nk1);
1389  if(nk1!=nKaons || np1!=nPionsForCheck){
1390  printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: read event does not match to the stored one\n");
1391  delete karray1;
1392  continue;
1393  }
1394  for(Int_t iEv2=0; iEv2<nEvents; iEv2++){
1395  if(iEv2==iEv1) continue;
1396  fEventBuffer[poolIndex]->GetEvent(iEv2);
1397  Double_t zVertex2=zVertex;
1398  Double_t mult2=mult;
1399  if(TMath::Abs(zVertex2-zVertex1)<0.0001 && TMath::Abs(mult2-mult1)<0.001){
1400  printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: same event in mixing??? %d %d %f %f %f %f\n",iEv1,iEv2,zVertex1,zVertex2,mult1,mult2);
1401  continue;
1402  }
1403  TObjArray* parray2=(TObjArray*)parray->Clone();
1404  Int_t nPions=parray2->GetEntries();
1405  Int_t nKaonsForCheck=karray->GetEntries();
1406  sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId2,&esdId2,&np2,&nk2);
1407  if(nk2!=nKaonsForCheck || np2!=nPions){
1408  printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: read event does not match to the stored one\n");
1409  delete parray2;
1410  continue;
1411  }
1412  if(evId2==evId1 && esdId2==esdId1){
1413  printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: same event in mixing??? %d %d nK=%d %d nPi=%d %d\n",evId1,evId2,nKaons,nKaonsForCheck,nPionsForCheck,nPions);
1414  delete parray2;
1415  continue;
1416  }
1417  TObjArray* parray3=0x0;
1418  Int_t nPions3=0;
1419  if(fMeson!=kDzero && fMeson!=kDs){
1420  Int_t iEv3=iEv2+1;
1421  if(iEv3==iEv1) iEv3=iEv2+2;
1422  if(iEv3>=nEvents) iEv3=iEv2-3;
1423  if(nEvents==2) iEv3=iEv1;
1424  if(iEv3<0) iEv3=iEv2-1;
1425  fEventBuffer[poolIndex]->GetEvent(iEv3);
1426  parray3=(TObjArray*)parray->Clone();
1427  nPions3=parray3->GetEntries();
1428  }
1429  for(Int_t iTr1=0; iTr1<nKaons; iTr1++){
1430  TLorentzVector* trK=(TLorentzVector*)karray1->At(iTr1);
1431  Double_t chargeK=trK->T();
1432  px[0] = trK->Px();
1433  py[0] = trK->Py();
1434  pz[0] = trK->Pz();
1435  for(Int_t iTr2=0; iTr2<nPions; iTr2++){
1436  TLorentzVector* trPi1=(TLorentzVector*)parray2->At(iTr2);
1437  Double_t chargePi1=trPi1->T();
1438  px[1] = trPi1->Px();
1439  py[1] = trPi1->Py();
1440  pz[1] = trPi1->Pz();
1441  if(chargePi1*chargeK<0){
1442  if(fMeson==kDzero){
1443  FillMEHistos(421,2,tmpRD2,px,py,pz,pdg0);
1444  }else if(fMeson==kDs) {
1445  for(Int_t iTr3=iTr1+1; iTr3<nKaons; iTr3++){
1446  TLorentzVector* trK2=(TLorentzVector*)karray1->At(iTr3);
1447  Double_t chargeK2=trK2->T();
1448  px[2] = trK2->Px();
1449  py[2] = trK2->Py();
1450  pz[2] = trK2->Pz();
1451  Double_t massKK=ComputeInvMassKK(trK,trK2);
1452  Double_t deltaMass=massKK-TDatabasePDG::Instance()->GetParticle(333)->Mass();
1453  Double_t cos1=CosPiKPhiRFrame(trK,trK2,trPi1);
1454  Double_t kincutPiKPhi=TMath::Abs(cos1*cos1*cos1);
1455  Double_t cosPiDsLabFrame=CosPiDsLabFrame(trK,trK2,trPi1);
1456  if(chargeK2*chargeK<0 && TMath::Abs(deltaMass)<fPhiMassCut && kincutPiKPhi>fCutCos3PiKPhiRFrame && cosPiDsLabFrame<fCutCosPiDsLabFrame){
1457  FillMEHistos(431,3,tmpRD3,px,py,pz,pdgs);
1458  }
1459  }
1460  } else{
1461  if(parray3){
1462  for(Int_t iTr3=iTr2+1; iTr3<nPions3; iTr3++){
1463  TLorentzVector* trPi2=(TLorentzVector*)parray3->At(iTr3);
1464  Double_t chargePi2=trPi2->T();
1465  px[2] = trPi2->Px();
1466  py[2] = trPi2->Py();
1467  pz[2] = trPi2->Pz();
1468  if(chargePi2*chargeK<0){
1469  if(fMeson==kDplus) FillMEHistos(411,3,tmpRD3,px,py,pz,pdgp);
1470  }
1471  }
1472  }
1473  }
1474  }else if(chargePi1*chargeK>0){
1475  if(fMeson==kDzero) FillMEHistosLS(421,2,tmpRD2,px,py,pz,pdg0,(Int_t)chargePi1);
1476  }
1477  }
1478  }
1479  delete parray3;
1480  delete parray2;
1481  }
1482  delete karray1;
1483  }
1484  delete tmpRD2;
1485  delete tmpRD3;
1486 }
1487 //_________________________________________________________________
1489 {
1491  if(fDoEventMixing==0) return;
1492  printf("AliAnalysisTaskCombinHF: FinishTaskOutput\n");
1493 
1494  if(fDoEventMixing==1){
1495  for(Int_t i=0; i<fNOfPools; i++){
1496  Int_t nEvents=fEventBuffer[i]->GetEntries();
1497  if(nEvents>1) DoMixingWithPools(i);
1498  }
1499  }else if(fDoEventMixing==2){
1500  DoMixingWithCuts();
1501  }
1502 }
1503 //_________________________________________________________________
1504 Double_t AliAnalysisTaskCombinHF::ComputeInvMassKK(AliAODTrack* tr1, AliAODTrack* tr2) const{
1506  Double_t massK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
1507  Double_t p1=tr1->P();
1508  Double_t p2=tr2->P();
1509  Double_t pxtot=tr1->Px()+tr2->Px();
1510  Double_t pytot=tr1->Py()+tr2->Py();
1511  Double_t pztot=tr1->Pz()+tr2->Pz();
1512  Double_t e1=TMath::Sqrt(massK*massK+p1*p1);
1513  Double_t e2=TMath::Sqrt(massK*massK+p2*p2);
1514  Double_t etot=e1+e2;
1515  Double_t m2=etot*etot-(pxtot*pxtot+pytot*pytot+pztot*pztot);
1516  return TMath::Sqrt(m2);
1517 }
1518 //_________________________________________________________________
1519 Double_t AliAnalysisTaskCombinHF::ComputeInvMassKK(TLorentzVector* tr1, TLorentzVector* tr2) const{
1521  Double_t massK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
1522  Double_t p1=tr1->P();
1523  Double_t p2=tr2->P();
1524  Double_t pxtot=tr1->Px()+tr2->Px();
1525  Double_t pytot=tr1->Py()+tr2->Py();
1526  Double_t pztot=tr1->Pz()+tr2->Pz();
1527  Double_t e1=TMath::Sqrt(massK*massK+p1*p1);
1528  Double_t e2=TMath::Sqrt(massK*massK+p2*p2);
1529  Double_t etot=e1+e2;
1530  Double_t m2=etot*etot-(pxtot*pxtot+pytot*pytot+pztot*pztot);
1531  return TMath::Sqrt(m2);
1532 }
1533 
1534 //----------------------------------------------------------------------
1535 Double_t AliAnalysisTaskCombinHF::CosPiKPhiRFrame(TLorentzVector* dauK1, TLorentzVector* dauK2, TLorentzVector* daupi) const {
1537 
1538  Double_t massK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
1539  Double_t masspi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
1540 
1541  Double_t eK1 = TMath::Sqrt(massK*massK+dauK1->P()*dauK1->P());
1542  Double_t eK2 = TMath::Sqrt(massK*massK+dauK2->P()*dauK2->P());
1543  Double_t epi = TMath::Sqrt(masspi*masspi+daupi->P()*daupi->P());
1544 
1545  Double_t ePhi=eK1+eK2;
1546  Double_t pxPhi=dauK1->Px()+dauK2->Px();
1547  Double_t pyPhi=dauK1->Py()+dauK2->Py();
1548  Double_t pzPhi=dauK1->Pz()+dauK2->Pz();
1549  Double_t bxPhi=pxPhi/ePhi;
1550  Double_t byPhi=pyPhi/ePhi;
1551  Double_t bzPhi=pzPhi/ePhi;
1552 
1553  TVector3 vecK1Phiframe;
1554  TLorentzVector vecK1(dauK1->Px(),dauK1->Py(),dauK1->Pz(),eK1);
1555  vecK1.Boost(-bxPhi,-byPhi,-bzPhi);
1556  vecK1.Boost(vecK1Phiframe);
1557  vecK1Phiframe=vecK1.BoostVector();
1558 
1559  TVector3 vecPiPhiframe;
1560  TLorentzVector vecPi(daupi->Px(),daupi->Py(),daupi->Pz(),epi);
1561  vecPi.Boost(-bxPhi,-byPhi,-bzPhi);
1562  vecPi.Boost(vecPiPhiframe);
1563  vecPiPhiframe=vecPi.BoostVector();
1564 
1565  Double_t innera=vecPiPhiframe.Dot(vecK1Phiframe);
1566  Double_t norm1a=TMath::Sqrt(vecPiPhiframe.Dot(vecPiPhiframe));
1567  Double_t norm2a=TMath::Sqrt(vecK1Phiframe.Dot(vecK1Phiframe));
1568  Double_t cosK1PhiFrame=innera/(norm1a*norm2a);
1569 
1570  return cosK1PhiFrame;
1571 }
1572 
1573 //----------------------------------------------------------------------
1574 Double_t AliAnalysisTaskCombinHF::CosPiDsLabFrame(TLorentzVector* dauK1, TLorentzVector* dauK2, TLorentzVector* daupi) const {
1576 
1577  Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1578  Double_t masspi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
1579 
1580  Double_t pxD=dauK1->Px()+dauK2->Px()+daupi->Px();
1581  Double_t pyD=dauK1->Px()+dauK2->Px()+daupi->Px();
1582  Double_t pzD=dauK1->Px()+dauK2->Px()+daupi->Px();
1583  Double_t pD=TMath::Sqrt(pxD*pxD+pyD*pyD+pzD*pzD);
1584  Double_t eD=TMath::Sqrt(massDs*massDs+pD*pD);
1585  Double_t epi = TMath::Sqrt(masspi*masspi+daupi->P()*daupi->P());
1586  Double_t bxD=pxD/eD;
1587  Double_t byD=pyD/eD;
1588  Double_t bzD=pzD/eD;
1589 
1590  TVector3 piDsframe;
1591  TLorentzVector vecPi(daupi->Px(),daupi->Py(),daupi->Pz(),epi);
1592  vecPi.Boost(-bxD,-byD,-bzD);
1593  vecPi.Boost(piDsframe);
1594  piDsframe=vecPi.BoostVector();
1595 
1596  TVector3 vecDs(pxD,pyD,pzD);
1597 
1598  Double_t inner=vecDs.Dot(piDsframe);
1599  Double_t norm1=TMath::Sqrt(vecDs.Dot(vecDs));
1600  Double_t norm2=TMath::Sqrt(piDsframe.Dot(piDsframe));
1601  Double_t cosPiDsFrame=inner/(norm1*norm2);
1602 
1603  return cosPiDsFrame;
1604 }
1605 
1606 //_________________________________________________________________
1608 {
1610  //
1611  if(fDebug > 1) printf("AliAnalysisTaskCombinHF: Terminate() \n");
1612  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1613  if (!fOutput) {
1614  printf("ERROR: fOutput not available\n");
1615  return;
1616  }
1617  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
1618  if(fHistNEvents){
1619  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1620  }else{
1621  printf("ERROR: fHistNEvents not available\n");
1622  return;
1623  }
1624  return;
1625 }
1626 
Int_t charge
Bool_t IsEventRejectedDueToCentrality() const
Definition: AliRDHFCuts.h:330
Double_t fEtaAccCut
width of pt bin (GeV/c)
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, Int_t *dgLabels)
Bool_t IsEventRejectedDueToZVertexOutsideFiducialRegion() const
Definition: AliRDHFCuts.h:324
Bool_t IsTrackSelected(AliAODTrack *track)
static Int_t CheckD0Decay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
Bool_t IsEventRejectedDueToNotRecoVertex() const
Definition: AliRDHFCuts.h:318
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
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 send on output slot 0
Bool_t fReadMC
mesonSpecies (see enum)
TH3F * fMassVsPtVsY
! hist. of Y vs. Pt vs. Mass (all cand)
TH2F * fHistEventMultZvEvSel
!hist. of evnt Mult vs. Zv for selected ev
Double_t fMinAngleForRot
number of rotations
Double_t mass
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
TH3F * fPtVsYVsMultRecoPrompt
! hist. of Y vs. Pt vs. Mult generated (Reco D)
Bool_t IsEventRejectedDueToVertexContributors() const
Definition: AliRDHFCuts.h:321
TH3F * fPtVsYVsMultGenLargeAccFeeddw
! hist. of Y vs. Pt vs. Mult generated (|y|<0.9)
Double_t fPtBinWidth
maximum pT value for inv. mass histograms
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
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
void FillGenHistos(TClonesArray *arrayMC, Bool_t isEvSel)
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:256
Double_t CosPiDsLabFrame(TLorentzVector *dauK1, TLorentzVector *dauK2, TLorentzVector *daupi) const
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
TH2F * fHistEventMultZv
!hist. of evnt Mult vs. Zv for all events
Bool_t CanBeMixed(Double_t zv1, Double_t zv2, Double_t mult1, Double_t mult2)
Double_t fPhiMassCut
kaon track selection
const Int_t nPtBins
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 * fMassVsPtVsYLSmm
! hist. of Y vs. Pt vs. Mass (like sign –)
int Int_t
Definition: External.C:63
Double_t CosPiKPhiRFrame(TLorentzVector *dauK1, TLorentzVector *dauK2, TLorentzVector *daupi) const
unsigned int UInt_t
Definition: External.C:33
AliPIDCombined * GetPidCombined() const
Definition: AliAODPidHF.h:161
Double_t fMaxAngleForRot
minimum angle for track rotation
Int_t GetPoolIndex(Double_t zvert, Double_t mult)
Double_t fMinAngleForRot3
number of rotations (3rd prong)
AliESDtrackCuts * fTrackCutsAll
FilterMask.
TH3F * fPtVsYVsMultGenAccEvSelFeeddw
! hist. of Y vs. Pt vs. Mult generated (D in acc, sel ev.)
AliRDHFCuts * fAnalysisCuts
PID configuration.
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
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.)
Bool_t IsEventRejectedDueToPileup() const
Definition: AliRDHFCuts.h:327
Double_t nsigma
Double_t nEvents
plot quality messages
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:160
virtual void Terminate(Option_t *option)
TH3F * fPtVsYVsMultGenFeeddw
! hist. of Y vs. Pt vs. Mult generated (all D)
Double_t ComputeInvMassKK(AliAODTrack *tr1, AliAODTrack *tr2) const
Double_t fMinMultiplicity
multiplicity
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 IsEventRejectedDuePhysicsSelection() const
Definition: AliRDHFCuts.h:345
Double_t fMinMass
Cuts for candidates.
Int_t fPIDstrategy
flag to set analysis level (0 is the fastest)
Double_t fVtxZ
unique event Id for event mixing checks
AliESDtrackCuts * fTrackCutsKaon
pion track selection
Bool_t fGoUpToQuark
flag for access to MC
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)
TH1F * fHistCheckDecChanAcc
!hist. of decay channel of D meson in acc.
Int_t fDoEventMixing
threshold for pion identification via Bayesian PID
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)
Double_t fmaxPforIDKaon
flag for upper p limit for id band for pion
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
TObjArray * fPionTracks
array of kaon-compatible tracks (TLorentzVectors)
void DoMixingWithPools(Int_t poolIndex)
Int_t fNRotations
pt limits for acceptance step
const char Option_t
Definition: External.C:48
Double_t fBayesThresKaon
flag to change PID strategy
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:315
TH1F * fDeltaMass
! hist. mass difference after rotations
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:310
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:279
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
Double_t fMaxAngleForRot3
minimum angle for track rotation (3rd prong)
void SetPidResponse(AliPIDResponse *pidResp)
Definition: AliAODPidHF.h:111
TH3F * fMassVsPtVsYRot
! hist. of Y vs. Pt vs. Mass (rotations)
TH1F * fHistCheckOrigin
!hist. of origin (c/b) of D meson