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