]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add event mixing option for background estimation in low pt D-meson task
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskCombinHF.cxx
CommitLineData
811e9cbd 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
23/////////////////////////////////////////////////////////////
24
25#include <TList.h>
26#include <TH1F.h>
ebc460de 27#include <TH2F.h>
811e9cbd 28#include <TH3F.h>
29#include <THnSparse.h>
30
31#include "AliAnalysisManager.h"
32#include "AliInputEventHandler.h"
33#include "AliPIDResponse.h"
34#include "AliAODHandler.h"
35#include "AliAODEvent.h"
36#include "AliAODMCParticle.h"
37#include "AliAODMCHeader.h"
38#include "AliAODVertex.h"
39#include "AliAODTrack.h"
ebc460de 40#include "AliVertexingHFUtils.h"
811e9cbd 41#include "AliAnalysisTaskCombinHF.h"
42
43ClassImp(AliAnalysisTaskCombinHF)
44
45
46//________________________________________________________________________
47AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF():
48 AliAnalysisTaskSE(),
f2221737 49 fOutput(0x0),
811e9cbd 50 fHistNEvents(0x0),
51 fHistTrackStatus(0x0),
ebc460de 52 fHistCheckOrigin(0x0),
53 fHistCheckOriginSel(0x0),
54 fHistCheckDecChan(0x0),
55 fHistCheckDecChanAcc(0x0),
56 fPtVsYGen(0x0),
f31ef968 57 fPtVsYGenLargeAcc(0x0),
ebc460de 58 fPtVsYGenLimAcc(0x0),
59 fPtVsYGenAcc(0x0),
60 fPtVsYReco(0x0),
811e9cbd 61 fMassVsPtVsY(0x0),
62 fMassVsPtVsYRot(0x0),
63 fMassVsPtVsYLSpp(0x0),
64 fMassVsPtVsYLSmm(0x0),
65 fMassVsPtVsYSig(0x0),
66 fMassVsPtVsYRefl(0x0),
67 fMassVsPtVsYBkg(0x0),
68 fNSelected(0x0),
69 fNormRotated(0x0),
70 fDeltaMass(0x0),
71 fDeltaMassFullAnalysis(0x0),
93f0ffc1 72 fMassVsPtVsYME(0x0),
811e9cbd 73 fFilterMask(BIT(4)),
74 fTrackCutsAll(0x0),
75 fTrackCutsPion(0x0),
76 fTrackCutsKaon(0x0),
77 fPidHF(new AliAODPidHF()),
78 fAnalysisCuts(0x0),
f2221737 79 fMinMass(1.720),
811e9cbd 80 fMaxMass(2.150),
f2221737 81 fMaxPt(10.),
adf098d9 82 fPtBinWidth(0.5),
ebc460de 83 fEtaAccCut(0.9),
84 fPtAccCut(0.1),
811e9cbd 85 fNRotations(9),
86 fMinAngleForRot(5*TMath::Pi()/6),
87 fMaxAngleForRot(7*TMath::Pi()/6),
88 fMinAngleForRot3(2*TMath::Pi()/6),
89 fMaxAngleForRot3(4*TMath::Pi()/6),
90 fCounter(0x0),
91 fMeson(kDzero),
92 fReadMC(kFALSE),
ebc460de 93 fPromptFeeddown(kPrompt),
94 fGoUpToQuark(kTRUE),
a25d1750 95 fFullAnalysis(0),
f2221737 96 fPIDstrategy(knSigma),
97 fmaxPforIDPion(0.8),
98 fmaxPforIDKaon(2.),
99 fKeepNegID(kFALSE),
100 fPIDselCaseZero(0),
101 fBayesThresKaon(0.4),
93f0ffc1 102 fBayesThresPion(0.4),
103 fDoEventMixing(kTRUE),
104 fMaxNumberOfEventsForMixing(20),
105 fMaxzVertDistForMix(20.),
106 fMaxMultDiffForMix(9999.),
107 fUsePoolsZ(kFALSE),
108 fUsePoolsM(kFALSE),
109 fNzVertPools(1),
110 fzVertPoolLims(0x0),
111 fNMultPools(1),
112 fMultPoolLims(0x0),
113 fEventBuffer(0x0),
114 fVtxZ(0),
115 fMultiplicity(0),
116 fKaonTracks(0x0),
117 fPionTracks(0x0)
811e9cbd 118{
119 // default constructor
120}
121
122//________________________________________________________________________
123AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analysiscuts):
124 AliAnalysisTaskSE("DmesonCombin"),
f2221737 125 fOutput(0x0),
811e9cbd 126 fHistNEvents(0x0),
127 fHistTrackStatus(0x0),
ebc460de 128 fHistCheckOrigin(0x0),
129 fHistCheckOriginSel(0x0),
130 fHistCheckDecChan(0x0),
131 fHistCheckDecChanAcc(0x0),
132 fPtVsYGen(0x0),
f31ef968 133 fPtVsYGenLargeAcc(0x0),
ebc460de 134 fPtVsYGenLimAcc(0x0),
135 fPtVsYGenAcc(0x0),
136 fPtVsYReco(0x0),
811e9cbd 137 fMassVsPtVsY(0x0),
138 fMassVsPtVsYRot(0x0),
139 fMassVsPtVsYLSpp(0x0),
140 fMassVsPtVsYLSmm(0x0),
141 fMassVsPtVsYSig(0x0),
142 fMassVsPtVsYRefl(0x0),
143 fMassVsPtVsYBkg(0x0),
144 fNSelected(0x0),
145 fNormRotated(0x0),
146 fDeltaMass(0x0),
147 fDeltaMassFullAnalysis(0x0),
93f0ffc1 148 fMassVsPtVsYME(0x0),
811e9cbd 149 fFilterMask(BIT(4)),
150 fTrackCutsAll(0x0),
151 fTrackCutsPion(0x0),
152 fTrackCutsKaon(0x0),
153 fPidHF(new AliAODPidHF()),
154 fAnalysisCuts(analysiscuts),
f2221737 155 fMinMass(1.720),
811e9cbd 156 fMaxMass(2.150),
f2221737 157 fMaxPt(10.),
adf098d9 158 fPtBinWidth(0.5),
ebc460de 159 fEtaAccCut(0.9),
160 fPtAccCut(0.1),
811e9cbd 161 fNRotations(9),
162 fMinAngleForRot(5*TMath::Pi()/6),
163 fMaxAngleForRot(7*TMath::Pi()/6),
164 fMinAngleForRot3(2*TMath::Pi()/6),
165 fMaxAngleForRot3(4*TMath::Pi()/6),
166 fCounter(0x0),
167 fMeson(meson),
168 fReadMC(kFALSE),
ebc460de 169 fPromptFeeddown(1),
170 fGoUpToQuark(kTRUE),
a25d1750 171 fFullAnalysis(0),
f2221737 172 fPIDstrategy(knSigma),
173 fmaxPforIDPion(0.8),
174 fmaxPforIDKaon(2.),
175 fKeepNegID(kFALSE),
176 fPIDselCaseZero(0),
177 fBayesThresKaon(0.4),
93f0ffc1 178 fBayesThresPion(0.4),
179 fDoEventMixing(kTRUE),
180 fMaxNumberOfEventsForMixing(20),
181 fMaxzVertDistForMix(20.),
182 fMaxMultDiffForMix(9999.),
183 fUsePoolsZ(kFALSE),
184 fUsePoolsM(kFALSE),
185 fNzVertPools(1),
186 fzVertPoolLims(0x0),
187 fNMultPools(1),
188 fMultPoolLims(0x0),
189 fEventBuffer(0x0),
190 fVtxZ(0),
191 fMultiplicity(0),
192 fKaonTracks(0x0),
193 fPionTracks(0x0)
811e9cbd 194{
195 // standard constructor
93f0ffc1 196
811e9cbd 197 DefineOutput(1,TList::Class()); //My private output
198 DefineOutput(2,AliNormalizationCounter::Class());
f2221737 199}
811e9cbd 200
201//________________________________________________________________________
202AliAnalysisTaskCombinHF::~AliAnalysisTaskCombinHF()
203{
204 //
205 // Destructor
206 //
93f0ffc1 207 if(!fOutput->IsOwner()){
208 delete fHistNEvents;
209 delete fHistTrackStatus;
210 delete fHistCheckOrigin;
211 delete fHistCheckOriginSel;
212 delete fHistCheckDecChan;
213 delete fHistCheckDecChanAcc;
214 delete fPtVsYGen;
215 delete fPtVsYGenLargeAcc;
216 delete fPtVsYGenLimAcc;
217 delete fPtVsYGenAcc;
218 delete fPtVsYReco;
219 delete fMassVsPtVsY;
220 delete fMassVsPtVsYLSpp;
221 delete fMassVsPtVsYLSmm;
222 delete fMassVsPtVsYRot;
223 delete fMassVsPtVsYSig;
224 delete fMassVsPtVsYRefl;
225 delete fMassVsPtVsYBkg;
226 delete fNSelected;
227 delete fNormRotated;
228 delete fDeltaMass;
229 delete fDeltaMassFullAnalysis;
230 delete fMassVsPtVsYME;
231 }
232
811e9cbd 233 delete fOutput;
811e9cbd 234 delete fCounter;
235 delete fTrackCutsAll;
236 delete fTrackCutsPion;
237 delete fTrackCutsKaon;
238 delete fPidHF;
239 delete fAnalysisCuts;
93f0ffc1 240 fKaonTracks->Delete();
241 fPionTracks->Delete();
242 delete fKaonTracks;
243 delete fPionTracks;
244 delete fEventBuffer;
811e9cbd 245}
246
93f0ffc1 247//________________________________________________________________________
248void AliAnalysisTaskCombinHF::ConfigureZVertPools(Int_t nPools, Double_t* zVertLimits)
249{
250 // sets the pools for event mizing in zvertex
251 if(fzVertPoolLims) delete [] fzVertPoolLims;
252 fNzVertPools=nPools;
253 fzVertPoolLims = new Double_t[fNzVertPools+1];
254 for(Int_t ib=0; ib<nPools+1; ib++) fzVertPoolLims[ib]=zVertLimits[ib];
255 fUsePoolsZ=kTRUE;
256 return;
257}
258//________________________________________________________________________
259void AliAnalysisTaskCombinHF::ConfigureMultiplicityPools(Int_t nPools, Double_t* multLimits)
260{
261 // sets the pools for event mizing in zvertex
262 if(fMultPoolLims) delete [] fMultPoolLims;
263 fNMultPools=nPools;
264 fMultPoolLims = new Double_t[fNMultPools+1];
265 for(Int_t ib=0; ib<nPools+1; ib++) fMultPoolLims[ib]=multLimits[ib];
266 fUsePoolsM=kTRUE;
267 return;
268}
811e9cbd 269//________________________________________________________________________
270void AliAnalysisTaskCombinHF::UserCreateOutputObjects()
271{
272 // Create the output container
273 //
274 if(fDebug > 1) printf("AnalysisTaskCombinHF::UserCreateOutputObjects() \n");
f2221737 275
811e9cbd 276 fOutput = new TList();
277 fOutput->SetOwner();
278 fOutput->SetName("OutputHistos");
f2221737 279
811e9cbd 280 fHistNEvents = new TH1F("hNEvents", "number of events ",8,-0.5,7.5);
281 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
282 fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
283 fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
284 fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to not reco vertex");
285 fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected for contr vertex");
286 fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for vertex out of accept");
287 fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for pileup events");
288 fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
f2221737 289
811e9cbd 290 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
291 fHistNEvents->Sumw2();
292 fHistNEvents->SetMinimum(0);
293 fOutput->Add(fHistNEvents);
f2221737 294
811e9cbd 295 fHistTrackStatus = new TH1F("hTrackStatus", "",8,-0.5,7.5);
296 fHistTrackStatus->GetXaxis()->SetBinLabel(1,"Not OK");
297 fHistTrackStatus->GetXaxis()->SetBinLabel(2,"Track OK");
298 fHistTrackStatus->GetXaxis()->SetBinLabel(3,"Kaon, Not OK");
299 fHistTrackStatus->GetXaxis()->SetBinLabel(4,"Kaon OK");
300 fHistTrackStatus->GetXaxis()->SetBinLabel(5,"Pion, Not OK");
301 fHistTrackStatus->GetXaxis()->SetBinLabel(6,"Pion OK");
302 fHistTrackStatus->GetXaxis()->SetBinLabel(7,"Kaon||Pion, Not OK");
303 fHistTrackStatus->GetXaxis()->SetBinLabel(8,"Kaon||Pion OK");
2b6a376c 304 fHistTrackStatus->GetXaxis()->SetNdivisions(1,kFALSE);
811e9cbd 305 fHistTrackStatus->Sumw2();
306 fHistTrackStatus->SetMinimum(0);
307 fOutput->Add(fHistTrackStatus);
f2221737 308
adf098d9 309 Int_t nPtBins = (Int_t)(fMaxPt/fPtBinWidth+0.001);
310 Double_t maxPt=fPtBinWidth*nPtBins;
311
ebc460de 312 if(fReadMC){
f2221737 313
ebc460de 314 fHistCheckOrigin=new TH1F("hCheckOrigin","",7,-1.5,5.5);
315 fHistCheckOrigin->Sumw2();
316 fHistCheckOrigin->SetMinimum(0);
317 fOutput->Add(fHistCheckOrigin);
f2221737 318
ebc460de 319 fHistCheckOriginSel=new TH1F("hCheckOriginSel","",7,-1.5,5.5);
320 fHistCheckOriginSel->Sumw2();
321 fHistCheckOriginSel->SetMinimum(0);
322 fOutput->Add(fHistCheckOriginSel);
f2221737 323
ebc460de 324 fHistCheckDecChan=new TH1F("hCheckDecChan","",7,-2.5,4.5);
325 fHistCheckDecChan->Sumw2();
326 fHistCheckDecChan->SetMinimum(0);
327 fOutput->Add(fHistCheckDecChan);
f2221737 328
ebc460de 329 fHistCheckDecChanAcc=new TH1F("hCheckDecChanAcc","",7,-2.5,4.5);
330 fHistCheckDecChanAcc->Sumw2();
331 fHistCheckDecChanAcc->SetMinimum(0);
332 fOutput->Add(fHistCheckDecChanAcc);
f2221737 333
adf098d9 334 fPtVsYGen= new TH2F("hPtVsYGen","",nPtBins,0.,maxPt,20,-1.,1.);
ebc460de 335 fPtVsYGen->Sumw2();
336 fPtVsYGen->SetMinimum(0);
337 fOutput->Add(fPtVsYGen);
f2221737 338
adf098d9 339 fPtVsYGenLargeAcc= new TH2F("hPtVsYGenLargeAcc","",nPtBins,0.,maxPt,20,-1.,1.);
f31ef968 340 fPtVsYGenLargeAcc->Sumw2();
341 fPtVsYGenLargeAcc->SetMinimum(0);
342 fOutput->Add(fPtVsYGenLargeAcc);
343
adf098d9 344 fPtVsYGenLimAcc= new TH2F("hPtVsYGenLimAcc","",nPtBins,0.,maxPt,20,-1.,1.);
ebc460de 345 fPtVsYGenLimAcc->Sumw2();
346 fPtVsYGenLimAcc->SetMinimum(0);
347 fOutput->Add(fPtVsYGenLimAcc);
f2221737 348
adf098d9 349 fPtVsYGenAcc= new TH2F("hPtVsYGenAcc","",nPtBins,0.,maxPt,20,-1.,1.);
ebc460de 350 fPtVsYGenAcc->Sumw2();
351 fPtVsYGenAcc->SetMinimum(0);
352 fOutput->Add(fPtVsYGenAcc);
f2221737 353
adf098d9 354 fPtVsYReco= new TH2F("hPtVsYReco","",nPtBins,0.,maxPt,20,-1.,1.);
ebc460de 355 fPtVsYReco->Sumw2();
356 fPtVsYReco->SetMinimum(0);
357 fOutput->Add(fPtVsYReco);
358 }
f2221737 359
360
361 Int_t nMassBins=fMaxMass*1000.-fMinMass*1000.;
811e9cbd 362 Double_t maxm=fMinMass+nMassBins*0.001;
adf098d9 363 fMassVsPtVsY=new TH3F("hMassVsPtVsY","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
811e9cbd 364 fMassVsPtVsY->Sumw2();
365 fMassVsPtVsY->SetMinimum(0);
366 fOutput->Add(fMassVsPtVsY);
f2221737 367
adf098d9 368 fMassVsPtVsYRot=new TH3F("hMassVsPtVsYRot","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
811e9cbd 369 fMassVsPtVsYRot->Sumw2();
370 fMassVsPtVsYRot->SetMinimum(0);
371 fOutput->Add(fMassVsPtVsYRot);
f2221737 372
811e9cbd 373 if(fMeson==kDzero){
adf098d9 374 fMassVsPtVsYLSpp=new TH3F("hMassVsPtVsYLSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
811e9cbd 375 fMassVsPtVsYLSpp->Sumw2();
376 fMassVsPtVsYLSpp->SetMinimum(0);
377 fOutput->Add(fMassVsPtVsYLSpp);
adf098d9 378 fMassVsPtVsYLSmm=new TH3F("hMassVsPtVsYLSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
811e9cbd 379 fMassVsPtVsYLSmm->Sumw2();
380 fMassVsPtVsYLSmm->SetMinimum(0);
381 fOutput->Add(fMassVsPtVsYLSmm);
382 }
f2221737 383
adf098d9 384 fMassVsPtVsYSig=new TH3F("hMassVsPtVsYSig","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
811e9cbd 385 fMassVsPtVsYSig->Sumw2();
386 fMassVsPtVsYSig->SetMinimum(0);
387 fOutput->Add(fMassVsPtVsYSig);
f2221737 388
adf098d9 389 fMassVsPtVsYRefl=new TH3F("hMassVsPtVsYRefl","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
811e9cbd 390 fMassVsPtVsYRefl->Sumw2();
391 fMassVsPtVsYRefl->SetMinimum(0);
392 fOutput->Add(fMassVsPtVsYRefl);
f2221737 393
adf098d9 394 fMassVsPtVsYBkg=new TH3F("hMassVsPtVsYBkg","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
811e9cbd 395 fMassVsPtVsYBkg->Sumw2();
396 fMassVsPtVsYBkg->SetMinimum(0);
397 fOutput->Add(fMassVsPtVsYBkg);
f2221737 398
811e9cbd 399 fNSelected=new TH1F("hNSelected","",100,-0.5,99.5);
400 fNSelected->Sumw2();
401 fNSelected->SetMinimum(0);
402 fOutput->Add(fNSelected);
f2221737 403
811e9cbd 404 fNormRotated=new TH1F("hNormRotated","",11,-0.5,10.5);
405 fNormRotated->Sumw2();
406 fNormRotated->SetMinimum(0);
407 fOutput->Add(fNormRotated);
f2221737 408
811e9cbd 409 fDeltaMass=new TH1F("hDeltaMass","",100,-0.4,0.4);
410 fDeltaMass->Sumw2();
411 fDeltaMass->SetMinimum(0);
412 fOutput->Add(fDeltaMass);
f2221737 413
811e9cbd 414 Int_t binSparseDMassRot[5]={nMassBins,100,24,40,20};
415 Double_t edgeLowSparseDMassRot[5]={fMinMass,-0.4,0.,-4.,0};
416 Double_t edgeHighSparseDMassRot[5]={maxm,0.4,12.,4.,3.14};
417 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);
418 fOutput->Add(fDeltaMassFullAnalysis);
f2221737 419
93f0ffc1 420 fMassVsPtVsYME=new TH3F("hMassVsPtVsYME","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
421 fMassVsPtVsYME->Sumw2();
422 fMassVsPtVsYME->SetMinimum(0);
423 fOutput->Add(fMassVsPtVsYME);
424
425
811e9cbd 426 //Counter for Normalization
427 fCounter = new AliNormalizationCounter("NormalizationCounter");
428 fCounter->Init();
f2221737 429
93f0ffc1 430 fKaonTracks = new TObjArray();
431 fPionTracks=new TObjArray();
432 fKaonTracks->SetOwner();
433 fPionTracks->SetOwner();
434 fEventBuffer = new TTree("EventBuffer", "Temporary buffer for event mixing");
435 fEventBuffer->Branch("zVertex", &fVtxZ);
436 fEventBuffer->Branch("multiplicity", &fMultiplicity);
437 fEventBuffer->Branch("karray", "TObjArray", &fKaonTracks);
438 fEventBuffer->Branch("parray", "TObjArray", &fPionTracks);
439
f2221737 440 PostData(1,fOutput);
441 PostData(2,fCounter);
811e9cbd 442}
443
444//________________________________________________________________________
445void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){
446 //Build the 3-track combinatorics (+-+ and -+-) for D+->Kpipi decays
f2221737 447
811e9cbd 448 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
449 if(!aod && AODEvent() && IsStandardAOD()) {
f2221737 450 // In case there is an AOD handler writing a standard AOD, use the AOD
451 // event in memory rather than the input (ESD) event.
811e9cbd 452 aod = dynamic_cast<AliAODEvent*> (AODEvent());
453 }
454 if(!aod){
455 printf("AliAnalysisTaskCombinHF::UserExec: AOD not found!\n");
456 return;
457 }
f2221737 458
459 // fix for temporary bug in ESDfilter
811e9cbd 460 // the AODs with null vertex pointer didn't pass the PhysSel
461 if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
462 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
463 AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
464 AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
465 fPidHF->SetPidResponse(pidResp);
466
f2221737 467
811e9cbd 468 fHistNEvents->Fill(0); // count event
469 // Post the data already here
470 PostData(1,fOutput);
f2221737 471
811e9cbd 472 fCounter->StoreEvent(aod,fAnalysisCuts,fReadMC);
f2221737 473
811e9cbd 474 Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
475 if(fAnalysisCuts->IsEventRejectedDueToTrigger())fHistNEvents->Fill(2);
476 if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(3);
477 if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(4);
478 if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(5);
479 // if(fAnalysisCuts->IsEventRejectedDueToPileup())fHistNEvents->Fill(6);
f2221737 480 if(fAnalysisCuts->IsEventRejectedDueToCentrality()) fHistNEvents->Fill(7);
481
811e9cbd 482 if(!isEvSel)return;
f2221737 483
811e9cbd 484 fHistNEvents->Fill(1);
485
486 TClonesArray *arrayMC=0;
487 AliAODMCHeader *mcHeader=0;
488 if(fReadMC){
489 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
490 if(!arrayMC) {
491 printf("AliAnalysisTaskCombinHF::UserExec: MC particles branch not found!\n");
492 return;
493 }
494
495 // load MC header
496 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
497 if(!mcHeader) {
498 printf("AliAnalysisTaskCombinHF::UserExec: MC header branch not found!\n");
499 return;
500 }
ebc460de 501 FillGenHistos(arrayMC);
811e9cbd 502 }
f2221737 503
504
811e9cbd 505 Int_t ntracks=aod->GetNTracks();
93f0ffc1 506 fVtxZ = aod->GetPrimaryVertex()->GetZ();
507 fMultiplicity = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
508
811e9cbd 509 // select and flag tracks
510 UChar_t* status = new UChar_t[ntracks];
511 for(Int_t iTr=0; iTr<ntracks; iTr++){
512 status[iTr]=0;
513 AliAODTrack* track=aod->GetTrack(iTr);
514 if(IsTrackSelected(track)) status[iTr]+=1;
f2221737 515
516 // PID
517 if (fPIDstrategy == knSigma) {
518 // nsigma PID
519 if(IsKaon(track)) status[iTr]+=2;
520 if(IsPion(track)) status[iTr]+=4;
521 }
522 else if (fPIDstrategy == kBayesianMaxProb || fPIDstrategy == kBayesianThres) {
523 // Bayesian PID
524 Double_t *weights = new Double_t[AliPID::kSPECIES];
525 fPidHF->GetPidCombined()->ComputeProbabilities(track, fPidHF->GetPidResponse(), weights);
526 if (fPIDstrategy == kBayesianMaxProb) {
527 if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kKaon]) status[iTr] += 2;
528 if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kPion]) status[iTr] += 4;
529 }
530 if (fPIDstrategy == kBayesianThres) {
531 if (weights[AliPID::kKaon] > fBayesThresKaon) status[iTr] += 2;
532 if (weights[AliPID::kPion] > fBayesThresPion) status[iTr] += 4;
533 }
534 delete[] weights;
535 }
536
811e9cbd 537 fHistTrackStatus->Fill(status[iTr]);
538 }
f2221737 539
811e9cbd 540 // build the combinatorics
541 Int_t nSelected=0;
542 Int_t nFiltered=0;
04731395 543 Double_t dummypos[3]={0.,0.,0.};
544 AliAODVertex* v2=new AliAODVertex(dummypos,999.,-1,2);
545 AliAODVertex* v3=new AliAODVertex(dummypos,999.,-1,3);
811e9cbd 546 // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
547 Double_t d02[2]={0.,0.};
548 Double_t d03[3]={0.,0.,0.};
549 AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
550 AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
551 UInt_t pdg0[2]={321,211};
552 UInt_t pdgp[3]={321,211,211};
553 // UInt_t pdgs[3]={321,321,211};
554 Double_t tmpp[3];
555 Double_t px[3],py[3],pz[3];
556 Int_t dgLabels[3];
93f0ffc1 557 fKaonTracks->Delete();
558 fPionTracks->Delete();
559
811e9cbd 560 for(Int_t iTr1=0; iTr1<ntracks; iTr1++){
561 AliAODTrack* trK=aod->GetTrack(iTr1);
562 if((status[iTr1] & 1)==0) continue;
93f0ffc1 563 if(fDoEventMixing){
564 if(status[iTr1] & 2) fKaonTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
565 if(status[iTr1] & 4) fPionTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
566 }
811e9cbd 567 if((status[iTr1] & 2)==0) continue;
568 Int_t chargeK=trK->Charge();
569 trK->GetPxPyPz(tmpp);
f2221737 570 px[0] = tmpp[0];
571 py[0] = tmpp[1];
572 pz[0] = tmpp[2];
811e9cbd 573 dgLabels[0]=trK->GetLabel();
574 for(Int_t iTr2=0; iTr2<ntracks; iTr2++){
575 if((status[iTr2] & 1)==0) continue;
576 if((status[iTr2] & 4)==0) continue;
577 if(iTr1==iTr2) continue;
578 AliAODTrack* trPi1=aod->GetTrack(iTr2);
579 Int_t chargePi1=trPi1->Charge();
580 trPi1->GetPxPyPz(tmpp);
f2221737 581 px[1] = tmpp[0];
582 py[1] = tmpp[1];
583 pz[1] = tmpp[2];
811e9cbd 584 dgLabels[1]=trPi1->GetLabel();
585 if(chargePi1==chargeK){
f2221737 586 if(fMeson==kDzero) FillLSHistos(421,2,tmpRD2,px,py,pz,pdg0,chargePi1);
587 continue;
811e9cbd 588 }
589 if(fMeson==kDzero){
f2221737 590 nFiltered++;
591 v2->AddDaughter(trK);
592 v2->AddDaughter(trPi1);
593 tmpRD2->SetSecondaryVtx(v2);
594 Bool_t ok=FillHistos(421,2,tmpRD2,px,py,pz,pdg0,arrayMC,dgLabels);
595 v2->RemoveDaughters();
596 if(ok) nSelected++;
811e9cbd 597 }else{
f2221737 598 for(Int_t iTr3=iTr2+1; iTr3<ntracks; iTr3++){
599 if((status[iTr3] & 1)==0) continue;
600 if((status[iTr3] & 4)==0) continue;
601 if(iTr1==iTr3) continue;
602 AliAODTrack* trPi2=aod->GetTrack(iTr3);
603 Int_t chargePi2=trPi2->Charge();
604 if(chargePi2==chargeK) continue;
605 nFiltered++;
606 trPi2->GetPxPyPz(tmpp);
607 px[2] = tmpp[0];
608 py[2] = tmpp[1];
609 pz[2] = tmpp[2];
610 dgLabels[2]=trPi2->GetLabel();
611 v3->AddDaughter(trK);
612 v3->AddDaughter(trPi1);
613 v3->AddDaughter(trPi2);
614 tmpRD3->SetSecondaryVtx(v3);
615 Bool_t ok=FillHistos(411,3,tmpRD3,px,py,pz,pdgp,arrayMC,dgLabels);
616 v3->RemoveDaughters();
617 if(ok) nSelected++;
618 }
811e9cbd 619 }
620 }
621 }
f2221737 622
811e9cbd 623 delete [] status;
624 delete v2;
625 delete v3;
626 delete tmpRD2;
627 delete tmpRD3;
f2221737 628
811e9cbd 629 fNSelected->Fill(nSelected);
f2221737 630
811e9cbd 631 fCounter->StoreCandidates(aod,nFiltered,kTRUE);
632 fCounter->StoreCandidates(aod,nSelected,kFALSE);
93f0ffc1 633
634 if(fDoEventMixing) fEventBuffer->Fill();
f2221737 635 PostData(1,fOutput);
636 PostData(2,fCounter);
637
811e9cbd 638 return;
639}
640
641//________________________________________________________________________
642void 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){
643 // Fill histos for LS candidates
f2221737 644
811e9cbd 645 tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
646 Double_t pt = tmpRD->Pt();
647 Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
648 if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
649 Double_t rapid = tmpRD->Y(pdgD);
650 if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
651 if(charge>0) fMassVsPtVsYLSpp->Fill(TMath::Sqrt(minv2),pt,rapid);
652 else fMassVsPtVsYLSmm->Fill(TMath::Sqrt(minv2),pt,rapid);
653 }
654 }
655 return;
656}
657
ebc460de 658//________________________________________________________________________
659void AliAnalysisTaskCombinHF::FillGenHistos(TClonesArray* arrayMC){
660 // Fill histos with generated quantities
661 Int_t totPart=arrayMC->GetEntriesFast();
662 Int_t thePDG=411;
663 Int_t nProng=3;
664 if(fMeson==kDzero){
665 thePDG=421;
666 nProng=2;
667 }
ebc460de 668 for(Int_t ip=0; ip<totPart; ip++){
669 AliAODMCParticle *part = (AliAODMCParticle*)arrayMC->At(ip);
670 if(TMath::Abs(part->GetPdgCode())==thePDG){
671 Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,part,fGoUpToQuark);
672 fHistCheckOrigin->Fill(orig);
673 if(fPromptFeeddown==kFeeddown && orig!=5) continue;
674 else if(fPromptFeeddown==kPrompt && orig!=4) continue;
675 else if(fPromptFeeddown==kBoth && orig<4) continue;
676 fHistCheckOriginSel->Fill(orig);
677 Int_t deca=0;
678 Bool_t isGoodDecay=kFALSE;
f1f4f53a 679 Int_t labDau[4]={-1,-1,-1,-1};
ebc460de 680 if(fMeson==kDzero){
f2221737 681 deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,part,labDau);
682 if(part->GetNDaughters()!=2) continue;
683 if(deca==1) isGoodDecay=kTRUE;
684 }else if(fMeson==kDplus){
685 deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,part,labDau);
686 if(deca>0) isGoodDecay=kTRUE;
ebc460de 687 }
ebc460de 688 fHistCheckDecChan->Fill(deca);
f1f4f53a 689 if(labDau[0]==-1){
f2221737 690 // printf(Form("Meson %d Label of daughters not filled correctly -- %d\n",fMeson,isGoodDecay));
691 continue; //protection against unfilled array of labels
f1f4f53a 692 }
693 Bool_t isInAcc=CheckAcceptance(arrayMC,nProng,labDau);
ebc460de 694 if(isInAcc) fHistCheckDecChanAcc->Fill(deca);
695 if(isGoodDecay){
f2221737 696 Double_t ptgen=part->Pt();
697 Double_t ygen=part->Y();
ebc460de 698 if(fAnalysisCuts->IsInFiducialAcceptance(ptgen,ygen)){
699 fPtVsYGen->Fill(ptgen,ygen);
700 if(TMath::Abs(ygen)<0.5) fPtVsYGenLimAcc->Fill(ptgen,ygen);
701 if(isInAcc) fPtVsYGenAcc->Fill(ptgen,ygen);
702 }
f31ef968 703 if(TMath::Abs(ygen)<0.9) fPtVsYGenLargeAcc->Fill(ptgen,ygen);
ebc460de 704 }
705 }
706 }
707}
708
811e9cbd 709//________________________________________________________________________
710Bool_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){
711 // Fill histos for candidates with proper charge sign
f2221737 712
811e9cbd 713 Bool_t accept=kFALSE;
f2221737 714
811e9cbd 715 tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
716 Double_t pt = tmpRD->Pt();
717 Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
718 Double_t mass=TMath::Sqrt(minv2);
f2221737 719
811e9cbd 720 if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
721 Double_t rapid = tmpRD->Y(pdgD);
722 if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
723 fMassVsPtVsY->Fill(mass,pt,rapid);
724 accept=kTRUE;
725 if(fReadMC){
f2221737 726 Int_t signPdg[3]={0,0,0};
727 for(Int_t iii=0; iii<nProngs; iii++) signPdg[iii]=pdgdau[iii];
728 Int_t labD = tmpRD->MatchToMC(pdgD,arrayMC,nProngs,signPdg);
729 if(labD>=0){
730 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(TMath::Abs(dgLabels[0])));
731 if(part){
732 Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,part,fGoUpToQuark);
733 if((fPromptFeeddown==kFeeddown && orig==5)|| (fPromptFeeddown==kPrompt && orig==4) || (fPromptFeeddown==kBoth && orig>=4)) {
734
735 Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
736 if(pdgCode==321){
737 fMassVsPtVsYSig->Fill(mass,pt,rapid);
738 AliAODMCParticle* dmes = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labD));
739 fPtVsYReco->Fill(dmes->Pt(),dmes->Y());
740 }else{
741 fMassVsPtVsYRefl->Fill(mass,pt,rapid);
742 }
743 }
744 }
745 }else{
746 fMassVsPtVsYBkg->Fill(mass,pt,rapid);
747 }
811e9cbd 748 }
749 }
750 }
f2221737 751
811e9cbd 752 Int_t nRotated=0;
753 Double_t massRot=0;// calculated later only if candidate is acceptable
754 Double_t angleProngXY;
755 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])));
756 else {
757 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])));
758 }
759 Double_t ptOrig=pt;
f2221737 760
761
811e9cbd 762 Double_t rotStep=(fMaxAngleForRot-fMinAngleForRot)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
763 Double_t rotStep3=(fMaxAngleForRot3-fMinAngleForRot3)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
f2221737 764
811e9cbd 765 for(Int_t irot=0; irot<fNRotations; irot++){
766 Double_t phirot=fMinAngleForRot+rotStep*irot;
767 Double_t tmpx=px[0];
768 Double_t tmpy=py[0];
769 Double_t tmpx2=px[2];
770 Double_t tmpy2=py[2];
771 px[0]=tmpx*TMath::Cos(phirot)-tmpy*TMath::Sin(phirot);
772 py[0]=tmpx*TMath::Sin(phirot)+tmpy*TMath::Cos(phirot);
773 if(pdgD==411 || pdgD==431){
774 Double_t phirot2=fMinAngleForRot3+rotStep3*irot;
775 px[2]=tmpx*TMath::Cos(phirot2)-tmpy*TMath::Sin(phirot2);
776 py[2]=tmpx*TMath::Sin(phirot2)+tmpy*TMath::Cos(phirot2);
777 }
778 tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
779 pt = tmpRD->Pt();
780 minv2 = tmpRD->InvMass2(nProngs,pdgdau);
781 if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
782 Double_t rapid = tmpRD->Y(pdgD);
783 if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
f2221737 784 massRot=TMath::Sqrt(minv2);
785 fMassVsPtVsYRot->Fill(massRot,pt,rapid);
786 nRotated++;
787 fDeltaMass->Fill(massRot-mass);
788 if(fFullAnalysis){
789 Double_t pointRot[5]={mass,massRot-mass,ptOrig,pt-ptOrig,angleProngXY};
790 fDeltaMassFullAnalysis->Fill(pointRot);
791 }
811e9cbd 792 }
793 }
794 px[0]=tmpx;
795 py[0]=tmpy;
796 if(pdgD==411 || pdgD==431){
797 px[2]=tmpx2;
798 py[2]=tmpy2;
799 }
800 }
801 fNormRotated->Fill(nRotated);
f2221737 802
811e9cbd 803 return accept;
f2221737 804
811e9cbd 805}
93f0ffc1 806//________________________________________________________________________
807void AliAnalysisTaskCombinHF::FillMEHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau){
808 // Fill histos for candidates in MixedEvents
809
810 tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
811 Double_t pt = tmpRD->Pt();
812 Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
813 Double_t mass=TMath::Sqrt(minv2);
814
815 if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
816 Double_t rapid = tmpRD->Y(pdgD);
817 if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
818 fMassVsPtVsYME->Fill(mass,pt,rapid);
819 }
820 }
821 return;
822}
823
811e9cbd 824//________________________________________________________________________
825Bool_t AliAnalysisTaskCombinHF::IsTrackSelected(AliAODTrack* track){
826 // track selection cuts
f2221737 827
811e9cbd 828 if(track->Charge()==0) return kFALSE;
a25d1750 829 if(track->GetID()<0&&!fKeepNegID)return kFALSE;
811e9cbd 830 if(!(track->TestFilterMask(fFilterMask))) return kFALSE;
f2221737 831 if(!SelectAODTrack(track,fTrackCutsAll)) return kFALSE;
811e9cbd 832 return kTRUE;
833}
f2221737 834
811e9cbd 835//________________________________________________________________________
836Bool_t AliAnalysisTaskCombinHF::IsKaon(AliAODTrack* track){
837 // kaon selection cuts
f2221737 838
811e9cbd 839 if(!fPidHF) return kTRUE;
f2221737 840 Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
841 Double_t mom=track->P();
a25d1750 842 if(SelectAODTrack(track,fTrackCutsKaon)) {
843 if(isKaon>=1) return kTRUE;
f2221737 844 if(isKaon<=-1) return kFALSE;
845 switch(fPIDselCaseZero){// isKaon=0
846 case 0:
847 {
848 return kTRUE;// always accept
849 }
850 break;
851 case 1:
852 {
853 if(isKaon>=0 && track->P()>fmaxPforIDKaon)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDKaon
854 }
855 break;
856 case 2:
857 {
858 if(track->P()>fmaxPforIDKaon)return kTRUE;
859 AliPIDResponse *pidResp=fPidHF->GetPidResponse();// more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
860 Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kKaon);
861 if(nsigma>-2.&& nsigma<3. && mom<0.6)isKaon=1;
862 else if(nsigma>-1.&& nsigma<3.&& mom<0.8)isKaon=1;
863 if(isKaon==1)return kTRUE;
864 }
865 break;
866 default:
867 {
868 AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
869 return kFALSE;// actually case 0 could be set as the default and return kTRUE
870 }
871 }
a25d1750 872 }
f2221737 873
811e9cbd 874 return kFALSE;
875}
f2221737 876//_______________________________________________________________________
811e9cbd 877Bool_t AliAnalysisTaskCombinHF::IsPion(AliAODTrack* track){
878 // pion selection cuts
f2221737 879
811e9cbd 880 if(!fPidHF) return kTRUE;
881 Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
f2221737 882 Double_t mom=track->P();
a25d1750 883 if(SelectAODTrack(track,fTrackCutsPion)) {
884 if(isPion>=1) return kTRUE;
f2221737 885 if(isPion<=-1) return kFALSE;
886 switch(fPIDselCaseZero){// isPion=0
887 case 0:
888 {
889 return kTRUE;// always accept
890 }
891 break;
892 case 1:
893 {
894 if(track->P()>fmaxPforIDPion)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDPion
895 }
896 break;
897 case 2:
898 {
899 // more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
900 if(track->P()>fmaxPforIDPion)return kTRUE;
901 AliPIDResponse *pidResp=fPidHF->GetPidResponse();
902 Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kPion);
903 if(nsigma<2.&& nsigma>-3. && mom<0.6)isPion=1;
904 else if(nsigma<1. && nsigma> -3. && mom<0.8)isPion=1;
905 if(isPion==1)return kTRUE;
906 }
907 break;
908 default:
909 {
910 AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
911 return kFALSE;// actually case 0 could be set as the default and return kTRUE
912 }
913 }
a25d1750 914 }
f2221737 915
811e9cbd 916 return kFALSE;
917}
f2221737 918
811e9cbd 919//________________________________________________________________________
920Bool_t AliAnalysisTaskCombinHF::SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts){
921 // AOD track selection
f2221737 922
811e9cbd 923 if(!cuts) return kTRUE;
f2221737 924
811e9cbd 925 AliESDtrack esdTrack(track);
926 // set the TPC cluster info
927 esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
928 esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
929 esdTrack.SetTPCPointsF(track->GetTPCNclsF());
f2221737 930 if(!cuts->IsSelected(&esdTrack)) return kFALSE;
931
932 return kTRUE;
811e9cbd 933}
934
ebc460de 935//_________________________________________________________________
936Bool_t AliAnalysisTaskCombinHF::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
937 // check if the decay products are in the good eta and pt range
938 for (Int_t iProng = 0; iProng<nProng; iProng++){
939 AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
940 if(!mcPartDaughter) return kFALSE;
941 Double_t eta = mcPartDaughter->Eta();
942 Double_t pt = mcPartDaughter->Pt();
943 if (TMath::Abs(eta) > fEtaAccCut || pt < fPtAccCut) return kFALSE;
944 }
945 return kTRUE;
946}
93f0ffc1 947//_________________________________________________________________
948Bool_t AliAnalysisTaskCombinHF::CanBeMixed(Double_t zv1, Double_t zv2, Double_t mult1, Double_t mult2){
949 if(fUsePoolsZ && fzVertPoolLims){
950 Int_t theBin1=TMath::BinarySearch(fNzVertPools+1,fzVertPoolLims,zv1);
951 Int_t theBin2=TMath::BinarySearch(fNzVertPools+1,fzVertPoolLims,zv2);
952 if(theBin1!=theBin2) return kFALSE;
953 if(theBin1<0 || theBin2<0) return kFALSE;
954 if(theBin1>=fNzVertPools || theBin2>=fNzVertPools) return kFALSE;
955 }else{
956 if(TMath::Abs(zv2-zv1)>fMaxzVertDistForMix) return kFALSE;
957 }
958
959 if(fUsePoolsM && fMultPoolLims){
960 Int_t theBin1=TMath::BinarySearch(fNMultPools+1,fMultPoolLims,mult1);
961 Int_t theBin2=TMath::BinarySearch(fNMultPools+1,fMultPoolLims,mult2);
962 if(theBin1!=theBin2) return kFALSE;
963 if(theBin1<0 || theBin2<0) return kFALSE;
964 if(theBin1>=fNMultPools || theBin2>=fNMultPools) return kFALSE;
965 }else{
966 if(TMath::Abs(mult2-mult1)>fMaxMultDiffForMix) return kFALSE;
967 }
968 return kTRUE;
969}
970
971//_________________________________________________________________
972void AliAnalysisTaskCombinHF::FinishTaskOutput()
973{
974 // perform mixed event analysis
975 if(!fDoEventMixing) return;
976
977 Int_t nEvents=fEventBuffer->GetEntries();
978 printf("Start Event Mixing of %d events\n",nEvents);
979
980 TObjArray* karray=0x0;
981 TObjArray* parray=0x0;
982 Double_t zVertex,mult;
983 fEventBuffer->SetBranchAddress("karray", &karray);
984 fEventBuffer->SetBranchAddress("parray", &parray);
985 fEventBuffer->SetBranchAddress("zVertex", &zVertex);
986 fEventBuffer->SetBranchAddress("multiplicity", &mult);
987
988 // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
989 Double_t d02[2]={0.,0.};
990 Double_t d03[3]={0.,0.,0.};
991 AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
992 AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
993 UInt_t pdg0[2]={321,211};
994 UInt_t pdgp[3]={321,211,211};
995 Double_t px[3],py[3],pz[3];
996
997 for(Int_t iEv1=0; iEv1<nEvents; iEv1++){
998 fEventBuffer->GetEvent(iEv1);
999 Double_t zVertex1=zVertex;
1000 Double_t mult1=mult;
1001 TObjArray* karray1=new TObjArray(*karray);
1002 for(Int_t iEv2=0; iEv2<fMaxNumberOfEventsForMixing; iEv2++){
1003 Int_t iToMix=iEv1+iEv2+1;
1004 if(iEv1>=(nEvents-fMaxNumberOfEventsForMixing)) iToMix=iEv1-iEv2-1;
1005 if(iToMix<0) continue;
1006 if(iToMix==iEv1) continue;
1007 fEventBuffer->GetEvent(iToMix);
1008 Double_t zVertex2=zVertex;
1009 Double_t mult2=mult;
1010 if(CanBeMixed(zVertex1,zVertex2,mult1,mult2)){
1011 TObjArray* parray2=new TObjArray(*parray);
1012 Int_t nKaons=karray1->GetEntries();
1013 Int_t nPions=parray2->GetEntries();
1014 Bool_t doThird=kTRUE;
1015 Int_t iToMix3=iEv1+2*fMaxNumberOfEventsForMixing+1;
1016 if(iToMix3>=nEvents) iToMix3=iEv1-2*fMaxNumberOfEventsForMixing-1;
1017 if(iToMix3<0 || iToMix3==iEv1 || iToMix3==iToMix) doThird=kFALSE;
1018 TObjArray* parray3=0x0;
1019 Double_t zVertex3=-999.;
1020 Double_t mult3=999999;
1021 if(fMeson!=kDzero && doThird){
1022 fEventBuffer->GetEvent(iToMix3);
1023 zVertex3=zVertex;
1024 mult3=mult;
1025 if(CanBeMixed(zVertex1,zVertex3,mult1,mult3)){
1026 parray3=new TObjArray(*parray);
1027 }else{
1028 doThird=kFALSE;
1029 }
1030 }
ebc460de 1031
93f0ffc1 1032 for(Int_t iTr1=0; iTr1<nKaons; iTr1++){
1033 TLorentzVector* trK=(TLorentzVector*)karray1->At(iTr1);
1034 Double_t chargeK=trK->T();
1035 px[0] = trK->Px();
1036 py[0] = trK->Py();
1037 pz[0] = trK->Pz();
1038 for(Int_t iTr2=0; iTr2<nPions; iTr2++){
1039 TLorentzVector* trPi1=(TLorentzVector*)parray2->At(iTr2);
1040 Double_t chargePi1=trPi1->T();
1041 px[1] = trPi1->Px();
1042 py[1] = trPi1->Py();
1043 pz[1] = trPi1->Pz();
1044 if(chargePi1*chargeK<0){
1045 if(fMeson==kDzero){
1046 FillMEHistos(421,2,tmpRD2,px,py,pz,pdg0);
1047 }else{
1048 if(doThird && parray3){
1049 Int_t nPions3=parray3->GetEntries();
1050 for(Int_t iTr3=iTr2+1; iTr3<nPions3; iTr3++){
1051 TLorentzVector* trPi2=(TLorentzVector*)parray3->At(iTr3);
1052 Double_t chargePi2=trPi2->T();
1053 px[2] = trPi2->Px();
1054 py[2] = trPi2->Py();
1055 pz[2] = trPi2->Pz();
1056 if(chargePi2*chargeK<0){
1057 FillMEHistos(411,3,tmpRD3,px,py,pz,pdgp);
1058 }
1059 }
1060 }
1061 }
1062 }
1063 }
1064 }
1065 delete parray3;
1066 delete parray2;
1067 }
1068 }
1069 delete karray1;
1070 }
1071 delete tmpRD2;
1072 delete tmpRD3;
93f0ffc1 1073}
811e9cbd 1074//_________________________________________________________________
1075void AliAnalysisTaskCombinHF::Terminate(Option_t */*option*/)
1076{
1077 // Terminate analysis
1078 //
1079 if(fDebug > 1) printf("AliAnalysisTaskCombinHF: Terminate() \n");
1080 fOutput = dynamic_cast<TList*> (GetOutputData(1));
f2221737 1081 if (!fOutput) {
811e9cbd 1082 printf("ERROR: fOutput not available\n");
1083 return;
1084 }
1085 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
1086 if(fHistNEvents){
1087 printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1088 }else{
1089 printf("ERROR: fHistNEvents not available\n");
1090 return;
1091 }
1092 return;
1093}
1094