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