]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliAnalysisTaskSEDs.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Possibility to use PID in the AOD filtering for Lc and Ds
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSEDs.cxx
CommitLineData
1f796bc2 1/**************************************************************************
2 * Copyright(c) 2007-2009, 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// //
20// Analysis task to produce Ds candidates mass spectra //
21// Origin: F.Prino, Torino, prino@to.infn.it //
22// //
23///////////////////////////////////////////////////////////////////
24
25#include <TClonesArray.h>
26#include <TNtuple.h>
27#include <TList.h>
28#include <TString.h>
29#include <TH1F.h>
30#include <TMath.h>
31#include <TDatabasePDG.h>
32
33#include "AliAnalysisManager.h"
34#include "AliAODHandler.h"
35#include "AliAODEvent.h"
36#include "AliAODVertex.h"
37#include "AliAODTrack.h"
38#include "AliAODMCHeader.h"
39#include "AliAODMCParticle.h"
40#include "AliAODRecoDecayHF3Prong.h"
41#include "AliAnalysisVertexingHF.h"
42#include "AliRDHFCutsDstoKKpi.h"
43#include "AliAnalysisTaskSE.h"
44#include "AliNormalizationCounter.h"
45#include "AliAnalysisTaskSEDs.h"
46
47ClassImp(AliAnalysisTaskSEDs)
48
49
50//________________________________________________________________________
51AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():
52 AliAnalysisTaskSE(),
53 fOutput(0),
54 fHistNEvents(0),
55 fPtVsMass(0),
811a34ce 56 fPtVsMassPhi(0),
57 fPtVsMassK0st(0),
1f796bc2 58 fYVsPt(0),
59 fYVsPtSig(0),
60 fNtupleDs(0),
61 fFillNtuple(0),
62 fReadMC(kFALSE),
63 fDoCutVarHistos(kTRUE),
64 fUseSelectionBit(kFALSE),
65 fNPtBins(0),
66 fListCuts(0),
67 fMassRange(0.8),
68 fMassBinSize(0.002),
69 fCounter(0),
70 fAnalysisCuts(0)
71{
72 // Default constructor
73
74 for(Int_t i=0;i<4;i++) {
75 if(fChanHist[i]) fChanHist[i]=0;
76 }
77
78 for(Int_t i=0;i<4*kMaxPtBins;i++){
79
80 if(fPtCandHist[i]) fPtCandHist[i]=0;
81 if(fMassHist[i]) fMassHist[i]=0;
82 if(fCosPHist[i]) fCosPHist[i]=0;
83 if(fDLenHist[i]) fDLenHist[i]=0;
84 if(fSumd02Hist[i]) fSumd02Hist[i]=0;
85 if(fSigVertHist[i]) fSigVertHist[i]=0;
86 if(fPtMaxHist[i]) fPtMaxHist[i]=0;
87 if(fDCAHist[i]) fDCAHist[i]=0;
88 if(fPtProng0Hist[i]) fPtProng0Hist[i]=0;
89 if(fPtProng1Hist[i]) fPtProng1Hist[i]=0;
90 if(fPtProng2Hist[i]) fPtProng2Hist[i]=0;
91 if(fDalitz[i]) fDalitz[i]=0;
92 if(fDalitzPhi[i]) fDalitzPhi[i]=0;
93 if(fDalitzK0st[i]) fDalitzK0st[i]=0;
94
95 }
96 for(Int_t i=0;i<3*kMaxPtBins;i++){
97 if(fMassHistPhi[i]) fMassHistPhi[i]=0;
98 if(fMassHistK0st[i]) fMassHistK0st[i]=0;
99
100 }
101
102 for(Int_t i=0;i<kMaxPtBins+1;i++){
103 fPtLimits[i]=0;
104 }
105
106}
107
108//________________________________________________________________________
109AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name,AliRDHFCutsDstoKKpi* analysiscuts,Int_t fillNtuple):
110 AliAnalysisTaskSE(name),
111 fOutput(0),
112 fHistNEvents(0),
113 fPtVsMass(0),
811a34ce 114 fPtVsMassPhi(0),
115 fPtVsMassK0st(0),
1f796bc2 116 fYVsPt(0),
117 fYVsPtSig(0),
118 fNtupleDs(0),
119 fFillNtuple(fillNtuple),
120 fReadMC(kFALSE),
121 fDoCutVarHistos(kTRUE),
122 fUseSelectionBit(kFALSE),
123 fNPtBins(0),
124 fListCuts(0),
125 fMassRange(0.8),
126 fMassBinSize(0.002),
127 fCounter(0),
128 fAnalysisCuts(analysiscuts)
129{
130 // Default constructor
131 // Output slot #1 writes into a TList container
132
133 for(Int_t i=0;i<4;i++) {
134 if(fChanHist[i]) fChanHist[i]=0;
135 }
136
137 for(Int_t i=0;i<4*kMaxPtBins;i++){
138
139 if(fPtCandHist[i]) fPtCandHist[i]=0;
140 if(fMassHist[i]) fMassHist[i]=0;
141 if(fCosPHist[i]) fCosPHist[i]=0;
142 if(fDLenHist[i]) fDLenHist[i]=0;
143 if(fSumd02Hist[i]) fSumd02Hist[i]=0;
144 if(fSigVertHist[i]) fSigVertHist[i]=0;
145 if(fPtMaxHist[i]) fPtMaxHist[i]=0;
146 if(fDCAHist[i]) fDCAHist[i]=0;
147 if(fPtProng0Hist[i]) fPtProng0Hist[i]=0;
148 if(fPtProng1Hist[i]) fPtProng1Hist[i]=0;
149 if(fPtProng2Hist[i]) fPtProng2Hist[i]=0;
150 if(fDalitz[i]) fDalitz[i]=0;
151 if(fDalitzPhi[i]) fDalitzPhi[i]=0;
152 if(fDalitzK0st[i]) fDalitzK0st[i]=0;
153
154 }
155 for(Int_t i=0;i<3*kMaxPtBins;i++){
156 if(fMassHistPhi[i]) fMassHistPhi[i]=0;
157 if(fMassHistK0st[i]) fMassHistK0st[i]=0;
158
159 }
160
161 for(Int_t i=0;i<kMaxPtBins+1;i++){
162 fPtLimits[i]=0;
163 }
164
165 Int_t nptbins=fAnalysisCuts->GetNPtBins();
166 Float_t *ptlim=fAnalysisCuts->GetPtBinLimits();
167 SetPtBins(nptbins,ptlim);
168
169 DefineOutput(1,TList::Class()); //My private output
170
171 DefineOutput(2,TList::Class());
172
173 DefineOutput(3,AliNormalizationCounter::Class());
174
175 if(fFillNtuple>0){
176 // Output slot #4 writes into a TNtuple container
177 DefineOutput(4,TNtuple::Class()); //My private output
178 }
179
180}
181
182//________________________________________________________________________
183void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Float_t* lim){
184 // define pt bins for analysis
185 if(n>kMaxPtBins){
186 printf("Max. number of Pt bins = %d\n",kMaxPtBins);
187 fNPtBins=kMaxPtBins;
188 fPtLimits[0]=0.;
189 fPtLimits[1]=1.;
190 fPtLimits[2]=3.;
191 fPtLimits[3]=5.;
192 fPtLimits[4]=10.;
193 for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
194 }else{
195 fNPtBins=n;
196 for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
197 for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
198 }
199 if(fDebug > 1){
200 printf("Number of Pt bins = %d\n",fNPtBins);
201 for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);
202 }
203}
204//________________________________________________________________________
205AliAnalysisTaskSEDs::~AliAnalysisTaskSEDs()
206{
207 // Destructor
208 delete fOutput;
209 delete fHistNEvents;
210 delete fListCuts;
211
212 for(Int_t i=0;i<4*fNPtBins;i++){
213
214 if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}
215 if(fMassHistPhi[i]){ delete fMassHistPhi[i]; fMassHistPhi[i]=0;}
216 if(fMassHistK0st[i]){ delete fMassHistK0st[i]; fMassHistK0st[i]=0;}
217 if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}
218 if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}
219 if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}
220 if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}
221 if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}
222 if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}
223 if(fPtProng0Hist[i]){ delete fPtProng0Hist[i]; fPtProng0Hist[i]=0;}
224 if(fPtProng1Hist[i]){ delete fPtProng1Hist[i]; fPtProng1Hist[i]=0;}
225 if(fPtProng2Hist[i]){ delete fPtProng2Hist[i]; fPtProng2Hist[i]=0;}
226 if(fDalitz[i]){ delete fDalitz[i]; fDalitz[i]=0;}
227 if(fDalitzPhi[i]){ delete fDalitzPhi[i]; fDalitzPhi[i]=0;}
228 if(fDalitzK0st[i]){ delete fDalitzK0st[i]; fDalitzK0st[i]=0;}
229
230 }
231
232 delete fPtVsMass;
811a34ce 233 delete fPtVsMassPhi;
234 delete fPtVsMassK0st;
1f796bc2 235 delete fYVsPt;
236 delete fYVsPtSig;
237 delete fNtupleDs;
238 delete fCounter;
239 delete fAnalysisCuts;
240
241}
242
243//________________________________________________________________________
244void AliAnalysisTaskSEDs::Init()
245{
246 // Initialization
247
248 if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
249
250 fListCuts=new TList();
251 fListCuts->SetOwner();
252 fListCuts->SetName("CutObjects");
253
254 AliRDHFCutsDstoKKpi *analysis = new AliRDHFCutsDstoKKpi(*fAnalysisCuts);
255 analysis->SetName("AnalysisCuts");
256
257 fListCuts->Add(analysis);
258 PostData(2,fListCuts);
259 return;
260}
261
262//________________________________________________________________________
263void AliAnalysisTaskSEDs::UserCreateOutputObjects()
264{
265 // Create the output container
266 //
267 if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
268
269 // Several histograms are more conveniently managed in a TList
270 fOutput = new TList();
271 fOutput->SetOwner();
272 fOutput->SetName("OutputHistos");
273
274 Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
275
276 Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);
277 if(nInvMassBins%2==1) nInvMassBins++;
278 // Double_t minMass=massDs-1.0*nInvMassBins*fMassBinSize;
279 Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
280 // Double_t maxMass=massDs+1.0*nInvMassBins*fMassBinSize;
281 Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
282
283 TString hisname;
284 TString htype;
285 Int_t index;
286 for(Int_t iType=0; iType<4; iType++){
287 for(Int_t i=0;i<fNPtBins;i++){
288 if(iType==0){
289 htype="All";
290 index=GetHistoIndex(i);
291 }else if(iType==1){
292 htype="Sig";
293 index=GetSignalHistoIndex(i);
294 }else if(iType==2){
295 htype="Bkg";
296 index=GetBackgroundHistoIndex(i);
297 }else{
298 htype="ReflSig";
299 index=GetReflSignalHistoIndex(i);
300 }
301 hisname.Form("hMass%sPt%d",htype.Data(),i);
302 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
303 fMassHist[index]->Sumw2();
304 hisname.Form("hMass%sPt%dphi",htype.Data(),i);
305 fMassHistPhi[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
306 fMassHistPhi[index]->Sumw2();
307 hisname.Form("hMass%sPt%dk0st",htype.Data(),i);
308 fMassHistK0st[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
309 fMassHistK0st[index]->Sumw2();
310 hisname.Form("hCosP%sPt%d",htype.Data(),i);
311 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
312 fCosPHist[index]->Sumw2();
313 hisname.Form("hDLen%sPt%d",htype.Data(),i);
314 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
315 fDLenHist[index]->Sumw2();
316 hisname.Form("hSumd02%sPt%d",htype.Data(),i);
317 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
318 fSumd02Hist[index]->Sumw2();
319 hisname.Form("hSigVert%sPt%d",htype.Data(),i);
320 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
321 fSigVertHist[index]->Sumw2();
322 hisname.Form("hPtMax%sPt%d",htype.Data(),i);
323 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
324 fPtMaxHist[index]->Sumw2();
325 hisname.Form("hPtCand%sPt%d",htype.Data(),i);
326 fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
327 fPtCandHist[index]->Sumw2();
328 hisname.Form("hDCA%sPt%d",htype.Data(),i);
329 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
330 fDCAHist[index]->Sumw2();
331 hisname.Form("hPtProng0%sPt%d",htype.Data(),i);
332 fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
333 fPtProng0Hist[index]->Sumw2();
334 hisname.Form("hPtProng1%sPt%d",htype.Data(),i);
335 fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
336 fPtProng1Hist[index]->Sumw2();
337 hisname.Form("hPtProng2%sPt%d",htype.Data(),i);
338 fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
339 fPtProng2Hist[index]->Sumw2();
340 hisname.Form("hDalitz%sPt%d",htype.Data(),i);
341 fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
342 fDalitz[index]->Sumw2();
343 hisname.Form("hDalitz%sPt%dphi",htype.Data(),i);
344 fDalitzPhi[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
345 fDalitzPhi[index]->Sumw2();
346 hisname.Form("hDalitz%sPt%dk0st",htype.Data(),i);
347 fDalitzK0st[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
348 fDalitzK0st[index]->Sumw2();
349 }
350 }
351
352 for(Int_t i=0; i<4*fNPtBins; i++){
353 fOutput->Add(fMassHist[i]);
354 fOutput->Add(fMassHistPhi[i]);
355 fOutput->Add(fMassHistK0st[i]);
356 fOutput->Add(fCosPHist[i]);
357 fOutput->Add(fDLenHist[i]);
358 fOutput->Add(fSumd02Hist[i]);
359 fOutput->Add(fSigVertHist[i]);
360 fOutput->Add(fPtMaxHist[i]);
361 fOutput->Add(fPtCandHist[i]);
362 fOutput->Add(fDCAHist[i]);
363 fOutput->Add(fPtProng0Hist[i]);
364 fOutput->Add(fPtProng1Hist[i]);
365 fOutput->Add(fPtProng2Hist[i]);
366 fOutput->Add(fDalitz[i]);
367 fOutput->Add(fDalitzPhi[i]);
368 fOutput->Add(fDalitzK0st[i]);
369 }
370
371 fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",64,-0.5,63.5);
372 fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",64,-0.5,63.5);
373 fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",64,-0.5,63.5);
374 fChanHist[3] = new TH1F("hChanReflSig", "KKpi and piKK candidates",64,-0.5,63.5);
375 for(Int_t i=0;i<4;i++){
376 fChanHist[i]->Sumw2();
377 fChanHist[i]->SetMinimum(0);
378 fOutput->Add(fChanHist[i]);
379 }
380
381 fHistNEvents = new TH1F("hNEvents", "number of events ",11,-0.5,10.5);
382 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
383 fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
384 fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
385 fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to not reco vertex");
386 fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected for contr vertex");
387 fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for vertex out of accept");
388 fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for pileup events");
389 fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
390 fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of 3 prong candidates");
391 fHistNEvents->GetXaxis()->SetBinLabel(10,"no. of Ds after filtering cuts");
392 fHistNEvents->GetXaxis()->SetBinLabel(11,"no. of Ds after selection cuts");
393
394 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
395
396 fHistNEvents->Sumw2();
397 fHistNEvents->SetMinimum(0);
398 fOutput->Add(fHistNEvents);
399
400 fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nInvMassBins,minMass,maxMass,40,0.,20.);
811a34ce 401 fPtVsMassPhi=new TH2F("hPtVsMassPhi","PtVsMass (phi selection)",nInvMassBins,minMass,maxMass,200,0.,20.);
402 fPtVsMassK0st=new TH2F("hPtVsMassK0st","PtVsMass (K0* selection)",nInvMassBins,minMass,maxMass,200,0.,20.);
1f796bc2 403 fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
404 fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
405
406 fOutput->Add(fPtVsMass);
811a34ce 407 fOutput->Add(fPtVsMassPhi);
408 fOutput->Add(fPtVsMassK0st);
1f796bc2 409 fOutput->Add(fYVsPt);
410 fOutput->Add(fYVsPtSig);
411
412 //Counter for Normalization
413 fCounter = new AliNormalizationCounter("NormalizationCounter");
414 fCounter->Init();
415
416 PostData(1,fOutput);
417 PostData(3,fCounter);
418
419 if(fFillNtuple>0){
420 OpenFile(4); // 4 is the slot number of the ntuple
421
8d4de239 422 fNtupleDs = new TNtuple("fNtupleDs","Ds","labDs:retcode:pdgcode0:Pt0:Pt1:Pt2:PtRec:P0:P1:P2:PidTrackBit0:PidTrackBit1:PidTrackBit2:PointingAngle:PointingAngleXY:DecLeng:DecLengXY:NorDecLeng:NorDecLengXY:InvMassKKpi:InvMasspiKK:sigvert:d00:d01:d02:dca:d0square:InvMassPhiKKpi:InvMassPhipiKK:InvMassK0starKKpi:InvMassK0starpiKK:cosinePiDsFrameKKpi:cosinePiDsFramepiKK:cosineKPhiFrameKKpi:cosineKPhiFramepiKK:centrality:runNumber");
1f796bc2 423
424 }
425
426
427 return;
428}
429
430//________________________________________________________________________
431void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
432{
433 // Ds selection for current event, fill mass histos and selecetion variable histo
434 // separate signal and backgound if fReadMC is activated
435
436 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
437
438 TClonesArray *array3Prong = 0;
439 if(!aod && AODEvent() && IsStandardAOD()) {
440 // In case there is an AOD handler writing a standard AOD, use the AOD
441 // event in memory rather than the input (ESD) event.
442 aod = dynamic_cast<AliAODEvent*> (AODEvent());
443 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
444 // have to taken from the AOD event hold by the AliAODExtension
445 AliAODHandler* aodHandler = (AliAODHandler*)
446 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
447 if(aodHandler->GetExtensions()) {
448 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
449 AliAODEvent *aodFromExt = ext->GetAOD();
450 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
451 }
452 } else if(aod) {
453 array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
454 }
455
456 if(!aod || !array3Prong) {
457 printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
458 return;
459 }
460
461
462 // fix for temporary bug in ESDfilter
463 // the AODs with null vertex pointer didn't pass the PhysSel
464 if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
465
466
467 fHistNEvents->Fill(0); // count event
468 // Post the data already here
469 PostData(1,fOutput);
470
471 fCounter->StoreEvent(aod,fAnalysisCuts,fReadMC);
472
473
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->IsEventRejectedDueToPileupSPD())fHistNEvents->Fill(6);
480 if(fAnalysisCuts->IsEventRejectedDueToCentrality())fHistNEvents->Fill(7);
481
482 Float_t centrality=fAnalysisCuts->GetCentrality(aod);
8d4de239 483 Int_t runNumber=aod->GetRunNumber();
1f796bc2 484
485 if(!isEvSel)return;
486
487
488
489
490 fHistNEvents->Fill(1);
491
492 TClonesArray *arrayMC=0;
493 AliAODMCHeader *mcHeader=0;
494
495 // AOD primary vertex
496 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
497 // vtx1->Print();
498
499 // load MC particles
500 if(fReadMC){
501
502 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
503 if(!arrayMC) {
504 printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
505 return;
506 }
507
508 // load MC header
509 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
510 if(!mcHeader) {
511 printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
512 return;
513 }
514 }
515
516 Int_t n3Prong = array3Prong->GetEntriesFast();
517 if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
518
519
520 Int_t pdgDstoKKpi[3]={321,321,211};
521 Int_t nSelected=0;
522 Int_t nFiltered=0;
523 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
524
525 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
526 fHistNEvents->Fill(8);
527
528 if(fUseSelectionBit && !(d->HasSelectionBit(AliRDHFCuts::kDsCuts))){
529 continue;
530 }
531 nFiltered++;
532 fHistNEvents->Fill(9);
533
534 Bool_t unsetvtx=kFALSE;
535 if(!d->GetOwnPrimaryVtx()){
536 d->SetOwnPrimaryVtx(vtx1);
537 unsetvtx=kTRUE;
538 }
539
540 Bool_t recVtx=kFALSE;
541 AliAODVertex *origownvtx=0x0;
542
543
544 Double_t ptCand = d->Pt();
545 Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
48abf62f 546 Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
1f796bc2 547 Double_t rapid=d->YDs();
548 fYVsPt->Fill(ptCand,rapid);
549
550 if(retCodeAnalysisCuts<=0) continue;
551 Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);
552 if(!isFidAcc) continue;
553
554 if(fAnalysisCuts->GetIsPrimaryWithoutDaughters()){
555 if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
556 if(fAnalysisCuts->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
557 else fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
558 }
559
560
561 fHistNEvents->Fill(10);
562 nSelected++;
563
564 Int_t index=GetHistoIndex(iPtBin);
565 fPtCandHist[index]->Fill(ptCand);
566
567 Int_t isKKpi=retCodeAnalysisCuts&1;
568 Int_t ispiKK=retCodeAnalysisCuts&2;
569 Int_t isPhiKKpi=retCodeAnalysisCuts&4;
570 Int_t isPhipiKK=retCodeAnalysisCuts&8;
571 Int_t isK0starKKpi=retCodeAnalysisCuts&16;
572 Int_t isK0starpiKK=retCodeAnalysisCuts&32;
573
574 fChanHist[0]->Fill(retCodeAnalysisCuts);
575
576 Int_t indexMCKKpi=-1;
577 Int_t indexMCpiKK=-1;
578 Int_t labDs=-1;
579 Int_t pdgCode0=-999;
580
581 if(fReadMC){
582 labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
583 if(labDs>=0){
584 Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
585 AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
586 pdgCode0=TMath::Abs(p->GetPdgCode());
587
588 if(isKKpi){
589 if(pdgCode0==321) {
590 indexMCKKpi=GetSignalHistoIndex(iPtBin);
591 fYVsPtSig->Fill(ptCand,rapid);
592 fChanHist[1]->Fill(retCodeAnalysisCuts);
593 }else{
594 indexMCKKpi=GetReflSignalHistoIndex(iPtBin);
595 fChanHist[3]->Fill(retCodeAnalysisCuts);
596 }
597 }
598 if(ispiKK){
599 if(pdgCode0==211) {
600 indexMCpiKK=GetSignalHistoIndex(iPtBin);
601 fYVsPtSig->Fill(ptCand,rapid);
602 fChanHist[1]->Fill(retCodeAnalysisCuts);
603 }else{
604 indexMCpiKK=GetReflSignalHistoIndex(iPtBin);
605 fChanHist[3]->Fill(retCodeAnalysisCuts);
606 }
607 }
608 }else{
609 indexMCpiKK=GetBackgroundHistoIndex(iPtBin);
610 indexMCKKpi=GetBackgroundHistoIndex(iPtBin);
611 fChanHist[2]->Fill(retCodeAnalysisCuts);
612 }
613 }
614
615 if(isKKpi){
616 Double_t invMass=d->InvMassDsKKpi();
617 fMassHist[index]->Fill(invMass);
618 fPtVsMass->Fill(invMass,ptCand);
811a34ce 619 if(isPhiKKpi){
620 fMassHistPhi[index]->Fill(invMass);
621 fPtVsMassPhi->Fill(invMass,ptCand);
622 }
623 if(isK0starKKpi){
624 fMassHistK0st[index]->Fill(invMass);
625 fPtVsMassK0st->Fill(invMass,ptCand);
626 }
1f796bc2 627 if(fReadMC && indexMCKKpi!=-1){
628 fMassHist[indexMCKKpi]->Fill(invMass);
629 if(isPhiKKpi) fMassHistPhi[indexMCKKpi]->Fill(invMass);
630 if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass);
631 }
632 }
633 if(ispiKK){
634 Double_t invMass=d->InvMassDspiKK();
635 fMassHist[index]->Fill(invMass);
636 fPtVsMass->Fill(invMass,ptCand);
811a34ce 637 if(isPhipiKK){
638 fMassHistPhi[index]->Fill(invMass);
639 fPtVsMassPhi->Fill(invMass,ptCand);
640 }
641 if(isK0starpiKK){
642 fMassHistK0st[index]->Fill(invMass);
643 fPtVsMassK0st->Fill(invMass,ptCand);
644 }
1f796bc2 645 if(fReadMC && indexMCpiKK!=-1){
646 fMassHist[indexMCpiKK]->Fill(invMass);
647 if(isPhipiKK) fMassHistPhi[indexMCpiKK]->Fill(invMass);
648 if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass);
649 }
650 }
651
652 if(fDoCutVarHistos){
653 Double_t dlen=d->DecayLength();
654 Double_t cosp=d->CosPointingAngle();
655 Double_t pt0=d->PtProng(0);
656 Double_t pt1=d->PtProng(1);
657 Double_t pt2=d->PtProng(2);
658 Double_t sigvert=d->GetSigmaVert();
659 Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
660 Double_t dca=d->GetDCA();
661 Double_t ptmax=0;
662 for(Int_t i=0;i<3;i++){
663 if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
664 }
665 fCosPHist[index]->Fill(cosp);
666 fDLenHist[index]->Fill(dlen);
667 fSigVertHist[index]->Fill(sigvert);
668 fSumd02Hist[index]->Fill(sumD02);
669 fPtMaxHist[index]->Fill(ptmax);
670 fDCAHist[index]->Fill(dca);
671 fPtProng0Hist[index]->Fill(pt0);
672 fPtProng1Hist[index]->Fill(pt1);
673 fPtProng2Hist[index]->Fill(pt2);
674 if(isKKpi){
675 Double_t massKK=d->InvMass2Prongs(0,1,321,321);
676 Double_t massKp=d->InvMass2Prongs(1,2,321,211);
677 fDalitz[index]->Fill(massKK,massKp);
678 if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);
679 if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);
680 if(fReadMC && indexMCKKpi!=-1){
681 fDalitz[indexMCKKpi]->Fill(massKK,massKp);
682 if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);
683 if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);
684 fCosPHist[indexMCKKpi]->Fill(cosp);
685 fDLenHist[indexMCKKpi]->Fill(dlen);
686 fSigVertHist[indexMCKKpi]->Fill(sigvert);
687 fSumd02Hist[indexMCKKpi]->Fill(sumD02);
688 fPtMaxHist[indexMCKKpi]->Fill(ptmax);
689 fPtCandHist[indexMCKKpi]->Fill(ptCand);
690 fDCAHist[indexMCKKpi]->Fill(dca);
691 fPtProng0Hist[indexMCKKpi]->Fill(pt0);
692 fPtProng1Hist[indexMCKKpi]->Fill(pt1);
693 fPtProng2Hist[indexMCKKpi]->Fill(pt2);
694 }
695 }
696 if(ispiKK){
697 Double_t massKK=d->InvMass2Prongs(1,2,321,321);
698 Double_t massKp=d->InvMass2Prongs(0,1,211,321);
699 fDalitz[index]->Fill(massKK,massKp);
700 if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);
701 if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);
702
703
704 if(fReadMC && indexMCpiKK!=-1){
705 fDalitz[indexMCpiKK]->Fill(massKK,massKp);
706 if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);
707 if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);
708 fCosPHist[indexMCpiKK]->Fill(cosp);
709 fDLenHist[indexMCpiKK]->Fill(dlen);
710 fSigVertHist[indexMCpiKK]->Fill(sigvert);
711 fSumd02Hist[indexMCpiKK]->Fill(sumD02);
712 fPtMaxHist[indexMCpiKK]->Fill(ptmax);
713 fPtCandHist[indexMCpiKK]->Fill(ptCand);
714 fDCAHist[indexMCpiKK]->Fill(dca);
715 fPtProng0Hist[indexMCpiKK]->Fill(pt0);
716 fPtProng1Hist[indexMCpiKK]->Fill(pt1);
717 fPtProng2Hist[indexMCpiKK]->Fill(pt2);
718 }
719 }
720
721
722 }
723
8d4de239 724 Float_t tmp[37];
1f796bc2 725 if(fFillNtuple>0){
726
727 if ((fFillNtuple==1 && (isPhiKKpi || isPhipiKK)) || (fFillNtuple==2 && (isK0starKKpi || isK0starpiKK)) || (fFillNtuple==3 && (isKKpi || ispiKK))){
728
729 AliAODTrack *track0=(AliAODTrack*)d->GetDaughter(0);
730 AliAODTrack *track1=(AliAODTrack*)d->GetDaughter(1);
731 AliAODTrack *track2=(AliAODTrack*)d->GetDaughter(2);
732
733 UInt_t BitMapPIDTrack0=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track0);
734 UInt_t BitMapPIDTrack1=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track1);
735 UInt_t BitMapPIDTrack2=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track2);
736
737 tmp[0]=Float_t(labDs);
738 tmp[1]=Float_t(retCodeAnalysisCuts);
739 tmp[2]=Float_t(pdgCode0);
740 tmp[3]=d->PtProng(0);
741 tmp[4]=d->PtProng(1);
742 tmp[5]=d->PtProng(2);
743 tmp[6]=d->Pt();
744 tmp[7]=d->PProng(0);
745 tmp[8]=d->PProng(1);
746 tmp[9]=d->PProng(2);
747 tmp[10]=Int_t(BitMapPIDTrack0);
748 tmp[11]=Int_t(BitMapPIDTrack1);
749 tmp[12]=Int_t(BitMapPIDTrack2);
750 tmp[13]=d->CosPointingAngle();
751 tmp[14]=d->CosPointingAngleXY();
752 tmp[15]=d->DecayLength();
753 tmp[16]=d->DecayLengthXY();
754 tmp[17]=d->NormalizedDecayLength();
755 tmp[18]=d->NormalizedDecayLengthXY();
756 tmp[19]=d->InvMassDsKKpi();
757 tmp[20]=d->InvMassDspiKK();
758 tmp[21]=d->GetSigmaVert();
759 tmp[22]=d->Getd0Prong(0);
760 tmp[23]=d->Getd0Prong(1);
761 tmp[24]=d->Getd0Prong(2);
762 tmp[25]=d->GetDCA();
763 tmp[26]=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
764 tmp[27]=d->InvMass2Prongs(0,1,321,321);
765 tmp[28]=d->InvMass2Prongs(1,2,321,321);
766 tmp[29]=d->InvMass2Prongs(1,2,321,211);
767 tmp[30]=d->InvMass2Prongs(0,1,211,321);
768 tmp[31]=d->CosPiDsLabFrameKKpi();
769 tmp[32]=d->CosPiDsLabFramepiKK();
770 tmp[33]=d->CosPiKPhiRFrameKKpi();
771 tmp[34]=d->CosPiKPhiRFramepiKK();
772 tmp[35]=(Float_t)(centrality);
8d4de239 773 tmp[36]=(Float_t)(runNumber);
1f796bc2 774
1f796bc2 775 fNtupleDs->Fill(tmp);
776 PostData(4,fNtupleDs);
777 }
778 }
779
780 if(unsetvtx) d->UnsetOwnPrimaryVtx();
781 if(recVtx)fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
782 }
783
784 fCounter->StoreCandidates(aod,nFiltered,kTRUE);
785 fCounter->StoreCandidates(aod,nSelected,kFALSE);
786
787 PostData(1,fOutput);
788 PostData(3,fCounter);
789
790 return;
791}
792
793//_________________________________________________________________
794
795void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
796{
797 // Terminate analysis
798 //
799 if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
800 fOutput = dynamic_cast<TList*> (GetOutputData(1));
801 if (!fOutput) {
802 printf("ERROR: fOutput not available\n");
803 return;
804 }
805 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
806 if(fHistNEvents){
807 printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
808 }else{
809 printf("ERROR: fHistNEvents not available\n");
810 return;
811 }
812 return;
813}
814