]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliAnalysisTaskSEDs.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Faster calculation of invariant mass (Yifei)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDs.cxx
CommitLineData
25c94fc8 1/**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
c2a0a9d5 16/* $Id: $ */
25c94fc8 17
18///////////////////////////////////////////////////////////////////
19// //
20// Analysis task to produce Ds candidates mass spectra //
21// Origin: F.Prino, Torino, prino@to.infn.it //
22// //
23///////////////////////////////////////////////////////////////////
24
25#include <TClonesArray.h>
26#include <TNtuple.h>
27#include <TList.h>
28#include <TString.h>
29#include <TH1F.h>
30#include <TMath.h>
31#include <TDatabasePDG.h>
32
33#include "AliAnalysisManager.h"
34#include "AliAODHandler.h"
35#include "AliAODEvent.h"
36#include "AliAODVertex.h"
37#include "AliAODTrack.h"
38#include "AliAODMCHeader.h"
39#include "AliAODMCParticle.h"
40#include "AliAODRecoDecayHF3Prong.h"
41#include "AliAnalysisVertexingHF.h"
4c230f6d 42#include "AliRDHFCutsDstoKKpi.h"
25c94fc8 43#include "AliAnalysisTaskSE.h"
44#include "AliAnalysisTaskSEDs.h"
45
46ClassImp(AliAnalysisTaskSEDs)
47
48
49//________________________________________________________________________
50AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():
51 AliAnalysisTaskSE(),
52 fOutput(0),
53 fHistNEvents(0),
54 fReadMC(kFALSE),
55 fNPtBins(0),
4c230f6d 56 fListCuts(0),
25c94fc8 57 fMassRange(0.2),
a0bb2d06 58 fMassBinSize(0.002),
4c230f6d 59 fProdCuts(0),
60 fAnalysisCuts(0)
25c94fc8 61{
62 // Default constructor
63}
64
65//________________________________________________________________________
4c230f6d 66AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name, AliRDHFCutsDstoKKpi* productioncuts, AliRDHFCutsDstoKKpi* analysiscuts):
25c94fc8 67 AliAnalysisTaskSE(name),
68 fOutput(0),
69 fHistNEvents(0),
70 fReadMC(kFALSE),
71 fNPtBins(0),
4c230f6d 72 fListCuts(0),
25c94fc8 73 fMassRange(0.2),
a0bb2d06 74 fMassBinSize(0.002),
4c230f6d 75 fProdCuts(productioncuts),
76 fAnalysisCuts(analysiscuts)
25c94fc8 77{
78 // Default constructor
79 // Output slot #1 writes into a TList container
4c230f6d 80 Int_t nptbins=fAnalysisCuts->GetNPtBins();
81 Float_t *ptlim=fAnalysisCuts->GetPtBinLimits();
82 SetPtBins(nptbins,ptlim);
25c94fc8 83
84 DefineOutput(1,TList::Class()); //My private output
85
4c230f6d 86 DefineOutput(2,TList::Class());
25c94fc8 87}
88
89//________________________________________________________________________
4c230f6d 90void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Float_t* lim){
25c94fc8 91 // define pt bins for analysis
92 if(n>kMaxPtBins){
93 printf("Max. number of Pt bins = %d\n",kMaxPtBins);
94 fNPtBins=kMaxPtBins;
95 fPtLimits[0]=0.;
96 fPtLimits[1]=1.;
97 fPtLimits[2]=3.;
98 fPtLimits[3]=5.;
99 fPtLimits[4]=10.;
100 for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
101 }else{
102 fNPtBins=n;
103 for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
104 for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
105 }
106 if(fDebug > 1){
107 printf("Number of Pt bins = %d\n",fNPtBins);
108 for(Int_t i=0; i<fNPtBins+1; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);
109 }
110}
111//________________________________________________________________________
112AliAnalysisTaskSEDs::~AliAnalysisTaskSEDs()
113{
114 // Destructor
115 if (fOutput) {
116 delete fOutput;
117 fOutput = 0;
118 }
4c230f6d 119 if (fListCuts) {
120 delete fListCuts;
121 fListCuts = 0;
122 }
123
124 if (fProdCuts) {
125 delete fProdCuts;
126 fProdCuts = 0;
127 }
128 if (fAnalysisCuts) {
129 delete fAnalysisCuts;
130 fAnalysisCuts = 0;
25c94fc8 131 }
132}
133
134//________________________________________________________________________
135void AliAnalysisTaskSEDs::Init()
136{
137 // Initialization
138
139 if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
140
4c230f6d 141 fListCuts=new TList();
142 AliRDHFCutsDstoKKpi *production = new AliRDHFCutsDstoKKpi();
143 production=fProdCuts;
144 AliRDHFCutsDstoKKpi *analysis = new AliRDHFCutsDstoKKpi();
145 analysis=fAnalysisCuts;
146
147 fListCuts->Add(production);
148 fListCuts->Add(analysis);
149 PostData(2,fListCuts);
25c94fc8 150 return;
151}
152
153//________________________________________________________________________
154void AliAnalysisTaskSEDs::UserCreateOutputObjects()
155{
156 // Create the output container
157 //
158 if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
159
160 // Several histograms are more conveniently managed in a TList
161 fOutput = new TList();
162 fOutput->SetOwner();
163 fOutput->SetName("OutputHistos");
164
165 Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
a0bb2d06 166 Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);
167 if(nInvMassBins%2==1) nInvMassBins++;
168 Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
169 Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
170
25c94fc8 171 TString hisname;
172 Int_t index;
173 for(Int_t i=0;i<fNPtBins;i++){
174 index=GetHistoIndex(i);
175 hisname.Form("hMassAllPt%d",i);
a0bb2d06 176 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
25c94fc8 177 fMassHist[index]->Sumw2();
178 hisname.Form("hMassAllPt%dCuts",i);
a0bb2d06 179 fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
25c94fc8 180 fMassHistCuts[index]->Sumw2();
181 hisname.Form("hCosPAllPt%d",i);
182 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
183 fCosPHist[index]->Sumw2();
184 hisname.Form("hDLenAllPt%d",i);
185 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
186 fDLenHist[index]->Sumw2();
c2a0a9d5 187 hisname.Form("hSumd02AllPt%d",i);
188 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
189 fSumd02Hist[index]->Sumw2();
190 hisname.Form("hSigVertAllPt%d",i);
191 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
192 fSigVertHist[index]->Sumw2();
193 hisname.Form("hPtMaxAllPt%d",i);
194 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
195 fPtMaxHist[index]->Sumw2();
196 hisname.Form("hPtCandAllPt%d",i);
197 fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
198 fPtCandHist[index]->Sumw2();
199 hisname.Form("hDCAAllPt%d",i);
200 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
201 fDCAHist[index]->Sumw2();
202 hisname.Form("hPtProng0AllPt%d",i);
203 fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
204 fPtProng0Hist[index]->Sumw2();
205 hisname.Form("hPtProng1AllPt%d",i);
206 fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
207 fPtProng1Hist[index]->Sumw2();
208 hisname.Form("hPtProng2AllPt%d",i);
209 fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
210 fPtProng2Hist[index]->Sumw2();
211
25c94fc8 212 hisname.Form("hDalitzAllPt%d",i);
213 fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
214 fDalitz[index]->Sumw2();
215
216 index=GetSignalHistoIndex(i);
217 hisname.Form("hMassSigPt%d",i);
a0bb2d06 218 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
25c94fc8 219 fMassHist[index]->Sumw2();
220 hisname.Form("hMassSigPt%dCuts",i);
a0bb2d06 221 fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
25c94fc8 222 fMassHistCuts[index]->Sumw2();
223 hisname.Form("hCosPSigPt%d",i);
224 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
225 fCosPHist[index]->Sumw2();
226 hisname.Form("hDLenSigPt%d",i);
227 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
228 fDLenHist[index]->Sumw2();
c2a0a9d5 229 hisname.Form("hSumd02SigPt%d",i);
230 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
231 fSumd02Hist[index]->Sumw2();
232 hisname.Form("hSigVertSigPt%d",i);
233 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
234 fSigVertHist[index]->Sumw2();
235 hisname.Form("hPtMaxSigPt%d",i);
236 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
237 fPtMaxHist[index]->Sumw2();
238 hisname.Form("hPtCandSigPt%d",i);
239 fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
240 fPtCandHist[index]->Sumw2();
241 hisname.Form("hDCASigPt%d",i);
242 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
243 fDCAHist[index]->Sumw2();
244 hisname.Form("hPtProng0SigPt%d",i);
245 fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
246 fPtProng0Hist[index]->Sumw2();
247 hisname.Form("hPtProng1SigPt%d",i);
248 fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
249 fPtProng1Hist[index]->Sumw2();
250 hisname.Form("hPtProng2SigPt%d",i);
251 fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
252 fPtProng2Hist[index]->Sumw2();
25c94fc8 253 hisname.Form("hDalitzSigPt%d",i);
254 fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
255 fDalitz[index]->Sumw2();
256
25c94fc8 257 index=GetBackgroundHistoIndex(i);
c2a0a9d5 258 hisname.Form("hMassBkgPt%d",i);
a0bb2d06 259 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
25c94fc8 260 fMassHist[index]->Sumw2();
261 hisname.Form("hMassBkgPt%dCuts",i);
a0bb2d06 262 fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
25c94fc8 263 fMassHistCuts[index]->Sumw2();
264 hisname.Form("hCosPBkgPt%d",i);
265 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
266 fCosPHist[index]->Sumw2();
267 hisname.Form("hDLenBkgPt%d",i);
268 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
269 fDLenHist[index]->Sumw2();
c2a0a9d5 270 hisname.Form("hSumd02BkgPt%d",i);
271 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
272 fSumd02Hist[index]->Sumw2();
273 hisname.Form("hSigVertBkgPt%d",i);
274 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
275 fSigVertHist[index]->Sumw2();
276 hisname.Form("hPtMaxSigBkg%d",i);
277 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
278 fPtMaxHist[index]->Sumw2();
279 hisname.Form("hPtCandBkgPt%d",i);
280 fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
281 fPtCandHist[index]->Sumw2();
282 hisname.Form("hDCABkgPt%d",i);
283 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
284 fDCAHist[index]->Sumw2();
285 hisname.Form("hPtProng0BkgPt%d",i);
286 fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
287 fPtProng0Hist[index]->Sumw2();
288 hisname.Form("hPtProng1BkgPt%d",i);
289 fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
290 fPtProng1Hist[index]->Sumw2();
291 hisname.Form("hPtProng2BkgPt%d",i);
292 fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
293 fPtProng2Hist[index]->Sumw2();
25c94fc8 294 hisname.Form("hDalitzBkgPt%d",i);
295 fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
296 fDalitz[index]->Sumw2();
c2a0a9d5 297
298 index=GetReflSignalHistoIndex(i);
299 hisname.Form("hMassReflSigPt%d",i);
300 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
301 fMassHist[index]->Sumw2();
302 hisname.Form("hMassReflSigPt%dCuts",i);
303 fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
304 fMassHistCuts[index]->Sumw2();
305 hisname.Form("hCosPReflSigPt%d",i);
306 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
307 fCosPHist[index]->Sumw2();
308 hisname.Form("hDLenReflSigPt%d",i);
309 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
310 fDLenHist[index]->Sumw2();
311 hisname.Form("hSumd02ReflSigPt%d",i);
312 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
313 fSumd02Hist[index]->Sumw2();
314 hisname.Form("hSigVertReflSigPt%d",i);
315 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
316 fSigVertHist[index]->Sumw2();
317 hisname.Form("hPtMaxReflSigPt%d",i);
318 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
319 fPtMaxHist[index]->Sumw2();
320 hisname.Form("hPtCandReflSigPt%d",i);
321 fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
322 fPtCandHist[index]->Sumw2();
323 hisname.Form("hDCAReflSigPt%d",i);
324 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
325 fDCAHist[index]->Sumw2();
326 hisname.Form("hPtProng0ReflSigPt%d",i);
327 fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
328 fPtProng0Hist[index]->Sumw2();
329 hisname.Form("hPtProng1ReflSigPt%d",i);
330 fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
331 fPtProng1Hist[index]->Sumw2();
332 hisname.Form("hPtProng2ReflSigPt%d",i);
333 fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
334 fPtProng2Hist[index]->Sumw2();
335 hisname.Form("hDalitzReflSigPt%d",i);
336 fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
337 fDalitz[index]->Sumw2();
25c94fc8 338 }
339
c2a0a9d5 340 for(Int_t i=0; i<4*fNPtBins; i++){
25c94fc8 341 fOutput->Add(fMassHist[i]);
342 fOutput->Add(fMassHistCuts[i]);
343 fOutput->Add(fCosPHist[i]);
344 fOutput->Add(fDLenHist[i]);
c2a0a9d5 345 fOutput->Add(fSumd02Hist[i]);
346 fOutput->Add(fSigVertHist[i]);
347 fOutput->Add(fPtMaxHist[i]);
348 fOutput->Add(fPtCandHist[i]);
349 fOutput->Add(fDCAHist[i]);
350 fOutput->Add(fPtProng0Hist[i]);
351 fOutput->Add(fPtProng1Hist[i]);
352 fOutput->Add(fPtProng2Hist[i]);
25c94fc8 353 fOutput->Add(fDalitz[i]);
354 }
355
356 fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",4,-0.5,3.5);
357 fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",4,-0.5,3.5);
358 fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",4,-0.5,3.5);
c2a0a9d5 359 fChanHist[3] = new TH1F("hChanReflSig", "KKpi and piKK candidates",4,-0.5,3.5);
25c94fc8 360 fChanHistCuts[0] = new TH1F("hChanAllCuts", "KKpi and piKK candidates",4,-0.5,3.5);
361 fChanHistCuts[1] = new TH1F("hChanSigCuts", "KKpi and piKK candidates",4,-0.5,3.5);
362 fChanHistCuts[2] = new TH1F("hChanBkgCuts", "KKpi and piKK candidates",4,-0.5,3.5);
c2a0a9d5 363 fChanHistCuts[3] = new TH1F("hChanReflSigCuts", "KKpi and piKK candidates",4,-0.5,3.5);
364 for(Int_t i=0;i<4;i++){
25c94fc8 365 fChanHist[i]->Sumw2();
366 fChanHist[i]->SetMinimum(0);
367 fChanHistCuts[i]->Sumw2();
368 fChanHistCuts[i]->SetMinimum(0);
369 fOutput->Add(fChanHist[i]);
370 fOutput->Add(fChanHistCuts[i]);
371 }
372
373 fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-1.5,1.5);
374 fHistNEvents->Sumw2();
375 fHistNEvents->SetMinimum(0);
376 fOutput->Add(fHistNEvents);
377
378
379 return;
380}
381
382//________________________________________________________________________
383void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
384{
385 // Ds selection for current event, fill mass histos and selecetion variable histo
386 // separate signal and backgound if fReadMC is activated
387
388 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
389 fHistNEvents->Fill(0); // count event
390 // Post the data already here
391 PostData(1,fOutput);
392
393
394 TClonesArray *array3Prong = 0;
395 if(!aod && AODEvent() && IsStandardAOD()) {
396 // In case there is an AOD handler writing a standard AOD, use the AOD
397 // event in memory rather than the input (ESD) event.
398 aod = dynamic_cast<AliAODEvent*> (AODEvent());
399 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
400 // have to taken from the AOD event hold by the AliAODExtension
401 AliAODHandler* aodHandler = (AliAODHandler*)
402 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
403 if(aodHandler->GetExtensions()) {
404 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
405 AliAODEvent *aodFromExt = ext->GetAOD();
406 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
407 }
408 } else {
409 array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
410 }
411
412 if(!array3Prong) {
413 printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
414 return;
415 }
416
417
418 TClonesArray *arrayMC=0;
419 AliAODMCHeader *mcHeader=0;
420
421 // AOD primary vertex
422 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
423 // vtx1->Print();
424
425 // load MC particles
426 if(fReadMC){
427
428 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
429 if(!arrayMC) {
430 printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
431 return;
432 }
433
434 // load MC header
435 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
436 if(!mcHeader) {
437 printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
438 return;
439 }
440 }
441
442 Int_t n3Prong = array3Prong->GetEntriesFast();
443 if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
444
445
446 Int_t pdgDstoKKpi[3]={321,321,211};
25c94fc8 447 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
448 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
449
450
451 Bool_t unsetvtx=kFALSE;
452 if(!d->GetOwnPrimaryVtx()){
453 d->SetOwnPrimaryVtx(vtx1);
454 unsetvtx=kTRUE;
455 }
4c230f6d 456
457 Int_t retCodeProductionCuts=fProdCuts->IsSelected(d,AliRDHFCuts::kCandidate);
458 if(retCodeProductionCuts>0){
459 Int_t isKKpi=retCodeProductionCuts&1;
460 Int_t ispiKK=retCodeProductionCuts&2;
461// Int_t isPhi=retCodeProductionCuts&4;
462// Int_t isK0star=retCodeProductionCuts&8;
25c94fc8 463 Double_t ptCand = d->Pt();
4c230f6d 464 Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
465 Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kCandidate);
466 Int_t isKKpiAC=retCodeAnalysisCuts&1;
467 Int_t ispiKKAC=retCodeAnalysisCuts&2;
468// Int_t isPhiAC=retCodeAnalysisCuts&4;
469// Int_t isK0starAC=retCodeAnalysisCuts&8;
470
25c94fc8 471 Int_t labDs=-1;
472 if(fReadMC){
473 labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
474 }
475
476 Double_t dlen=d->DecayLength();
477 Double_t cosp=d->CosPointingAngle();
c2a0a9d5 478 Double_t pt0=d->PtProng(0);
479 Double_t pt1=d->PtProng(1);
480 Double_t pt2=d->PtProng(2);
481 Double_t sigvert=d->GetSigmaVert();
482 Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
483 Double_t dca=d->GetDCA();
484 Double_t ptmax=0;
485 for(Int_t i=0;i<3;i++){
486 if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
487 }
488
25c94fc8 489 Int_t index=GetHistoIndex(iPtBin);
490 Int_t type=0;
491 if(isKKpi) type+=1;
492 if(ispiKK) type+=2;
4c230f6d 493 Int_t typeAC=0;
494 if(isKKpiAC) typeAC+=1;
495 if(ispiKKAC) typeAC+=2;
25c94fc8 496 fCosPHist[index]->Fill(cosp);
497 fDLenHist[index]->Fill(dlen);
c2a0a9d5 498 fSigVertHist[index]->Fill(sigvert);
499 fSumd02Hist[index]->Fill(sumD02);
500 fPtMaxHist[index]->Fill(ptmax);
501 fPtCandHist[index]->Fill(ptCand);
502 fDCAHist[index]->Fill(dca);
25c94fc8 503 fChanHist[0]->Fill(type);
c2a0a9d5 504 fPtProng0Hist[index]->Fill(pt0);
505 fPtProng1Hist[index]->Fill(pt1);
506 fPtProng2Hist[index]->Fill(pt2);
507
4c230f6d 508 if(retCodeAnalysisCuts>0) fChanHistCuts[0]->Fill(typeAC);
25c94fc8 509 if(fReadMC){
510 if(labDs>=0) {
511 index=GetSignalHistoIndex(iPtBin);
25c94fc8 512 fChanHist[1]->Fill(type);
4c230f6d 513 if(retCodeAnalysisCuts>0) fChanHistCuts[1]->Fill(typeAC);
25c94fc8 514 }else{
515 index=GetBackgroundHistoIndex(iPtBin);
25c94fc8 516 fChanHist[2]->Fill(type);
4c230f6d 517 if(retCodeAnalysisCuts>0) fChanHistCuts[2]->Fill(typeAC);
25c94fc8 518 }
519 }
520 if(isKKpi){
521 index=GetHistoIndex(iPtBin);
522 Double_t invMass=d->InvMassDsKKpi();
523 fMassHist[index]->Fill(invMass);
524 Double_t mass01=d->InvMass2Prongs(0,1,321,321);
525 Double_t mass12=d->InvMass2Prongs(1,2,321,211);
526 fDalitz[index]->Fill(mass01,mass12);
4c230f6d 527 if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
25c94fc8 528 if(fReadMC){
c2a0a9d5 529 Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
530 AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
531 Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
532 if(labDs>=0){
533 if(pdgCode0==321) {
534 index=GetSignalHistoIndex(iPtBin);
535 }else{
536 index=GetReflSignalHistoIndex(iPtBin);
537 }
25c94fc8 538 }else{
539 index=GetBackgroundHistoIndex(iPtBin);
25c94fc8 540 }
c2a0a9d5 541 fMassHist[index]->Fill(invMass);
542 fCosPHist[index]->Fill(cosp);
543 fDLenHist[index]->Fill(dlen);
544 fSigVertHist[index]->Fill(sigvert);
545 fSumd02Hist[index]->Fill(sumD02);
546 fPtMaxHist[index]->Fill(ptmax);
547 fPtCandHist[index]->Fill(ptCand);
548 fDCAHist[index]->Fill(dca);
549 fPtProng0Hist[index]->Fill(pt0);
550 fPtProng1Hist[index]->Fill(pt1);
551 fPtProng2Hist[index]->Fill(pt2);
552 fDalitz[index]->Fill(mass01,mass12);
553 if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
25c94fc8 554 }
555 }
556 if(ispiKK){
557 index=GetHistoIndex(iPtBin);
558 Double_t invMass=d->InvMassDspiKK();
559 fMassHist[index]->Fill(invMass);
4c230f6d 560 if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
25c94fc8 561 if(fReadMC){
c2a0a9d5 562 Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
563 AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
564 Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
25c94fc8 565 if(labDs>=0) {
c2a0a9d5 566 if(pdgCode0==211) {
567 index=GetSignalHistoIndex(iPtBin);
568 }else{
569 index=GetReflSignalHistoIndex(iPtBin);
570 }
25c94fc8 571 }else{
572 index=GetBackgroundHistoIndex(iPtBin);
25c94fc8 573 }
c2a0a9d5 574 fMassHist[index]->Fill(invMass);
575 fCosPHist[index]->Fill(cosp);
576 fDLenHist[index]->Fill(dlen);
577 fSigVertHist[index]->Fill(sigvert);
578 fSumd02Hist[index]->Fill(sumD02);
579 fPtMaxHist[index]->Fill(ptmax);
580 fPtCandHist[index]->Fill(ptCand);
581 fDCAHist[index]->Fill(dca);
582 fPtProng0Hist[index]->Fill(pt0);
583 fPtProng1Hist[index]->Fill(pt1);
584 fPtProng2Hist[index]->Fill(pt2);
585 if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
25c94fc8 586 }
587 }
588 }
589 if(unsetvtx) d->UnsetOwnPrimaryVtx();
590 }
591
592
593 PostData(1,fOutput);
594 return;
595}
596
597//_________________________________________________________________
598
599void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
600{
601 // Terminate analysis
602 //
603 if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
604 fOutput = dynamic_cast<TList*> (GetOutputData(1));
605 if (!fOutput) {
606 printf("ERROR: fOutput not available\n");
607 return;
608 }
609 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
610 fChanHist[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAll"));
611 fChanHist[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSig"));
612 fChanHist[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkg"));
c2a0a9d5 613 fChanHist[3] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanReflSig"));
25c94fc8 614 fChanHistCuts[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAllCuts"));
615 fChanHistCuts[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSigCuts"));
616 fChanHistCuts[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkgCuts"));
c2a0a9d5 617 fChanHistCuts[3] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanReflSigCuts"));
25c94fc8 618
619
620 TString hisname;
621 Int_t index;
622 for(Int_t i=0;i<fNPtBins;i++){
623
624 index=GetHistoIndex(i);
625 hisname.Form("hMassAllPt%d",i);
626 fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
627 hisname.Form("hMassAllPt%dCuts",i);
628 fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
629 hisname.Form("hCosPAllPt%d",i);
630 fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
631 hisname.Form("hDLenAllPt%d",i);
632 fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
c2a0a9d5 633 hisname.Form("hSumd02AllPt%d",i);
634 fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
635 hisname.Form("hSigVertAllPt%d",i);
636 fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
637 hisname.Form("hPtMaxAllPt%d",i);
638 fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
639 hisname.Form("hPtCandAllPt%d",i);
640 fPtCandHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
641 hisname.Form("hDCAAllPt%d",i);
642 fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
643 hisname.Form("hPtProng0AllPt%d",i);
644 fPtProng0Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
645 hisname.Form("hPtProng1AllPt%d",i);
646 fPtProng1Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
647 hisname.Form("hPtProng2AllPt%d",i);
648 fPtProng2Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
25c94fc8 649 hisname.Form("hDalitzAllPt%d",i);
650 fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
651
652 index=GetSignalHistoIndex(i);
653 hisname.Form("hMassSigPt%d",i);
654 fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
655 hisname.Form("hMassSigPt%dCuts",i);
656 fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
657 hisname.Form("hCosPSigPt%d",i);
658 fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
659 hisname.Form("hDLenSigPt%d",i);
660 fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
c2a0a9d5 661 hisname.Form("hSumd02SigPt%d",i);
662 fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
663 hisname.Form("hSigVertSigPt%d",i);
664 fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
665 hisname.Form("hPtMaxSigPt%d",i);
666 fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
667 hisname.Form("hPtCandSigPt%d",i);
668 fPtCandHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
669 hisname.Form("hDCASigPt%d",i);
670 fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
671 hisname.Form("hPtProng0SigPt%d",i);
672 fPtProng0Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
673 hisname.Form("hPtProng1SigPt%d",i);
674 fPtProng1Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
675 hisname.Form("hPtProng2SigPt%d",i);
676 fPtProng2Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
25c94fc8 677 hisname.Form("hDalitzSigPt%d",i);
678 fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
679
680 index=GetBackgroundHistoIndex(i);
681 hisname.Form("hMassBkgPt%d",i);
682 fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
683 hisname.Form("hMassBkgPt%dCuts",i);
684 fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
685 hisname.Form("hCosPBkgPt%d",i);
686 fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
687 hisname.Form("hDLenBkgPt%d",i);
688 fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
c2a0a9d5 689 hisname.Form("hSumd02BkgPt%d",i);
690 fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
691 hisname.Form("hSigVertBkgPt%d",i);
692 fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
693 hisname.Form("hPtMaxBkgPt%d",i);
694 fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
695 hisname.Form("hPtCandBkgPt%d",i);
696 fPtCandHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
697 hisname.Form("hDCABkgPt%d",i);
698 fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
699 hisname.Form("hPtProng0BkgPt%d",i);
700 fPtProng0Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
701 hisname.Form("hPtProng1BkgPt%d",i);
702 fPtProng1Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
703 hisname.Form("hPtProng2BkgPt%d",i);
704 fPtProng2Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
25c94fc8 705 hisname.Form("hDalitzBkgPt%d",i);
706 fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
707
c2a0a9d5 708 index=GetReflSignalHistoIndex(i);
709 hisname.Form("hMassReflSigPt%d",i);
710 fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
711 hisname.Form("hMassReflSigPt%dCuts",i);
712 fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
713 hisname.Form("hCosPReflSigPt%d",i);
714 fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
715 hisname.Form("hDLenReflSigPt%d",i);
716 fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
717 hisname.Form("hSumd02ReflSigPt%d",i);
718 fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
719 hisname.Form("hSigVertReflSigPt%d",i);
720 fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
721 hisname.Form("hPtMaxReflSigPt%d",i);
722 fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
723 hisname.Form("hPtCandReflSigPt%d",i);
724 fPtCandHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
725 hisname.Form("hDCAReflSigPt%d",i);
726 fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
727 hisname.Form("hPtProng0ReflSigPt%d",i);
728 fPtProng0Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
729 hisname.Form("hPtProng1ReflSigPt%d",i);
730 fPtProng1Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
731 hisname.Form("hPtProng2ReflSigPt%d",i);
732 fPtProng2Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
733 hisname.Form("hDalitzReflSigPt%d",i);
734 fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
735
25c94fc8 736 }
737 return;
738}