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