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