AliPhysics  aaf9c62 (aaf9c62)
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 ",10,-0.5,9.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 zSPD-zTrack");
349  fHistNEvents->GetXaxis()->SetBinLabel(8,"n. rejected for vertex out of accept");
350  fHistNEvents->GetXaxis()->SetBinLabel(9,"n. rejected for pileup");
351  fHistNEvents->GetXaxis()->SetBinLabel(10,"no. of out centrality events");
352 
353  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
354  fHistNEvents->SetMinimum(0);
355  fOutput->Add(fHistNEvents);
356 
357  fHistEventMultZv = new TH2F("hEventMultZv","",30,-15.,15.,200,fMinMultiplicity,fMaxMultiplicity);
359 
360  fHistEventMultZvEvSel = new TH2F("hEventMultZvEvSel","",30,-15.,15.,200,fMinMultiplicity,fMaxMultiplicity);
362 
363  fHistTrackStatus = new TH1F("hTrackStatus", "",8,-0.5,7.5);
364  fHistTrackStatus->GetXaxis()->SetBinLabel(1,"Not OK");
365  fHistTrackStatus->GetXaxis()->SetBinLabel(2,"Track OK");
366  fHistTrackStatus->GetXaxis()->SetBinLabel(3,"Kaon, Not OK");
367  fHistTrackStatus->GetXaxis()->SetBinLabel(4,"Kaon OK");
368  fHistTrackStatus->GetXaxis()->SetBinLabel(5,"Pion, Not OK");
369  fHistTrackStatus->GetXaxis()->SetBinLabel(6,"Pion OK");
370  fHistTrackStatus->GetXaxis()->SetBinLabel(7,"Kaon||Pion, Not OK");
371  fHistTrackStatus->GetXaxis()->SetBinLabel(8,"Kaon||Pion OK");
372  fHistTrackStatus->GetXaxis()->SetNdivisions(1,kFALSE);
373  fHistTrackStatus->SetMinimum(0);
375 
376  fHistTrackEtaMultZv = new TH3F("hTrackEtaMultZv","",40,-1.,1.,30,-15.,15.,200,fMinMultiplicity,fMaxMultiplicity);
378 
379  Int_t nPtBins = (Int_t)(fMaxPt/fPtBinWidth+0.001);
381 
382  if(fReadMC){
383 
384  fHistCheckOrigin=new TH1F("hCheckOrigin","",7,-1.5,5.5);
385  fHistCheckOrigin->SetMinimum(0);
387 
388  fHistCheckDecChan=new TH1F("hCheckDecChan","",7,-2.5,4.5);
389  fHistCheckDecChan->SetMinimum(0);
391 
392  fHistCheckDecChanAcc=new TH1F("hCheckDecChanAcc","",7,-2.5,4.5);
393  fHistCheckDecChanAcc->SetMinimum(0);
395 
396  fPtVsYVsMultGenPrompt = new TH3F("hPtVsYVsMultGenPrompt","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
397  fPtVsYVsMultGenPrompt->Sumw2();
398  fPtVsYVsMultGenPrompt->SetMinimum(0);
400 
401  fPtVsYVsMultGenLargeAccPrompt = new TH3F("hPtVsYVsMultGenLargeAccPrompt","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
403  fPtVsYVsMultGenLargeAccPrompt->SetMinimum(0);
405 
406  fPtVsYVsMultGenLimAccPrompt = new TH3F("hPtVsYVsMultGenLimAccPrompt","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
408  fPtVsYVsMultGenLimAccPrompt->SetMinimum(0);
410 
411  fPtVsYVsMultGenAccPrompt = new TH3F("hPtVsYVsMultGenAccPrompt","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
412  fPtVsYVsMultGenAccPrompt->Sumw2();
413  fPtVsYVsMultGenAccPrompt->SetMinimum(0);
415 
416  fPtVsYVsMultGenAccEvSelPrompt = new TH3F("hPtVsYVsMultGenAccEvSelPrompt","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
418  fPtVsYVsMultGenAccEvSelPrompt->SetMinimum(0);
420 
421  fPtVsYVsMultRecoPrompt = new TH3F("hPtVsYVsMultRecoPrompt","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
422  fPtVsYVsMultRecoPrompt->Sumw2();
423  fPtVsYVsMultRecoPrompt->SetMinimum(0);
425 
426  fPtVsYVsMultGenFeeddw = new TH3F("hPtVsYVsMultGenFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
427  fPtVsYVsMultGenFeeddw->Sumw2();
428  fPtVsYVsMultGenFeeddw->SetMinimum(0);
430 
431  fPtVsYVsMultGenLargeAccFeeddw = new TH3F("hPtVsYVsMultGenLargeAccFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
433  fPtVsYVsMultGenLargeAccFeeddw->SetMinimum(0);
435 
436  fPtVsYVsMultGenLimAccFeeddw = new TH3F("hPtVsYVsMultGenLimAccFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
438  fPtVsYVsMultGenLimAccFeeddw->SetMinimum(0);
440 
441  fPtVsYVsMultGenAccFeeddw = new TH3F("hPtVsYVsMultGenAccFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
442  fPtVsYVsMultGenAccFeeddw->Sumw2();
443  fPtVsYVsMultGenAccFeeddw->SetMinimum(0);
445 
446  fPtVsYVsMultGenAccEvSelFeeddw = new TH3F("hPtVsYVsMultGenAccEvSelFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
448  fPtVsYVsMultGenAccEvSelFeeddw->SetMinimum(0);
450 
451  fPtVsYVsMultRecoFeeddw = new TH3F("hPtVsYVsMultRecoFeeddw","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
452  fPtVsYVsMultRecoFeeddw->Sumw2();
453  fPtVsYVsMultRecoFeeddw->SetMinimum(0);
455 
456  }
457 
458 
459  Int_t nMassBins=static_cast<Int_t>(fMaxMass*1000.-fMinMass*1000.);
460  Double_t maxm=fMinMass+nMassBins*0.001;
461  fMassVsPtVsY=new TH3F("hMassVsPtVsY","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
462  fMassVsPtVsY->Sumw2();
463  fMassVsPtVsY->SetMinimum(0);
464  fOutput->Add(fMassVsPtVsY);
465 
466  fMassVsPtVsYRot=new TH3F("hMassVsPtVsYRot","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
467  fMassVsPtVsYRot->Sumw2();
468  fMassVsPtVsYRot->SetMinimum(0);
469  fOutput->Add(fMassVsPtVsYRot);
470 
471  fMassVsPtVsYLSpp=new TH3F("hMassVsPtVsYLSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
472  fMassVsPtVsYLSpp->Sumw2();
473  fMassVsPtVsYLSpp->SetMinimum(0);
475  fMassVsPtVsYLSmm=new TH3F("hMassVsPtVsYLSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
476  fMassVsPtVsYLSmm->Sumw2();
477  fMassVsPtVsYLSmm->SetMinimum(0);
479 
480  fMassVsPtVsYSig=new TH3F("hMassVsPtVsYSig","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
481  fMassVsPtVsYSig->Sumw2();
482  fMassVsPtVsYSig->SetMinimum(0);
483  fOutput->Add(fMassVsPtVsYSig);
484 
485  fMassVsPtVsYRefl=new TH3F("hMassVsPtVsYRefl","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
486  fMassVsPtVsYRefl->Sumw2();
487  fMassVsPtVsYRefl->SetMinimum(0);
489 
490  fMassVsPtVsYBkg=new TH3F("hMassVsPtVsYBkg","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
491  fMassVsPtVsYBkg->Sumw2();
492  fMassVsPtVsYBkg->SetMinimum(0);
493  fOutput->Add(fMassVsPtVsYBkg);
494 
495  fNSelected=new TH1F("hNSelected","",100,-0.5,99.5);
496  fNSelected->Sumw2();
497  fNSelected->SetMinimum(0);
498  fOutput->Add(fNSelected);
499 
500  fNormRotated=new TH1F("hNormRotated","",11,-0.5,10.5);
501  fNormRotated->Sumw2();
502  fNormRotated->SetMinimum(0);
503  fOutput->Add(fNormRotated);
504 
505  fDeltaMass=new TH1F("hDeltaMass","",100,-0.4,0.4);
506  fDeltaMass->Sumw2();
507  fDeltaMass->SetMinimum(0);
508  fOutput->Add(fDeltaMass);
509 
510  Int_t binSparseDMassRot[5]={nMassBins,100,24,40,20};
511  Double_t edgeLowSparseDMassRot[5]={fMinMass,-0.4,0.,-4.,0};
512  Double_t edgeHighSparseDMassRot[5]={maxm,0.4,12.,4.,3.14};
513  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);
515 
516  fMassVsPtVsYME=new TH3F("hMassVsPtVsYME","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
517  fMassVsPtVsYME->Sumw2();
518  fMassVsPtVsYME->SetMinimum(0);
519  fOutput->Add(fMassVsPtVsYME);
520 
521  fMassVsPtVsYMELSpp=new TH3F("hMassVsPtVsYMELSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
522  fMassVsPtVsYMELSpp->Sumw2();
523  fMassVsPtVsYMELSpp->SetMinimum(0);
525 
526  fMassVsPtVsYMELSmm=new TH3F("hMassVsPtVsYMELSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
527  fMassVsPtVsYMELSmm->Sumw2();
528  fMassVsPtVsYMELSmm->SetMinimum(0);
530 
533  if(fDoEventMixing==2) fNOfPools=1;
535  fEventsPerPool=new TH2F("hEventsPerPool","hEventsPerPool",fNzVertPools,fzVertPoolLims,fNMultPools,fMultPoolLims);
536  fMixingsPerPool=new TH2F("hMixingsPerPool","hMixingsPerPool",fNzVertPools,fzVertPoolLims,fNMultPools,fMultPoolLims);
537  }else{
538  fEventsPerPool=new TH2F("hEventsPerPool","hEventsPerPool",1,-10.,10.,1,-0.5,2000.5);
539  fMixingsPerPool=new TH2F("hMixingsPerPool","hMixingsPerPool",1,-10.,10.,1,-0.5,2000.5);
540  }
541  fEventsPerPool->Sumw2();
542  fEventsPerPool->SetMinimum(0);
543  fOutput->Add(fEventsPerPool);
544  fMixingsPerPool->Sumw2();
545  fMixingsPerPool->SetMinimum(0);
546  fOutput->Add(fMixingsPerPool);
547 
548  //Counter for Normalization
549  fCounter = new AliNormalizationCounter("NormalizationCounter");
550  fCounter->Init();
551 
552  fKaonTracks = new TObjArray();
553  fPionTracks=new TObjArray();
554  fKaonTracks->SetOwner();
555  fPionTracks->SetOwner();
556 
557  fEventBuffer = new TTree*[fNOfPools];
558  for(Int_t i=0; i<fNOfPools; i++){
559  fEventBuffer[i]=new TTree(Form("EventBuffer_%d",i), "Temporary buffer for event mixing");
560  fEventBuffer[i]->Branch("zVertex", &fVtxZ);
561  fEventBuffer[i]->Branch("multiplicity", &fMultiplicity);
562  fEventBuffer[i]->Branch("eventInfo", "TObjString",&fEventInfo);
563  fEventBuffer[i]->Branch("karray", "TObjArray", &fKaonTracks);
564  fEventBuffer[i]->Branch("parray", "TObjArray", &fPionTracks);
565  }
566 
567  PostData(1,fOutput);
568  PostData(2,fCounter);
569 }
570 
571 //________________________________________________________________________
574 
575  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
576  if(!aod && AODEvent() && IsStandardAOD()) {
577  // In case there is an AOD handler writing a standard AOD, use the AOD
578  // event in memory rather than the input (ESD) event.
579  aod = dynamic_cast<AliAODEvent*> (AODEvent());
580  }
581  if(!aod){
582  printf("AliAnalysisTaskCombinHF::UserExec: AOD not found!\n");
583  return;
584  }
585 
586  // fix for temporary bug in ESDfilter
587  // the AODs with null vertex pointer didn't pass the PhysSel
588  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
589 
590  // Reject events with trigger mask 0 of the LHC13d3 production
591  // For these events the ITS layers are skipped in the trakcing
592  // and the vertex reconstruction efficiency from tracks is biased
593  if(fReadMC){
594  Int_t runnumber = aod->GetRunNumber();
595  if(aod->GetTriggerMask()==0 &&
596  (runnumber>=195344 && runnumber<=195677)){
597  return;
598  }
599  }
600 
601  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
602  AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
603  AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
604  fPidHF->SetPidResponse(pidResp);
605 
606 
607  fHistNEvents->Fill(0); // count event
608  // Post the data already here
609  PostData(1,fOutput);
610 
612 
613  Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
617  }else{
619  fHistNEvents->Fill(9);
620  }else{
625  }else{
628  }
629  }
630  }
631 
633  // events not passing the centrality selection can be removed immediately. For the others we must count the generated D mesons
634 
635  Int_t ntracks=aod->GetNumberOfTracks();
636  fVtxZ = aod->GetPrimaryVertex()->GetZ();
638 
639  TClonesArray *arrayMC=0;
640  AliAODMCHeader *mcHeader=0;
641  if(fReadMC){
642  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
643  if(!arrayMC) {
644  printf("AliAnalysisTaskCombinHF::UserExec: MC particles branch not found!\n");
645  return;
646  }
647 
648  // load MC header
649  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
650  if(!mcHeader) {
651  printf("AliAnalysisTaskCombinHF::UserExec: MC header branch not found!\n");
652  return;
653  }
654  Double_t zMCVertex = mcHeader->GetVtxZ();
655  if (TMath::Abs(zMCVertex) < fAnalysisCuts->GetMaxVtxZ()){ // only cut on zVertex applied to count the signal
656  FillGenHistos(arrayMC,isEvSel);
657  }
658  fHistEventMultZv->Fill(zMCVertex,fMultiplicity);
659  if(isEvSel) fHistEventMultZvEvSel->Fill(zMCVertex,fMultiplicity);
660  }else{
662  if(isEvSel) fHistEventMultZvEvSel->Fill(fVtxZ,fMultiplicity);
663  }
664 
665 
666  if(!isEvSel)return;
667 
668  fHistNEvents->Fill(1);
669 
670 
671  // select and flag tracks
672  UChar_t* status = new UChar_t[ntracks];
673  for(Int_t iTr=0; iTr<ntracks; iTr++){
674  status[iTr]=0;
675  AliAODTrack* track=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr));
676  if(!track){
677  AliWarning("Error in casting track to AOD track. Not a standard AOD?");
678  continue;
679  }
680  if(IsTrackSelected(track)) status[iTr]+=1;
681 
682  // PID
683  if (fPIDstrategy == knSigma) {
684  // nsigma PID
685  if(IsKaon(track)) status[iTr]+=2;
686  if(IsPion(track)) status[iTr]+=4;
687  }
689  // Bayesian PID
690  Double_t *weights = new Double_t[AliPID::kSPECIES];
691  fPidHF->GetPidCombined()->ComputeProbabilities(track, fPidHF->GetPidResponse(), weights);
693  if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kKaon]) status[iTr] += 2;
694  if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kPion]) status[iTr] += 4;
695  }
696  if (fPIDstrategy == kBayesianThres) {
697  if (weights[AliPID::kKaon] > fBayesThresKaon) status[iTr] += 2;
698  if (weights[AliPID::kPion] > fBayesThresPion) status[iTr] += 4;
699  }
700  delete[] weights;
701  }
702 
703  fHistTrackStatus->Fill(status[iTr]);
704  fHistTrackEtaMultZv->Fill(track->Eta(),fVtxZ,fMultiplicity);
705  }
706 
707  // build the combinatorics
708  Int_t nSelected=0;
709  Int_t nFiltered=0;
710  Double_t dummypos[3]={0.,0.,0.};
711  AliAODVertex* v2=new AliAODVertex(dummypos,999.,-1,2);
712  AliAODVertex* v3=new AliAODVertex(dummypos,999.,-1,3);
713  // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
714  Double_t d02[2]={0.,0.};
715  Double_t d03[3]={0.,0.,0.};
716  AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
717  AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
718  UInt_t pdg0[2]={321,211};
719  UInt_t pdgp[3]={321,211,211};
720  UInt_t pdgs[3]={321,211,321};
721  Double_t tmpp[3];
722  Double_t px[3],py[3],pz[3];
723  Int_t dgLabels[3];
724  fKaonTracks->Delete();
725  fPionTracks->Delete();
726 
727  for(Int_t iTr1=0; iTr1<ntracks; iTr1++){
728  AliAODTrack* trK=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr1));
729  if(!trK){
730  AliWarning("Error in casting track to AOD track. Not a standard AOD?");
731  continue;
732  }
733  if((status[iTr1] & 1)==0) continue;
734  if(fDoEventMixing>0){
735  if(status[iTr1] & 2) fKaonTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
736  if(status[iTr1] & 4) fPionTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
737  }
738  if((status[iTr1] & 2)==0) continue;
739  Int_t chargeK=trK->Charge();
740  trK->GetPxPyPz(tmpp);
741  px[0] = tmpp[0];
742  py[0] = tmpp[1];
743  pz[0] = tmpp[2];
744  dgLabels[0]=trK->GetLabel();
745  for(Int_t iTr2=0; iTr2<ntracks; iTr2++){
746  if((status[iTr2] & 1)==0) continue;
747  if((status[iTr2] & 4)==0) continue;
748  if(iTr1==iTr2) continue;
749  AliAODTrack* trPi1=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr2));
750  if(!trPi1){
751  AliWarning("Error in casting track to AOD track. Not a standard AOD?");
752  continue;
753  }
754  Int_t chargePi1=trPi1->Charge();
755  trPi1->GetPxPyPz(tmpp);
756  px[1] = tmpp[0];
757  py[1] = tmpp[1];
758  pz[1] = tmpp[2];
759  dgLabels[1]=trPi1->GetLabel();
760  if(fMeson==kDzero){
761  if(chargePi1==chargeK){
762  // LS candidate
763  FillLSHistos(421,2,tmpRD2,px,py,pz,pdg0,chargePi1);
764  }else{
765  // OS candidate
766  nFiltered++;
767  v2->AddDaughter(trK);
768  v2->AddDaughter(trPi1);
769  tmpRD2->SetSecondaryVtx(v2);
770  Bool_t ok=FillHistos(421,2,tmpRD2,px,py,pz,pdg0,arrayMC,dgLabels);
771  v2->RemoveDaughters();
772  if(ok) nSelected++;
773  }
774  }else{
775  for(Int_t iTr3=iTr2+1; iTr3<ntracks; iTr3++){
776  if((status[iTr3] & 1)==0) continue;
777  if(fMeson==kDplus && (status[iTr3] & 4)==0) continue;
778  if(fMeson==kDs && (status[iTr3] & 2)==0) continue;
779  if(iTr1==iTr3) continue;
780  AliAODTrack* trPi2=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr3));
781  if(!trPi2){
782  AliWarning("Error in casting track to AOD track. Not a standard AOD?");
783  continue;
784  }
785  Int_t chargePi2=trPi2->Charge();
786  trPi2->GetPxPyPz(tmpp);
787  px[2] = tmpp[0];
788  py[2] = tmpp[1];
789  pz[2] = tmpp[2];
790  dgLabels[2]=trPi2->GetLabel();
791  if(fMeson==kDs){
792  Double_t massKK=ComputeInvMassKK(trK,trPi2);
793  Double_t deltaMass=massKK-TDatabasePDG::Instance()->GetParticle(333)->Mass();
794  if(TMath::Abs(deltaMass)>fPhiMassCut) continue;
795  tmpRD3->SetPxPyPzProngs(3,px,py,pz);
796  Double_t cos1=((AliAODRecoDecayHF3Prong*)tmpRD3)->CosPiKPhiRFrameKpiK();
797  Double_t kincutPiKPhi=TMath::Abs(cos1*cos1*cos1);
798  if(kincutPiKPhi<fCutCos3PiKPhiRFrame) continue;
799  Double_t cosPiDsLabFrame=((AliAODRecoDecayHF3Prong*)tmpRD3)->CosPiDsLabFrameKpiK();
800  if(cosPiDsLabFrame>fCutCosPiDsLabFrame) continue;
801  }
802  Bool_t isThreeLS=kFALSE;
803  if(chargePi1==chargeK && chargePi2==chargeK){
804  isThreeLS=kTRUE;
805  if(fMeson==kDplus)FillLSHistos(411,3,tmpRD3,px,py,pz,pdgp,chargePi1);
806  else if(fMeson==kDs)FillLSHistos(431,3,tmpRD3,px,py,pz,pdgs,chargePi1);
807  }
808  Bool_t acceptOS=kFALSE;
809  if(fMeson==kDplus){
810  if(chargePi1!=chargeK && chargePi2!=chargeK)acceptOS=kTRUE;
811  }else if(fMeson==kDs){
812  if(chargePi2!=chargeK && !isThreeLS) acceptOS=kTRUE;
813  }
814  if(acceptOS){
815  nFiltered++;
816  v3->AddDaughter(trK);
817  v3->AddDaughter(trPi1);
818  v3->AddDaughter(trPi2);
819  tmpRD3->SetSecondaryVtx(v3);
820  Bool_t ok=kFALSE;
821  if(fMeson==kDplus) ok=FillHistos(411,3,tmpRD3,px,py,pz,pdgp,arrayMC,dgLabels);
822  else if(fMeson==kDs) ok=FillHistos(431,3,tmpRD3,px,py,pz,pdgs,arrayMC,dgLabels);
823  v3->RemoveDaughters();
824  if(ok) nSelected++;
825  }
826  }
827  }
828  }
829  }
830 
831  delete [] status;
832  delete v2;
833  delete v3;
834  delete tmpRD2;
835  delete tmpRD3;
836 
837  fNSelected->Fill(nSelected);
838 
839  fCounter->StoreCandidates(aod,nFiltered,kTRUE);
840  fCounter->StoreCandidates(aod,nSelected,kFALSE);
841  fEventInfo->SetString(Form("Ev%d_esd%d_Pi%d_K%d",mgr->GetNcalls(),((AliAODHeader*)aod->GetHeader())->GetEventNumberESDFile(),fPionTracks->GetEntries(),fKaonTracks->GetEntries()));
842  if(fDoEventMixing==1){
844  if(ind>=0 && ind<fNOfPools){
846  fEventBuffer[ind]->Fill();
847  if(fEventBuffer[ind]->GetEntries() >= fNumberOfEventsForMixing){
849  DoMixingWithPools(ind);
850  ResetPool(ind);
851  }
852  }
853  }else if(fDoEventMixing==2){ // mix with cuts, no pools
854  fEventBuffer[0]->Fill();
855  }
856  PostData(1,fOutput);
857  PostData(2,fCounter);
858 
859  return;
860 }
861 
862 //________________________________________________________________________
863 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){
865 
866  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
867  Double_t pt = tmpRD->Pt();
868  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
869  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
870  Double_t rapid = tmpRD->Y(pdgD);
871  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
872  if(charge>0) fMassVsPtVsYLSpp->Fill(TMath::Sqrt(minv2),pt,rapid);
873  else fMassVsPtVsYLSmm->Fill(TMath::Sqrt(minv2),pt,rapid);
874  }
875  }
876  return;
877 }
878 
879 //________________________________________________________________________
880 void AliAnalysisTaskCombinHF::FillGenHistos(TClonesArray* arrayMC, Bool_t isEvSel){
882  Int_t totPart=arrayMC->GetEntriesFast();
883  Int_t thePDG=411;
884  Int_t nProng=3;
885  if(fMeson==kDzero){
886  thePDG=421;
887  nProng=2;
888  }else if(fMeson==kDs){
889  thePDG=431;
890  nProng=3;
891  }
892  for(Int_t ip=0; ip<totPart; ip++){
893  AliAODMCParticle *part = (AliAODMCParticle*)arrayMC->At(ip);
894  if(TMath::Abs(part->GetPdgCode())==thePDG){
896  fHistCheckOrigin->Fill(orig);
897  Int_t deca=0;
898  Bool_t isGoodDecay=kFALSE;
899  Int_t labDau[4]={-1,-1,-1,-1};
900  if(fMeson==kDzero){
901  deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,part,labDau);
902  if(part->GetNDaughters()!=2) continue;
903  if(deca==1) isGoodDecay=kTRUE;
904  }else if(fMeson==kDplus){
905  deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,part,labDau);
906  if(deca>0) isGoodDecay=kTRUE;
907  }else if(fMeson==kDs){
908  deca=AliVertexingHFUtils::CheckDsDecay(arrayMC,part,labDau);
909  if(deca==1) isGoodDecay=kTRUE;
910  }
911  fHistCheckDecChan->Fill(deca);
912  if(labDau[0]==-1){
913  // printf(Form("Meson %d Label of daughters not filled correctly -- %d\n",fMeson,isGoodDecay));
914  continue; //protection against unfilled array of labels
915  }
916  Bool_t isInAcc=CheckAcceptance(arrayMC,nProng,labDau);
917  if(isInAcc) fHistCheckDecChanAcc->Fill(deca);
918  if(isGoodDecay){
919  Double_t ptgen=part->Pt();
920  Double_t ygen=part->Y();
921  if(fAnalysisCuts->IsInFiducialAcceptance(ptgen,ygen)){
922  if(orig==4){
923  fPtVsYVsMultGenPrompt->Fill(ptgen,ygen,fMultiplicity);
924  if(TMath::Abs(ygen)<0.5) fPtVsYVsMultGenLimAccPrompt->Fill(ptgen,ygen,fMultiplicity);
925  if(isInAcc) fPtVsYVsMultGenAccPrompt->Fill(ptgen,ygen,fMultiplicity);
926  if(isEvSel && isInAcc) fPtVsYVsMultGenAccEvSelPrompt->Fill(ptgen,ygen,fMultiplicity);
927  }else if(orig==5){
928  fPtVsYVsMultGenFeeddw->Fill(ptgen,ygen,fMultiplicity);
929  if(TMath::Abs(ygen)<0.5) fPtVsYVsMultGenLimAccFeeddw->Fill(ptgen,ygen,fMultiplicity);
930  if(isInAcc) fPtVsYVsMultGenAccFeeddw->Fill(ptgen,ygen,fMultiplicity);
931  if(isEvSel && isInAcc) fPtVsYVsMultGenAccEvSelFeeddw->Fill(ptgen,ygen,fMultiplicity);
932  }
933  }
934  if(TMath::Abs(ygen)<0.9){
935  if(orig==4) fPtVsYVsMultGenLargeAccPrompt->Fill(ptgen,ygen,fMultiplicity);
936  else if(orig==5) fPtVsYVsMultGenLargeAccFeeddw->Fill(ptgen,ygen,fMultiplicity);
937  }
938  }
939  }
940  }
941 }
942 
943 //________________________________________________________________________
944 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){
946 
947  Bool_t accept=kFALSE;
948 
949  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
950  Double_t pt = tmpRD->Pt();
951  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
952  Double_t mass=TMath::Sqrt(minv2);
953 
954  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
955  Double_t rapid = tmpRD->Y(pdgD);
956  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
957  fMassVsPtVsY->Fill(mass,pt,rapid);
958  accept=kTRUE;
959  if(fReadMC){
960  Int_t signPdg[3]={0,0,0};
961  for(Int_t iii=0; iii<nProngs; iii++) signPdg[iii]=pdgdau[iii];
962  Int_t labD = tmpRD->MatchToMC(pdgD,arrayMC,nProngs,signPdg);
963  if(labD>=0){
964  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(TMath::Abs(dgLabels[0])));
965  if(part){
967  Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
968  if(pdgCode==321){
969  fMassVsPtVsYSig->Fill(mass,pt,rapid);
970  AliAODMCParticle* dmes = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labD));
971  if(dmes){
972  if(orig==4) fPtVsYVsMultRecoPrompt->Fill(dmes->Pt(),dmes->Y(),fMultiplicity);
973  else if(orig==5) fPtVsYVsMultRecoFeeddw->Fill(dmes->Pt(),dmes->Y(),fMultiplicity);
974  }
975  }else{
976  fMassVsPtVsYRefl->Fill(mass,pt,rapid);
977  }
978  }
979  }else{
980  fMassVsPtVsYBkg->Fill(mass,pt,rapid);
981  }
982  }
983  }
984  }
985 
986  Int_t nRotated=0;
987  Double_t massRot=0;// calculated later only if candidate is acceptable
988  Double_t angleProngXY;
989  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])));
990  else if(TMath::Abs(pdgD)==431) {
991  Double_t px_phi = px[0]+px[2];
992  Double_t py_phi = py[0]+py[2];
993  Double_t pz_phi = pz[0]+pz[2];
994  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])));
995  }
996  else {//angle between pion and phi meson
997  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])));
998  }
999  Double_t ptOrig=pt;
1000 
1001 
1002  Double_t rotStep=0.;
1003  if(fNRotations>1) rotStep=(fMaxAngleForRot-fMinAngleForRot)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
1004  if(TMath::Abs(pdgD)==421 || TMath::Abs(pdgD)==431) fNRotations3=1;
1005  Double_t rotStep3=0.;
1006  if(fNRotations3>1) rotStep3=(fMaxAngleForRot3-fMinAngleForRot3)/(fNRotations3-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
1007 
1008  for(Int_t irot=0; irot<fNRotations; irot++){
1009  Double_t phirot=fMinAngleForRot+rotStep*irot;
1010  Double_t tmpx=px[0];
1011  Double_t tmpy=py[0];
1012  Double_t tmpx2=px[2];
1013  Double_t tmpy2=py[2];
1014  if(pdgD==431) {
1015  //rotate pion w.r.t. phi meson
1016  tmpx=px[1];
1017  tmpy=py[1];
1018  px[1]=tmpx*TMath::Cos(phirot)-tmpy*TMath::Sin(phirot);
1019  py[1]=tmpx*TMath::Sin(phirot)+tmpy*TMath::Cos(phirot);
1020  }
1021  else {
1022  px[0]=tmpx*TMath::Cos(phirot)-tmpy*TMath::Sin(phirot);
1023  py[0]=tmpx*TMath::Sin(phirot)+tmpy*TMath::Cos(phirot);
1024  }
1025  for(Int_t irot3=0; irot3<fNRotations3; irot3++){
1026  if(pdgD==411){
1027  Double_t phirot2=fMaxAngleForRot3-rotStep3*irot;
1028  px[2]=tmpx*TMath::Cos(phirot2)-tmpy*TMath::Sin(phirot2);
1029  py[2]=tmpx*TMath::Sin(phirot2)+tmpy*TMath::Cos(phirot2);
1030  }
1031  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
1032  pt = tmpRD->Pt();
1033  minv2 = tmpRD->InvMass2(nProngs,pdgdau);
1034  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
1035  Double_t rapid = tmpRD->Y(pdgD);
1036  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
1037  massRot=TMath::Sqrt(minv2);
1038  fMassVsPtVsYRot->Fill(massRot,pt,rapid);
1039  nRotated++;
1040  fDeltaMass->Fill(massRot-mass);
1041  if(fFullAnalysis){
1042  Double_t pointRot[5]={mass,massRot-mass,ptOrig,pt-ptOrig,angleProngXY};
1043  fDeltaMassFullAnalysis->Fill(pointRot);
1044  }
1045  }
1046  }
1047  }
1048  if(pdgD==431) {
1049  px[1]=tmpx;
1050  py[1]=tmpy;
1051  }
1052  else {
1053  px[0]=tmpx;
1054  py[0]=tmpy;
1055  if(pdgD==411){
1056  px[2]=tmpx2;
1057  py[2]=tmpy2;
1058  }
1059  }
1060  }
1061  fNormRotated->Fill(nRotated);
1062 
1063  return accept;
1064 
1065 }
1066 //________________________________________________________________________
1067 void AliAnalysisTaskCombinHF::FillMEHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau){
1069 
1070  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
1071  Double_t pt = tmpRD->Pt();
1072  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
1073  Double_t mass=TMath::Sqrt(minv2);
1074 
1075  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
1076  Double_t rapid = tmpRD->Y(pdgD);
1077  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
1078  fMassVsPtVsYME->Fill(mass,pt,rapid);
1079  }
1080  }
1081  return;
1082 }
1083 //________________________________________________________________________
1084 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){
1086 
1087  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
1088  Double_t pt = tmpRD->Pt();
1089  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
1090  Double_t mass=TMath::Sqrt(minv2);
1091 
1092  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
1093  Double_t rapid = tmpRD->Y(pdgD);
1094  if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
1095  if(charge>0) fMassVsPtVsYMELSpp->Fill(mass,pt,rapid);
1096  if(charge<0) fMassVsPtVsYMELSmm->Fill(mass,pt,rapid);
1097  }
1098  }
1099  return;
1100 }
1101 //________________________________________________________________________
1104 
1105  if(track->Charge()==0) return kFALSE;
1106  if(track->GetID()<0&&!fKeepNegID)return kFALSE;
1107  if(fFilterMask>=0){
1108  if(!(track->TestFilterMask(fFilterMask))) return kFALSE;
1109  }
1110  if(!SelectAODTrack(track,fTrackCutsAll)) return kFALSE;
1111  return kTRUE;
1112 }
1113 
1114 //________________________________________________________________________
1117 
1118  if(!fPidHF) return kTRUE;
1119  Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
1120  Double_t mom=track->P();
1121  if(SelectAODTrack(track,fTrackCutsKaon)) {
1122  if(isKaon>=1) return kTRUE;
1123  if(isKaon<=-1) return kFALSE;
1124  switch(fPIDselCaseZero){// isKaon=0
1125  case 0:
1126  {
1127  return kTRUE;// always accept
1128  }
1129  break;
1130  case 1:
1131  {
1132  if(isKaon>=0 && track->P()>fmaxPforIDKaon)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDKaon
1133  }
1134  break;
1135  case 2:
1136  {
1137  if(track->P()>fmaxPforIDKaon)return kTRUE;
1138  AliPIDResponse *pidResp=fPidHF->GetPidResponse();// more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
1139  Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kKaon);
1140  if(nsigma>-2.&& nsigma<3. && mom<0.6)isKaon=1;
1141  else if(nsigma>-1.&& nsigma<3.&& mom<0.8)isKaon=1;
1142  if(isKaon==1)return kTRUE;
1143  }
1144  break;
1145  default:
1146  {
1147  AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
1148  return kFALSE;// actually case 0 could be set as the default and return kTRUE
1149  }
1150  }
1151  }
1152 
1153  return kFALSE;
1154 }
1155 //_______________________________________________________________________
1158 
1159  if(!fPidHF) return kTRUE;
1160  Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
1161  Double_t mom=track->P();
1162  if(SelectAODTrack(track,fTrackCutsPion)) {
1163  if(isPion>=1) return kTRUE;
1164  if(isPion<=-1) return kFALSE;
1165  switch(fPIDselCaseZero){// isPion=0
1166  case 0:
1167  {
1168  return kTRUE;// always accept
1169  }
1170  break;
1171  case 1:
1172  {
1173  if(track->P()>fmaxPforIDPion)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDPion
1174  }
1175  break;
1176  case 2:
1177  {
1178  // more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
1179  if(track->P()>fmaxPforIDPion)return kTRUE;
1180  AliPIDResponse *pidResp=fPidHF->GetPidResponse();
1181  Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kPion);
1182  if(nsigma<2.&& nsigma>-3. && mom<0.6)isPion=1;
1183  else if(nsigma<1. && nsigma> -3. && mom<0.8)isPion=1;
1184  if(isPion==1)return kTRUE;
1185  }
1186  break;
1187  default:
1188  {
1189  AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
1190  return kFALSE;// actually case 0 could be set as the default and return kTRUE
1191  }
1192  }
1193  }
1194 
1195  return kFALSE;
1196 }
1197 
1198 //________________________________________________________________________
1199 Bool_t AliAnalysisTaskCombinHF::SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts){
1201 
1202  if(!cuts) return kTRUE;
1203 
1204  AliESDtrack esdTrack(track);
1205  // set the TPC cluster info
1206  esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
1207  esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
1208  esdTrack.SetTPCPointsF(track->GetTPCNclsF());
1209  if(!cuts->IsSelected(&esdTrack)) return kFALSE;
1210 
1211  return kTRUE;
1212 }
1213 
1214 //_________________________________________________________________
1215 Bool_t AliAnalysisTaskCombinHF::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1217  for (Int_t iProng = 0; iProng<nProng; iProng++){
1218  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1219  if(!mcPartDaughter) return kFALSE;
1220  Double_t eta = mcPartDaughter->Eta();
1221  Double_t pt = mcPartDaughter->Pt();
1222  if (TMath::Abs(eta) > fEtaAccCut || pt < fPtAccCut) return kFALSE;
1223  }
1224  return kTRUE;
1225 }
1226 //_________________________________________________________________
1229  if(!fzVertPoolLims || !fMultPoolLims) return 0;
1230  Int_t theBinZ=TMath::BinarySearch(fNzVertPoolsLimSize,fzVertPoolLims,zvert);
1231  if(theBinZ<0 || theBinZ>=fNzVertPoolsLimSize) return -1;
1232  Int_t theBinM=TMath::BinarySearch(fNMultPoolsLimSize,fMultPoolLims,mult);
1233  if(theBinM<0 || theBinM>=fNMultPoolsLimSize) return -1;
1234  return fNMultPools*theBinZ+theBinM;
1235 }
1236 //_________________________________________________________________
1239  if(poolIndex<0 || poolIndex>=fNOfPools) return;
1240  delete fEventBuffer[poolIndex];
1241  fEventBuffer[poolIndex]=new TTree(Form("EventBuffer_%d",poolIndex), "Temporary buffer for event mixing");
1242  fEventBuffer[poolIndex]->Branch("zVertex", &fVtxZ);
1243  fEventBuffer[poolIndex]->Branch("multiplicity", &fMultiplicity);
1244  fEventBuffer[poolIndex]->Branch("eventInfo", "TObjString",&fEventInfo);
1245  fEventBuffer[poolIndex]->Branch("karray", "TObjArray", &fKaonTracks);
1246  fEventBuffer[poolIndex]->Branch("parray", "TObjArray", &fPionTracks);
1247  return;
1248 }
1249 //_________________________________________________________________
1252  if(TMath::Abs(zv2-zv1)>fMaxzVertDistForMix) return kFALSE;
1253  if(TMath::Abs(mult2-mult1)>fMaxMultDiffForMix) return kFALSE;
1254  return kTRUE;
1255 }
1256 //_________________________________________________________________
1259 
1260  if(fDoEventMixing==0) return;
1261  Int_t nEvents=fEventBuffer[0]->GetEntries();
1262  if(fDebug > 1) printf("AnalysisTaskCombinHF::DoMixingWithCuts Start Event Mixing of %d events\n",nEvents);
1263 
1264  TObjArray* karray=0x0;
1265  TObjArray* parray=0x0;
1266  Double_t zVertex,mult;
1267  TObjString* eventInfo=0x0;
1268  fEventBuffer[0]->SetBranchAddress("karray", &karray);
1269  fEventBuffer[0]->SetBranchAddress("parray", &parray);
1270  fEventBuffer[0]->SetBranchAddress("eventInfo",&eventInfo);
1271  fEventBuffer[0]->SetBranchAddress("zVertex", &zVertex);
1272  fEventBuffer[0]->SetBranchAddress("multiplicity", &mult);
1273  Double_t d02[2]={0.,0.};
1274  Double_t d03[3]={0.,0.,0.};
1275  AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
1276  AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
1277  UInt_t pdg0[2]={321,211};
1278  // UInt_t pdgp[3]={321,211,211};
1279  Double_t px[3],py[3],pz[3];
1280  Int_t evId1,esdId1,nk1,np1;
1281  Int_t evId2,esdId2,nk2,np2;
1282 
1283  for(Int_t iEv1=0; iEv1<nEvents; iEv1++){
1284  fEventBuffer[0]->GetEvent(iEv1);
1285  TObjArray* karray1=(TObjArray*)karray->Clone();
1286  Double_t zVertex1=zVertex;
1287  Double_t mult1=mult;
1288  Int_t nKaons=karray1->GetEntries();
1289  Int_t nPionsForCheck=parray->GetEntries();
1290  sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId1,&esdId1,&np1,&nk1);
1291  if(nk1!=nKaons || np1!=nPionsForCheck){
1292  printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: read event does not match to the stored one\n");
1293  delete karray1;
1294  continue;
1295  }
1296  for(Int_t iEv2=0; iEv2<fNumberOfEventsForMixing; iEv2++){
1297  Int_t iToMix=iEv1+iEv2+1;
1298  if(iEv1>=(nEvents-fNumberOfEventsForMixing)) iToMix=iEv1-iEv2-1;
1299  if(iToMix<0) continue;
1300  if(iToMix==iEv1) continue;
1301  if(iToMix<iEv1) continue;
1302  fEventBuffer[0]->GetEvent(iToMix);
1303  Double_t zVertex2=zVertex;
1304  Double_t mult2=mult;
1305  if(TMath::Abs(zVertex2-zVertex1)<0.0001 && TMath::Abs(mult2-mult1)<0.001){
1306  printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: same event in mixing??? %d %d %f %f %f %f\n",iEv1,iEv2,zVertex1,zVertex2,mult1,mult2);
1307  continue;
1308  }
1309  TObjArray* parray2=(TObjArray*)parray->Clone();
1310  Int_t nPions=parray2->GetEntries();
1311  Int_t nKaonsForCheck=karray->GetEntries();
1312  sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId2,&esdId2,&np2,&nk2);
1313  if(nk2!=nKaonsForCheck || np2!=nPions){
1314  printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: read event does not match to the stored one\n");
1315  delete parray2;
1316  continue;
1317  }
1318  if(evId2==evId1 && esdId2==esdId1){
1319  printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: same event in mixing??? %d %d nK=%d %d nPi=%d %d\n",evId1,evId2,nKaons,nKaonsForCheck,nPionsForCheck,nPions);
1320  delete parray2;
1321  continue;
1322  }
1323  if(CanBeMixed(zVertex1,zVertex2,mult1,mult2)){
1324  for(Int_t iTr1=0; iTr1<nKaons; iTr1++){
1325  TLorentzVector* trK=(TLorentzVector*)karray1->At(iTr1);
1326  Double_t chargeK=trK->T();
1327  px[0] = trK->Px();
1328  py[0] = trK->Py();
1329  pz[0] = trK->Pz();
1330  for(Int_t iTr2=0; iTr2<nPions; iTr2++){
1331  TLorentzVector* trPi1=(TLorentzVector*)parray2->At(iTr2);
1332  Double_t chargePi1=trPi1->T();
1333  px[1] = trPi1->Px();
1334  py[1] = trPi1->Py();
1335  pz[1] = trPi1->Pz();
1336  if(fMeson==kDzero && chargePi1*chargeK<0){
1337  FillMEHistos(421,2,tmpRD2,px,py,pz,pdg0);
1338  }
1339  if(fMeson==kDzero && chargePi1*chargeK>0){
1340  FillMEHistosLS(421,2,tmpRD2,px,py,pz,pdg0,(Int_t)chargePi1);
1341  }
1342  }
1343  }
1344  }
1345  delete parray2;
1346  }
1347  delete karray1;
1348  }
1349  delete tmpRD2;
1350  delete tmpRD3;
1351 }
1352 //_________________________________________________________________
1355 
1356  if(fDoEventMixing==0) return;
1357  if(poolIndex<0 || poolIndex>fNzVertPools*fNMultPools) return;
1358 
1359  Int_t nEvents=fEventBuffer[poolIndex]->GetEntries();
1360  if(fDebug > 1) printf("AliAnalysisTaskCombinHF::DoMixingWithPools Start Event Mixing of %d events\n",nEvents);
1361  TObjArray* karray=0x0;
1362  TObjArray* parray=0x0;
1363  Double_t zVertex,mult;
1364  TObjString* eventInfo=0x0;
1365  fEventBuffer[poolIndex]->SetBranchAddress("karray", &karray);
1366  fEventBuffer[poolIndex]->SetBranchAddress("parray", &parray);
1367  fEventBuffer[poolIndex]->SetBranchAddress("eventInfo",&eventInfo);
1368  fEventBuffer[poolIndex]->SetBranchAddress("zVertex", &zVertex);
1369  fEventBuffer[poolIndex]->SetBranchAddress("multiplicity", &mult);
1370 
1371  // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
1372  Double_t d02[2]={0.,0.};
1373  Double_t d03[3]={0.,0.,0.};
1374  AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
1375  AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
1376  UInt_t pdg0[2]={321,211};
1377  UInt_t pdgp[3]={321,211,211};
1378  UInt_t pdgs[3]={321,211,321};
1379  Double_t px[3],py[3],pz[3];
1380  Int_t evId1,esdId1,nk1,np1;
1381  Int_t evId2,esdId2,nk2,np2;
1382 
1383  for(Int_t iEv1=0; iEv1<nEvents; iEv1++){
1384  fEventBuffer[poolIndex]->GetEvent(iEv1);
1385  TObjArray* karray1=(TObjArray*)karray->Clone();
1386  Double_t zVertex1=zVertex;
1387  Double_t mult1=mult;
1388  Int_t nKaons=karray1->GetEntries();
1389  Int_t nPionsForCheck=parray->GetEntries();
1390  sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId1,&esdId1,&np1,&nk1);
1391  if(nk1!=nKaons || np1!=nPionsForCheck){
1392  printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: read event does not match to the stored one\n");
1393  delete karray1;
1394  continue;
1395  }
1396  for(Int_t iEv2=0; iEv2<nEvents; iEv2++){
1397  if(iEv2==iEv1) continue;
1398  fEventBuffer[poolIndex]->GetEvent(iEv2);
1399  Double_t zVertex2=zVertex;
1400  Double_t mult2=mult;
1401  if(TMath::Abs(zVertex2-zVertex1)<0.0001 && TMath::Abs(mult2-mult1)<0.001){
1402  printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: same event in mixing??? %d %d %f %f %f %f\n",iEv1,iEv2,zVertex1,zVertex2,mult1,mult2);
1403  continue;
1404  }
1405  TObjArray* parray2=(TObjArray*)parray->Clone();
1406  Int_t nPions=parray2->GetEntries();
1407  Int_t nKaonsForCheck=karray->GetEntries();
1408  sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId2,&esdId2,&np2,&nk2);
1409  if(nk2!=nKaonsForCheck || np2!=nPions){
1410  printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: read event does not match to the stored one\n");
1411  delete parray2;
1412  continue;
1413  }
1414  if(evId2==evId1 && esdId2==esdId1){
1415  printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: same event in mixing??? %d %d nK=%d %d nPi=%d %d\n",evId1,evId2,nKaons,nKaonsForCheck,nPionsForCheck,nPions);
1416  delete parray2;
1417  continue;
1418  }
1419  TObjArray* parray3=0x0;
1420  Int_t nPions3=0;
1421  if(fMeson!=kDzero && fMeson!=kDs){
1422  Int_t iEv3=iEv2+1;
1423  if(iEv3==iEv1) iEv3=iEv2+2;
1424  if(iEv3>=nEvents) iEv3=iEv2-3;
1425  if(nEvents==2) iEv3=iEv1;
1426  if(iEv3<0) iEv3=iEv2-1;
1427  fEventBuffer[poolIndex]->GetEvent(iEv3);
1428  parray3=(TObjArray*)parray->Clone();
1429  nPions3=parray3->GetEntries();
1430  }
1431  for(Int_t iTr1=0; iTr1<nKaons; iTr1++){
1432  TLorentzVector* trK=(TLorentzVector*)karray1->At(iTr1);
1433  Double_t chargeK=trK->T();
1434  px[0] = trK->Px();
1435  py[0] = trK->Py();
1436  pz[0] = trK->Pz();
1437  for(Int_t iTr2=0; iTr2<nPions; iTr2++){
1438  TLorentzVector* trPi1=(TLorentzVector*)parray2->At(iTr2);
1439  Double_t chargePi1=trPi1->T();
1440  px[1] = trPi1->Px();
1441  py[1] = trPi1->Py();
1442  pz[1] = trPi1->Pz();
1443  if(chargePi1*chargeK<0){
1444  if(fMeson==kDzero){
1445  FillMEHistos(421,2,tmpRD2,px,py,pz,pdg0);
1446  }else if(fMeson==kDs) {
1447  for(Int_t iTr3=iTr1+1; iTr3<nKaons; iTr3++){
1448  TLorentzVector* trK2=(TLorentzVector*)karray1->At(iTr3);
1449  Double_t chargeK2=trK2->T();
1450  px[2] = trK2->Px();
1451  py[2] = trK2->Py();
1452  pz[2] = trK2->Pz();
1453  Double_t massKK=ComputeInvMassKK(trK,trK2);
1454  Double_t deltaMass=massKK-TDatabasePDG::Instance()->GetParticle(333)->Mass();
1455  Double_t cos1=CosPiKPhiRFrame(trK,trK2,trPi1);
1456  Double_t kincutPiKPhi=TMath::Abs(cos1*cos1*cos1);
1457  Double_t cosPiDsLabFrame=CosPiDsLabFrame(trK,trK2,trPi1);
1458  if(chargeK2*chargeK<0 && TMath::Abs(deltaMass)<fPhiMassCut && kincutPiKPhi>fCutCos3PiKPhiRFrame && cosPiDsLabFrame<fCutCosPiDsLabFrame){
1459  FillMEHistos(431,3,tmpRD3,px,py,pz,pdgs);
1460  }
1461  }
1462  } else{
1463  if(parray3){
1464  for(Int_t iTr3=iTr2+1; iTr3<nPions3; iTr3++){
1465  TLorentzVector* trPi2=(TLorentzVector*)parray3->At(iTr3);
1466  Double_t chargePi2=trPi2->T();
1467  px[2] = trPi2->Px();
1468  py[2] = trPi2->Py();
1469  pz[2] = trPi2->Pz();
1470  if(chargePi2*chargeK<0){
1471  if(fMeson==kDplus) FillMEHistos(411,3,tmpRD3,px,py,pz,pdgp);
1472  }
1473  }
1474  }
1475  }
1476  }else if(chargePi1*chargeK>0){
1477  if(fMeson==kDzero) FillMEHistosLS(421,2,tmpRD2,px,py,pz,pdg0,(Int_t)chargePi1);
1478  }
1479  }
1480  }
1481  delete parray3;
1482  delete parray2;
1483  }
1484  delete karray1;
1485  }
1486  delete tmpRD2;
1487  delete tmpRD3;
1488 }
1489 //_________________________________________________________________
1491 {
1493  if(fDoEventMixing==0) return;
1494  printf("AliAnalysisTaskCombinHF: FinishTaskOutput\n");
1495 
1496  if(fDoEventMixing==1){
1497  for(Int_t i=0; i<fNOfPools; i++){
1498  Int_t nEvents=fEventBuffer[i]->GetEntries();
1499  if(nEvents>1) DoMixingWithPools(i);
1500  }
1501  }else if(fDoEventMixing==2){
1502  DoMixingWithCuts();
1503  }
1504 }
1505 //_________________________________________________________________
1506 Double_t AliAnalysisTaskCombinHF::ComputeInvMassKK(AliAODTrack* tr1, AliAODTrack* tr2) const{
1508  Double_t massK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
1509  Double_t p1=tr1->P();
1510  Double_t p2=tr2->P();
1511  Double_t pxtot=tr1->Px()+tr2->Px();
1512  Double_t pytot=tr1->Py()+tr2->Py();
1513  Double_t pztot=tr1->Pz()+tr2->Pz();
1514  Double_t e1=TMath::Sqrt(massK*massK+p1*p1);
1515  Double_t e2=TMath::Sqrt(massK*massK+p2*p2);
1516  Double_t etot=e1+e2;
1517  Double_t m2=etot*etot-(pxtot*pxtot+pytot*pytot+pztot*pztot);
1518  return TMath::Sqrt(m2);
1519 }
1520 //_________________________________________________________________
1521 Double_t AliAnalysisTaskCombinHF::ComputeInvMassKK(TLorentzVector* tr1, TLorentzVector* tr2) const{
1523  Double_t massK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
1524  Double_t p1=tr1->P();
1525  Double_t p2=tr2->P();
1526  Double_t pxtot=tr1->Px()+tr2->Px();
1527  Double_t pytot=tr1->Py()+tr2->Py();
1528  Double_t pztot=tr1->Pz()+tr2->Pz();
1529  Double_t e1=TMath::Sqrt(massK*massK+p1*p1);
1530  Double_t e2=TMath::Sqrt(massK*massK+p2*p2);
1531  Double_t etot=e1+e2;
1532  Double_t m2=etot*etot-(pxtot*pxtot+pytot*pytot+pztot*pztot);
1533  return TMath::Sqrt(m2);
1534 }
1535 
1536 //----------------------------------------------------------------------
1537 Double_t AliAnalysisTaskCombinHF::CosPiKPhiRFrame(TLorentzVector* dauK1, TLorentzVector* dauK2, TLorentzVector* daupi) const {
1539 
1540  Double_t massK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
1541  Double_t masspi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
1542 
1543  Double_t eK1 = TMath::Sqrt(massK*massK+dauK1->P()*dauK1->P());
1544  Double_t eK2 = TMath::Sqrt(massK*massK+dauK2->P()*dauK2->P());
1545  Double_t epi = TMath::Sqrt(masspi*masspi+daupi->P()*daupi->P());
1546 
1547  Double_t ePhi=eK1+eK2;
1548  Double_t pxPhi=dauK1->Px()+dauK2->Px();
1549  Double_t pyPhi=dauK1->Py()+dauK2->Py();
1550  Double_t pzPhi=dauK1->Pz()+dauK2->Pz();
1551  Double_t bxPhi=pxPhi/ePhi;
1552  Double_t byPhi=pyPhi/ePhi;
1553  Double_t bzPhi=pzPhi/ePhi;
1554 
1555  TVector3 vecK1Phiframe;
1556  TLorentzVector vecK1(dauK1->Px(),dauK1->Py(),dauK1->Pz(),eK1);
1557  vecK1.Boost(-bxPhi,-byPhi,-bzPhi);
1558  vecK1.Boost(vecK1Phiframe);
1559  vecK1Phiframe=vecK1.BoostVector();
1560 
1561  TVector3 vecPiPhiframe;
1562  TLorentzVector vecPi(daupi->Px(),daupi->Py(),daupi->Pz(),epi);
1563  vecPi.Boost(-bxPhi,-byPhi,-bzPhi);
1564  vecPi.Boost(vecPiPhiframe);
1565  vecPiPhiframe=vecPi.BoostVector();
1566 
1567  Double_t innera=vecPiPhiframe.Dot(vecK1Phiframe);
1568  Double_t norm1a=TMath::Sqrt(vecPiPhiframe.Dot(vecPiPhiframe));
1569  Double_t norm2a=TMath::Sqrt(vecK1Phiframe.Dot(vecK1Phiframe));
1570  Double_t cosK1PhiFrame=innera/(norm1a*norm2a);
1571 
1572  return cosK1PhiFrame;
1573 }
1574 
1575 //----------------------------------------------------------------------
1576 Double_t AliAnalysisTaskCombinHF::CosPiDsLabFrame(TLorentzVector* dauK1, TLorentzVector* dauK2, TLorentzVector* daupi) const {
1578 
1579  Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1580  Double_t masspi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
1581 
1582  Double_t pxD=dauK1->Px()+dauK2->Px()+daupi->Px();
1583  Double_t pyD=dauK1->Px()+dauK2->Px()+daupi->Px();
1584  Double_t pzD=dauK1->Px()+dauK2->Px()+daupi->Px();
1585  Double_t pD=TMath::Sqrt(pxD*pxD+pyD*pyD+pzD*pzD);
1586  Double_t eD=TMath::Sqrt(massDs*massDs+pD*pD);
1587  Double_t epi = TMath::Sqrt(masspi*masspi+daupi->P()*daupi->P());
1588  Double_t bxD=pxD/eD;
1589  Double_t byD=pyD/eD;
1590  Double_t bzD=pzD/eD;
1591 
1592  TVector3 piDsframe;
1593  TLorentzVector vecPi(daupi->Px(),daupi->Py(),daupi->Pz(),epi);
1594  vecPi.Boost(-bxD,-byD,-bzD);
1595  vecPi.Boost(piDsframe);
1596  piDsframe=vecPi.BoostVector();
1597 
1598  TVector3 vecDs(pxD,pyD,pzD);
1599 
1600  Double_t inner=vecDs.Dot(piDsframe);
1601  Double_t norm1=TMath::Sqrt(vecDs.Dot(vecDs));
1602  Double_t norm2=TMath::Sqrt(piDsframe.Dot(piDsframe));
1603  Double_t cosPiDsFrame=inner/(norm1*norm2);
1604 
1605  return cosPiDsFrame;
1606 }
1607 
1608 //_________________________________________________________________
1610 {
1612  //
1613  if(fDebug > 1) printf("AliAnalysisTaskCombinHF: Terminate() \n");
1614  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1615  if (!fOutput) {
1616  printf("ERROR: fOutput not available\n");
1617  return;
1618  }
1619  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
1620  if(fHistNEvents){
1621  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1622  }else{
1623  printf("ERROR: fHistNEvents not available\n");
1624  return;
1625  }
1626  return;
1627 }
1628 
Int_t charge
Bool_t IsEventRejectedDueToCentrality() const
Definition: AliRDHFCuts.h:335
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:329
Bool_t IsTrackSelected(AliAODTrack *track)
static Int_t CheckD0Decay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
Bool_t IsEventRejectedDueToNotRecoVertex() const
Definition: AliRDHFCuts.h:320
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:323
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:257
Double_t CosPiDsLabFrame(TLorentzVector *dauK1, TLorentzVector *dauK2, TLorentzVector *daupi) const
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.
Bool_t IsEventRejectedDueToBadTrackVertex() const
Definition: AliRDHFCuts.h:344
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:332
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 IsEventRejectedDueToBadPrimaryVertex() const
Definition: AliRDHFCuts.h:353
Bool_t IsEventRejectedDuePhysicsSelection() const
Definition: AliRDHFCuts.h:350
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:317
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:312
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:280
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