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