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