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