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