AliPhysics  29d4213 (29d4213)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
AliAnalysisTaskSEDs.cxx
Go to the documentation of this file.
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 
19 // //
20 // Analysis task to produce Ds candidates mass spectra //
21 // Origin: F.Prino, Torino, prino@to.infn.it //
22 // //
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 <THnSparse.h>
32 #include <TDatabasePDG.h>
33 #include <Riostream.h>
34 
35 #include "AliAnalysisManager.h"
36 #include "AliAODHandler.h"
37 #include "AliAODEvent.h"
38 #include "AliAODVertex.h"
39 #include "AliAODTrack.h"
40 #include "AliAODMCHeader.h"
41 #include "AliAODMCParticle.h"
42 #include "AliAODRecoDecay.h"
44 #include "AliAnalysisVertexingHF.h"
45 #include "AliRDHFCutsDstoKKpi.h"
46 #include "AliAnalysisTaskSE.h"
48 #include "AliAnalysisTaskSEDs.h"
49 #include "AliVertexingHFUtils.h"
50 
54 
55 //________________________________________________________________________
57 AliAnalysisTaskSE(),
58 fOutput(0),
59 fHistNEvents(0),
60 fPtVsMass(0),
61 fPtVsMassPhi(0),
62 fPtVsMassK0st(0),
63 fYVsPt(0),
64 fYVsPtSig(0),
65 fNtupleDs(0),
66 fFillNtuple(0),
67 fReadMC(kFALSE),
68 fWriteOnlySignal(kFALSE),
69 fDoCutVarHistos(kTRUE),
70 fUseSelectionBit(kFALSE),
71 fFillSparse(kTRUE),
72 fNPtBins(0),
73 fListCuts(0),
74 fMassRange(0.8),
75 fMassBinSize(0.002),
76 fCounter(0),
77 fAnalysisCuts(0),
78 fnSparse(0),
79 fnSparseIP(0)
80 {
82 
83  for(Int_t i=0;i<3;i++){
84  fHistCentrality[i]=0;
86  }
87  for(Int_t i=0;i<4;i++) {
88  fChanHist[i]=0;
89  }
90  for(Int_t i=0;i<4*kMaxPtBins;i++){
91  fPtCandHist[i]=0;
92  fMassHist[i]=0;
93  fMassHistPhi[i]=0;
94  fMassHistK0st[i]=0;
95  fCosPHist[i]=0;
96  fDLenHist[i]=0;
97  fSumd02Hist[i]=0;
98  fSigVertHist[i]=0;
99  fPtMaxHist[i]=0;
100  fDCAHist[i]=0;
101  fPtProng0Hist[i]=0;
102  fPtProng1Hist[i]=0;
103  fPtProng2Hist[i]=0;
104  fDalitz[i]=0;
105  fDalitzPhi[i]=0;
106  fDalitzK0st[i]=0;
107  }
108  for(Int_t i=0;i<kMaxPtBins;i++){
109  fMassHistKK[i]=0;
110  fMassHistKpi[i]=0;
111  }
112  for(Int_t i=0;i<kMaxPtBins+1;i++){
113  fPtLimits[i]=0;
114  }
115  for (Int_t i=0; i<4; i++) {
116  fnSparseMC[i]=0;
117  }
118 }
119 
120 //________________________________________________________________________
121 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name,AliRDHFCutsDstoKKpi* analysiscuts,Int_t fillNtuple):
122 AliAnalysisTaskSE(name),
123 fOutput(0),
124 fHistNEvents(0),
125 fPtVsMass(0),
126 fPtVsMassPhi(0),
127 fPtVsMassK0st(0),
128 fYVsPt(0),
129 fYVsPtSig(0),
130 fNtupleDs(0),
131 fFillNtuple(fillNtuple),
132 fReadMC(kFALSE),
133 fWriteOnlySignal(kFALSE),
134 fDoCutVarHistos(kTRUE),
135 fUseSelectionBit(kFALSE),
136 fFillSparse(kTRUE),
137 fNPtBins(0),
138 fListCuts(0),
139 fMassRange(0.8),
140 fMassBinSize(0.002),
141 fCounter(0),
142 fAnalysisCuts(analysiscuts),
143 fnSparse(0),
144 fnSparseIP(0)
145 {
148 
149  for(Int_t i=0;i<3;i++){
150  fHistCentrality[i]=0;
151  fHistCentralityMult[i]=0;
152  }
153  for(Int_t i=0;i<4;i++) {
154  fChanHist[i]=0;
155  }
156  for(Int_t i=0;i<4*kMaxPtBins;i++){
157  fPtCandHist[i]=0;
158  fMassHist[i]=0;
159  fMassHistPhi[i]=0;
160  fMassHistK0st[i]=0;
161  fCosPHist[i]=0;
162  fDLenHist[i]=0;
163  fSumd02Hist[i]=0;
164  fSigVertHist[i]=0;
165  fPtMaxHist[i]=0;
166  fDCAHist[i]=0;
167  fPtProng0Hist[i]=0;
168  fPtProng1Hist[i]=0;
169  fPtProng2Hist[i]=0;
170  fDalitz[i]=0;
171  fDalitzPhi[i]=0;
172  fDalitzK0st[i]=0;
173  }
174  for(Int_t i=0;i<kMaxPtBins;i++){
175  fMassHistKK[i]=0;
176  fMassHistKpi[i]=0;
177  }
178  for(Int_t i=0;i<kMaxPtBins+1;i++){
179  fPtLimits[i]=0;
180  }
181 
182  for (Int_t i=0; i<4; i++) {
183  fnSparseMC[i]=0;
184  }
185 
187  Float_t *ptlim=fAnalysisCuts->GetPtBinLimits();
188  SetPtBins(nptbins,ptlim);
189 
190  DefineOutput(1,TList::Class()); //My private output
191 
192  DefineOutput(2,TList::Class());
193 
194  DefineOutput(3,AliNormalizationCounter::Class());
195 
196  if(fFillNtuple>0){
197  // Output slot #4 writes into a TNtuple container
198  DefineOutput(4,TNtuple::Class()); //My private output
199  }
200 
201 }
202 
203 //________________________________________________________________________
204 void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Float_t* lim){
206  if(n>kMaxPtBins){
207  printf("Max. number of Pt bins = %d\n",kMaxPtBins);
209  fPtLimits[0]=0.;
210  fPtLimits[1]=1.;
211  fPtLimits[2]=3.;
212  fPtLimits[3]=5.;
213  fPtLimits[4]=10.;
214  for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
215  }else{
216  fNPtBins=n;
217  for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
218  for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
219  }
220  if(fDebug > 1){
221  printf("Number of Pt bins = %d\n",fNPtBins);
222  for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);
223  }
224 }
225 //________________________________________________________________________
227 {
228  // Destructor
229  if(fOutput && !fOutput->IsOwner()){
230  delete fHistNEvents;
231  for(Int_t i=0;i<4;i++){
232  delete fChanHist[i];
233  }
234  for(Int_t i=0;i<4*fNPtBins;i++){
235  delete fMassHist[i];
236  delete fMassHistPhi[i];
237  delete fMassHistK0st[i];
238  delete fCosPHist[i];
239  delete fDLenHist[i];
240  delete fSumd02Hist[i];
241  delete fSigVertHist[i];
242  delete fPtMaxHist[i];
243  delete fPtCandHist[i];
244  delete fDCAHist[i];
245  delete fPtProng0Hist[i];
246  delete fPtProng1Hist[i];
247  delete fPtProng2Hist[i];
248  delete fDalitz[i];
249  delete fDalitzPhi[i];
250  delete fDalitzK0st[i];
251  }
252  for(Int_t i=0;i<fNPtBins;i++){
253  delete fMassHistKK[i];
254  delete fMassHistKpi[i];
255  }
256  delete fPtVsMass;
257  delete fPtVsMassPhi;
258  delete fPtVsMassK0st;
259  delete fYVsPt;
260  delete fYVsPtSig;
261  for(Int_t i=0;i<3;i++){
262  delete fHistCentrality[i];
263  delete fHistCentralityMult[i];
264  }
265 
266  delete fnSparse;
267  delete fnSparseIP;
268  for (Int_t i=0; i<4; i++) {
269  delete fnSparseMC[i];
270  }
271  }
272  delete fOutput;
273  delete fListCuts;
274  delete fNtupleDs;
275  delete fCounter;
276  delete fAnalysisCuts;
277 
278 }
279 
280 //________________________________________________________________________
282 {
284 
285  if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
286 
287  fListCuts=new TList();
288  fListCuts->SetOwner();
289  fListCuts->SetName("CutObjects");
290 
292  analysis->SetName("AnalysisCuts");
293 
294  fListCuts->Add(analysis);
295  PostData(2,fListCuts);
296  return;
297 }
298 
299 //________________________________________________________________________
301 {
303  //
304  if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
305 
306  // Several histograms are more conveniently managed in a TList
307  fOutput = new TList();
308  fOutput->SetOwner();
309  fOutput->SetName("OutputHistos");
310 
311  fHistNEvents = new TH1F("hNEvents", "number of events ",11,-0.5,10.5);
312  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
313  fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
314  fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
315  fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to not reco vertex");
316  fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected for contr vertex");
317  fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for vertex out of accept");
318  fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for pileup events");
319  fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
320  fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of 3 prong candidates");
321  fHistNEvents->GetXaxis()->SetBinLabel(10,"no. of Ds after filtering cuts");
322  fHistNEvents->GetXaxis()->SetBinLabel(11,"no. of Ds after selection cuts");
323 
324  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
325 
326  fHistNEvents->Sumw2();
327  fHistNEvents->SetMinimum(0);
328  fOutput->Add(fHistNEvents);
329 
330  fHistCentrality[0]=new TH1F("hCentr","centrality",10000,0.,100.);
331  fHistCentrality[1]=new TH1F("hCentr(selectedCent)","centrality(selectedCent)",10000,0.,100.);
332  fHistCentrality[2]=new TH1F("hCentr(OutofCent)","centrality(OutofCent)",10000,0.,100.);
333  fHistCentralityMult[0]=new TH2F("hCentrMult","centrality vs mult",100,0.5,30000.5,40,0.,100.);
334  fHistCentralityMult[1]=new TH2F("hCentrMult(selectedCent)","centrality vs mult(selectedCent)",100,0.5,30000.5,40,0.,100.);
335  fHistCentralityMult[2]=new TH2F("hCentrMult(OutofCent)","centrality vs mult(OutofCent)",100,0.5,30000.5,40,0.,100.);
336  for(Int_t i=0;i<3;i++){
337  fHistCentrality[i]->Sumw2();
338  fOutput->Add(fHistCentrality[i]);
339  fHistCentralityMult[i]->Sumw2();
340  fOutput->Add(fHistCentralityMult[i]);
341  }
342 
343  Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
344 
345  Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);
346  if(nInvMassBins%2==1) nInvMassBins++;
347  // Double_t minMass=massDs-1.0*nInvMassBins*fMassBinSize;
348  Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
349  // Double_t maxMass=massDs+1.0*nInvMassBins*fMassBinSize;
350  Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
351 
352  TString hisname;
353  TString htype;
354  Int_t index;
355  for(Int_t iType=0; iType<4; iType++){
356  for(Int_t i=0;i<fNPtBins;i++){
357  if(iType==0){
358  htype="All";
359  index=GetHistoIndex(i);
360  }else if(iType==1){
361  htype="Sig";
362  index=GetSignalHistoIndex(i);
363  }else if(iType==2){
364  htype="Bkg";
365  index=GetBackgroundHistoIndex(i);
366  }else{
367  htype="ReflSig";
368  index=GetReflSignalHistoIndex(i);
369  }
370  hisname.Form("hMass%sPt%d",htype.Data(),i);
371  fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
372  fMassHist[index]->Sumw2();
373  hisname.Form("hMass%sPt%dphi",htype.Data(),i);
374  fMassHistPhi[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
375  fMassHistPhi[index]->Sumw2();
376  hisname.Form("hMass%sPt%dk0st",htype.Data(),i);
377  fMassHistK0st[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
378  fMassHistK0st[index]->Sumw2();
379  hisname.Form("hCosP%sPt%d",htype.Data(),i);
380  fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
381  fCosPHist[index]->Sumw2();
382  hisname.Form("hDLen%sPt%d",htype.Data(),i);
383  fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
384  fDLenHist[index]->Sumw2();
385  hisname.Form("hSumd02%sPt%d",htype.Data(),i);
386  fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
387  fSumd02Hist[index]->Sumw2();
388  hisname.Form("hSigVert%sPt%d",htype.Data(),i);
389  fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
390  fSigVertHist[index]->Sumw2();
391  hisname.Form("hPtMax%sPt%d",htype.Data(),i);
392  fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
393  fPtMaxHist[index]->Sumw2();
394  hisname.Form("hPtCand%sPt%d",htype.Data(),i);
395  fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
396  fPtCandHist[index]->Sumw2();
397  hisname.Form("hDCA%sPt%d",htype.Data(),i);
398  fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
399  fDCAHist[index]->Sumw2();
400  hisname.Form("hPtProng0%sPt%d",htype.Data(),i);
401  fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
402  fPtProng0Hist[index]->Sumw2();
403  hisname.Form("hPtProng1%sPt%d",htype.Data(),i);
404  fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
405  fPtProng1Hist[index]->Sumw2();
406  hisname.Form("hPtProng2%sPt%d",htype.Data(),i);
407  fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
408  fPtProng2Hist[index]->Sumw2();
409  hisname.Form("hDalitz%sPt%d",htype.Data(),i);
410  fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
411  fDalitz[index]->Sumw2();
412  hisname.Form("hDalitz%sPt%dphi",htype.Data(),i);
413  fDalitzPhi[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
414  fDalitzPhi[index]->Sumw2();
415  hisname.Form("hDalitz%sPt%dk0st",htype.Data(),i);
416  fDalitzK0st[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
417  fDalitzK0st[index]->Sumw2();
418  }
419  }
420 
421  for(Int_t i=0; i<4*fNPtBins; i++){
422  fOutput->Add(fMassHist[i]);
423  fOutput->Add(fMassHistPhi[i]);
424  fOutput->Add(fMassHistK0st[i]);
425  fOutput->Add(fPtCandHist[i]);
426  if(fDoCutVarHistos){
427  fOutput->Add(fCosPHist[i]);
428  fOutput->Add(fDLenHist[i]);
429  fOutput->Add(fSumd02Hist[i]);
430  fOutput->Add(fSigVertHist[i]);
431  fOutput->Add(fPtMaxHist[i]);
432  fOutput->Add(fDCAHist[i]);
433  fOutput->Add(fPtProng0Hist[i]);
434  fOutput->Add(fPtProng1Hist[i]);
435  fOutput->Add(fPtProng2Hist[i]);
436  fOutput->Add(fDalitz[i]);
437  fOutput->Add(fDalitzPhi[i]);
438  fOutput->Add(fDalitzK0st[i]);
439  }
440  }
441 
442  fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",64,-0.5,63.5);
443  fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",64,-0.5,63.5);
444  fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",64,-0.5,63.5);
445  fChanHist[3] = new TH1F("hChanReflSig", "KKpi and piKK candidates",64,-0.5,63.5);
446  for(Int_t i=0;i<4;i++){
447  fChanHist[i]->Sumw2();
448  fChanHist[i]->SetMinimum(0);
449  fOutput->Add(fChanHist[i]);
450  }
451 
452 
453  fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nInvMassBins,minMass,maxMass,40,0.,20.);
454  fPtVsMassPhi=new TH2F("hPtVsMassPhi","PtVsMass (phi selection)",nInvMassBins,minMass,maxMass,200,0.,20.);
455  fPtVsMassK0st=new TH2F("hPtVsMassK0st","PtVsMass (K0* selection)",nInvMassBins,minMass,maxMass,200,0.,20.);
456  fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
457  fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
458 
459  for(Int_t i=0;i<fNPtBins;i++){
460  hisname.Form("hMassKKPt%d",i);
461  fMassHistKK[i]=new TH1F(hisname.Data(),hisname.Data(),200,0.95,1.35);
462  fMassHistKK[i]->Sumw2();
463  fOutput->Add(fMassHistKK[i]);
464  hisname.Form("hMassKpiPt%d",i);
465  fMassHistKpi[i]=new TH1F(hisname.Data(),hisname.Data(),200,0.7,1.1);
466  fMassHistKpi[i]->Sumw2();
467  fOutput->Add(fMassHistKpi[i]);
468  }
469 
470  fOutput->Add(fPtVsMass);
471  fOutput->Add(fPtVsMassPhi);
472  fOutput->Add(fPtVsMassK0st);
473  fOutput->Add(fYVsPt);
474  fOutput->Add(fYVsPtSig);
475 
476  const Int_t nvarsReco = 16;
477  const Int_t nvarsAcc = 2;
478  Int_t nBinsReco[nvarsReco] = {500,20,60,1000,1000,500,500,50,50,100,100,100,400,400,400,400};
479  Int_t nBinsAcc[nvarsAcc] = {100,20};
480  Double_t xminAcc[nvarsAcc] = {0.,-1.};
481  Double_t xmaxAcc[nvarsAcc] = {100.,1.};
482  Double_t xminReco[nvarsReco] = {1.5,0.,0.,0.,0.,0.,0.,0.5,0.5,0.,0.,0.,-10.,-10.,-10.,-10.};
483  Double_t xmaxReco[nvarsReco] = {2.5,20.,0.03,5.,5.,500.,500.,1.,1.,0.1,1.,1.,10.,10.,10.,10.};
484  TString axis[nvarsReco] = {"invMassAllPhi","p_{T}","#Delta Mass(KK)","dlen","dlen_{xy}","normdl","normdl_{xy}","cosP","cosP_{xy}",
485  "sigVert","cosPiDs","|cosPiKPhi^{3}|","normIP","normIP_kaon1","normIP_kaon2","normIP_pion"};
486 
487  const Int_t nvarsIP = 6;
488  Int_t nBinsIP[nvarsIP] = {20,400,400,400,400,3};
489  Double_t xminIP[nvarsIP] = {0.,-10.,-10.,-10.,-10.,0.};
490  Double_t xmaxIP[nvarsIP] = {20.,10.,10.,10.,10.,3.};
491  TString axisIP[nvarsIP] = {"motherPt","maxNormImp","IP0","IP1","IP2","candType"};
492 
493  if(fFillSparse) {
494 
495  if(fReadMC) {
496  TString label[nvarsAcc] = {"fromC","fromB"};
497  for (Int_t i=0; i<2; i++) {
498  fnSparseMC[i] = new THnSparseF(Form("fnSparseAcc_%s",label[i].Data()),Form("MC nSparse (Acc.Step)- %s",label[i].Data()),
499  nvarsAcc, nBinsAcc, xminAcc, xmaxAcc);
500  fnSparseMC[i]->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
501  fnSparseMC[i]->GetAxis(1)->SetTitle("y");
502  fOutput->Add(fnSparseMC[i]);
503  }
504  for (Int_t i=2; i<4; i++) {
505  fnSparseMC[i] = new THnSparseF(Form("fnSparseReco_%s",label[i-2].Data()),Form("MC nSparse (Reco Step)- %s",label[i-2].Data()),
506  nvarsReco, nBinsReco, xminReco, xmaxReco);
507  for (Int_t j=0; j<nvarsReco; j++) {
508  fnSparseMC[i]->GetAxis(j)->SetTitle(Form("%s",axis[j].Data()));
509  }
510  fOutput->Add(fnSparseMC[i]);
511  }
512 
513  fnSparseIP = new THnSparseF("fnSparseIP","nSparseIP", nvarsIP, nBinsIP, xminIP, xmaxIP);
514  for (Int_t j=0; j<nvarsIP; j++) {
515  fnSparseIP->GetAxis(j)->SetTitle(Form("%s",axisIP[j].Data()));
516  }
517  fnSparseIP->GetAxis(5)->SetTitle("candType (0.5=bkg; 1.5=prompt; 2.5=FD)");
518  fOutput->Add(fnSparseIP);
519  }
520  else {
521  fnSparse = new THnSparseF("fnSparse","nSparse", nvarsReco, nBinsReco, xminReco, xmaxReco);
522  for (Int_t j=0; j<nvarsReco; j++) {
523  fnSparse->GetAxis(j)->SetTitle(Form("%s",axis[j].Data()));
524  }
525  fOutput->Add(fnSparse);
526  }
527  }
528 
529  //Counter for Normalization
530  fCounter = new AliNormalizationCounter("NormalizationCounter");
531  fCounter->Init();
532 
533  PostData(1,fOutput);
534  PostData(3,fCounter);
535 
536  if(fFillNtuple>0){
537  OpenFile(4); // 4 is the slot number of the ntuple
538 
539  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");
540 
541  }
542 
543 
544  return;
545 }
546 
547 //________________________________________________________________________
548 void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
549 {
552 
553  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
554 
555  TClonesArray *array3Prong = 0;
556  if(!aod && AODEvent() && IsStandardAOD()) {
557  // In case there is an AOD handler writing a standard AOD, use the AOD
558  // event in memory rather than the input (ESD) event.
559  aod = dynamic_cast<AliAODEvent*> (AODEvent());
560  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
561  // have to taken from the AOD event hold by the AliAODExtension
562  AliAODHandler* aodHandler = (AliAODHandler*)
563  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
564  if(aodHandler->GetExtensions()) {
565  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
566  AliAODEvent *aodFromExt = ext->GetAOD();
567  array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
568  }
569  } else if(aod) {
570  array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
571  }
572 
573  if(!aod || !array3Prong) {
574  printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
575  return;
576  }
577 
578 
579  // fix for temporary bug in ESDfilter
580  // the AODs with null vertex pointer didn't pass the PhysSel
581  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
582 
583 
584  fHistNEvents->Fill(0); // count event
585  // Post the data already here
586  PostData(1,fOutput);
587 
589 
590 
591  Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
592  Float_t ntracks=aod->GetNumberOfTracks();
593  Float_t evCentr=fAnalysisCuts->GetCentrality(aod);
594 
595  fHistCentrality[0]->Fill(evCentr);
596  fHistCentralityMult[0]->Fill(ntracks,evCentr);
603  fHistNEvents->Fill(7);
604  fHistCentrality[2]->Fill(evCentr);
605  fHistCentralityMult[2]->Fill(ntracks,evCentr);
606  }
607 
608  Float_t centrality=fAnalysisCuts->GetCentrality(aod);
609  Int_t runNumber=aod->GetRunNumber();
610 
611 
612 
613  TClonesArray *arrayMC=0;
614  AliAODMCHeader *mcHeader=0;
615 
616  // AOD primary vertex
617  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
618  // vtx1->Print();
619 
620  // load MC particles
621  if(fReadMC){
622 
623  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
624  if(!arrayMC) {
625  printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
626  return;
627  }
628 
629  // load MC header
630  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
631  if(!mcHeader) {
632  printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
633  return;
634  }
635  }
636 
637 
638  if(fReadMC && fFillSparse){
639  if(aod->GetTriggerMask()==0 && (runNumber>=195344 && runNumber<=195677))
640  // protection for events with empty trigger mask in p-Pb
641  return;
643  // events not passing the centrality selection can be removed immediately.
644  return;
645  Double_t zMCVertex = mcHeader->GetVtxZ();
646  if (TMath::Abs(zMCVertex) > fAnalysisCuts->GetMaxVtxZ())
647  return;
648  FillMCGenAccHistos(arrayMC, mcHeader);
649  }
650 
651  if(!isEvSel)return;
652  fHistNEvents->Fill(1);
653  fHistCentrality[1]->Fill(evCentr);
654  fHistCentralityMult[1]->Fill(ntracks,evCentr);
655 
656  Int_t n3Prong = array3Prong->GetEntriesFast();
657  if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
658 
659 
660  Int_t pdgDstoKKpi[3]={321,321,211};
661  Int_t nSelected=0;
662  Int_t nFiltered=0;
663  Double_t massPhi=TDatabasePDG::Instance()->GetParticle(333)->Mass();
664 
665  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
666 
667  AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
668  fHistNEvents->Fill(8);
669 
671  continue;
672  }
673  nFiltered++;
674  fHistNEvents->Fill(9);
675 
676  Bool_t unsetvtx=kFALSE;
677  if(!d->GetOwnPrimaryVtx()){
678  d->SetOwnPrimaryVtx(vtx1);
679  unsetvtx=kTRUE;
680  }
681 
682  Bool_t recVtx=kFALSE;
683  AliAODVertex *origownvtx=0x0;
684 
685  Double_t ptCand = d->Pt();
686  Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
687  Double_t rapid=d->YDs();
688  fYVsPt->Fill(ptCand,rapid);
689  Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);
690 
691  if(!isFidAcc) continue;
692 
693  Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
694  Int_t retCodeNoRes=retCodeAnalysisCuts;
695  Bool_t origRes=fAnalysisCuts->IsCutOnResonancesApplied();
696  if(origRes){
698  retCodeNoRes=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
700  }
701  Double_t massKK = 0.;
702 
703  if(retCodeNoRes&1){ //KKpi
704  massKK=d->InvMass2Prongs(0,1,321,321);
705  Double_t massKp=d->InvMass2Prongs(1,2,321,211);
706  fMassHistKK[iPtBin]->Fill(massKK);
707  fMassHistKpi[iPtBin]->Fill(massKp);
708  }
709  if(retCodeNoRes&2){ //piKK
710  massKK=d->InvMass2Prongs(1,2,321,321);
711  Double_t massKp=d->InvMass2Prongs(0,1,211,321);
712  fMassHistKK[iPtBin]->Fill(massKK);
713  fMassHistKpi[iPtBin]->Fill(massKp);
714  }
715 
716  Int_t isKKpi=retCodeAnalysisCuts&1;
717  Int_t ispiKK=retCodeAnalysisCuts&2;
718  Int_t isPhiKKpi=retCodeAnalysisCuts&4;
719  Int_t isPhipiKK=retCodeAnalysisCuts&8;
720  Int_t isK0starKKpi=retCodeAnalysisCuts&16;
721  Int_t isK0starpiKK=retCodeAnalysisCuts&32;
722 
723  if(retCodeAnalysisCuts<=0) continue;
725  if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
726  if(fAnalysisCuts->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
727  else fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
728  }
729 
730 
731  fHistNEvents->Fill(10);
732  nSelected++;
733 
734  Int_t index=GetHistoIndex(iPtBin);
735  fPtCandHist[index]->Fill(ptCand);
736 
737 
738  Double_t weightKKpi=1.;
739  Double_t weightpiKK=1.;
741  weightKKpi=fAnalysisCuts->GetWeightForKKpi();
742  weightpiKK=fAnalysisCuts->GetWeightForpiKK();
743  if(weightKKpi>1. || weightKKpi<0.) weightKKpi=0.;
744  if(weightpiKK>1. || weightpiKK<0.) weightpiKK=0.;
745  }
746 
747  fChanHist[0]->Fill(retCodeAnalysisCuts);
748 
749 
750  Double_t invMass = 0.;
751  Int_t indexMCKKpi=-1;
752  Int_t indexMCpiKK=-1;
753  Int_t labDs=-1;
754  Int_t pdgCode0=-999;
755  Int_t isMCSignal=-1;
756 
757 
758  if(fReadMC){
759 
760  labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
761  if(labDs>=0){
762  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
763  AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
764  pdgCode0=TMath::Abs(p->GetPdgCode());
765 
766  if(isKKpi){
767  if(pdgCode0==321) {
768  indexMCKKpi=GetSignalHistoIndex(iPtBin);
769  fYVsPtSig->Fill(ptCand,rapid);
770  fChanHist[1]->Fill(retCodeAnalysisCuts);
771  isMCSignal=1;
772  }else{
773  indexMCKKpi=GetReflSignalHistoIndex(iPtBin);
774  fChanHist[3]->Fill(retCodeAnalysisCuts);
775  isMCSignal=0;
776  }
777  }
778  if(ispiKK){
779  if(pdgCode0==211) {
780  indexMCpiKK=GetSignalHistoIndex(iPtBin);
781  fYVsPtSig->Fill(ptCand,rapid);
782  fChanHist[1]->Fill(retCodeAnalysisCuts);
783  isMCSignal=1;
784  }else{
785  indexMCpiKK=GetReflSignalHistoIndex(iPtBin);
786  fChanHist[3]->Fill(retCodeAnalysisCuts);
787  isMCSignal=0;
788  }
789  }
790  }else{
791  indexMCpiKK=GetBackgroundHistoIndex(iPtBin);
792  indexMCKKpi=GetBackgroundHistoIndex(iPtBin);
793  fChanHist[2]->Fill(retCodeAnalysisCuts);
794  }
795  }
796 
797  Double_t candType = 0.5; //for bkg
798 
799  if(isKKpi){
800  invMass=d->InvMassDsKKpi();
801  fMassHist[index]->Fill(invMass,weightKKpi);
802  fPtVsMass->Fill(invMass,ptCand,weightKKpi);
803  if(isPhiKKpi){
804  fMassHistPhi[index]->Fill(invMass,weightKKpi);
805  fPtVsMassPhi->Fill(invMass,ptCand,weightKKpi);
806  }
807  if(isK0starKKpi){
808  fMassHistK0st[index]->Fill(invMass,weightKKpi);
809  fPtVsMassK0st->Fill(invMass,ptCand,weightKKpi);
810  }
811  if(fReadMC && indexMCKKpi!=-1){
812  fMassHist[indexMCKKpi]->Fill(invMass,weightKKpi);
813  if(isPhiKKpi) {
814  fMassHistPhi[indexMCKKpi]->Fill(invMass,weightKKpi);
815  if(fFillSparse) {
816  if(indexMCKKpi==GetSignalHistoIndex(iPtBin)) {
817  AliAODMCParticle *partDs = (AliAODMCParticle*)arrayMC->At(labDs);
818  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partDs,kTRUE);
819  if(orig==4) {
820  candType = 1.5;
821  }
822  if(orig==5) {
823  candType = 2.5;
824  }
825  }
826  }
827  }
828  if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass,weightKKpi);
829  }
830  }
831  if(ispiKK){
832  invMass=d->InvMassDspiKK();
833  fMassHist[index]->Fill(invMass,weightpiKK);
834  fPtVsMass->Fill(invMass,ptCand,weightpiKK);
835  if(isPhipiKK){
836  fMassHistPhi[index]->Fill(invMass,weightpiKK);
837  fPtVsMassPhi->Fill(invMass,ptCand,weightpiKK);
838  }
839  if(isK0starpiKK){
840  fMassHistK0st[index]->Fill(invMass,weightpiKK);
841  fPtVsMassK0st->Fill(invMass,ptCand,weightpiKK);
842  }
843  if(fReadMC && indexMCpiKK!=-1){
844  fMassHist[indexMCpiKK]->Fill(invMass,weightpiKK);
845  if(isPhipiKK) {
846  fMassHistPhi[indexMCpiKK]->Fill(invMass,weightpiKK);
847  if(fFillSparse) {
848  if(indexMCpiKK==GetSignalHistoIndex(iPtBin)) {
849  AliAODMCParticle *partDs = (AliAODMCParticle*)arrayMC->At(labDs);
850  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partDs,kTRUE);
851  if(orig==4) {
852  candType = 1.5;
853  }
854  if(orig==5) {
855  candType = 2.5;
856  }
857  }
858  }
859  }
860  if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass,weightpiKK);
861  }
862  }
863 
864 
866 
867  if(fFillSparse) {
868 
869  const Int_t nProng = 3;
870  Double_t deltaMassKK=999.;
871  Double_t dlen=d->DecayLength();
872  Double_t dlenxy=d->DecayLengthXY();
873  Double_t normdl=d->NormalizedDecayLength();
874  Double_t normdlxy=d->NormalizedDecayLengthXY();
875  Double_t cosp=d->CosPointingAngle();
876  Double_t cospxy=d->CosPointingAngleXY();
877  Double_t pt0=d->PtProng(0);
878  Double_t pt1=d->PtProng(1);
879  Double_t pt2=d->PtProng(2);
880  Double_t sigvert=d->GetSigmaVert();
881  Double_t cosPiDs=-99.;
882  Double_t cosPiKPhi=-99.;
883  Double_t normIP; //to store the maximum topomatic var. among the 3 prongs
884  Double_t normIPprong[nProng]; //to store IP of k,k,pi
885  if(isPhiKKpi) {
886  Double_t tmpNormIP[nProng];
887 
888  invMass = d->InvMassDsKKpi();
889  massKK = d->InvMass2Prongs(0,1,321,321);
890  deltaMassKK = TMath::Abs(massKK-massPhi);
891  cosPiDs = d->CosPiDsLabFrameKKpi();
892  cosPiKPhi = d->CosPiKPhiRFrameKKpi();
893  cosPiKPhi = TMath::Abs(cosPiKPhi*cosPiKPhi*cosPiKPhi);
894 
895  for(Int_t ip=0; ip<nProng; ip++) {
896  Double_t diffIP, errdiffIP;
897  d->Getd0MeasMinusExpProng(ip,aod->GetMagneticField(),diffIP,errdiffIP);
898  tmpNormIP[ip] = diffIP/errdiffIP;
899  if(ip==0) normIP = tmpNormIP[ip];
900  else if(TMath::Abs(tmpNormIP[ip])>TMath::Abs(normIP)) normIP = tmpNormIP[ip];
901  }
902  normIPprong[0] = tmpNormIP[0];
903  normIPprong[1] = tmpNormIP[1];
904  normIPprong[2] = tmpNormIP[2];
905 
906  Double_t var4nSparse[16] = {invMass,ptCand,deltaMassKK,dlen,dlenxy,normdl,normdlxy,cosp,cospxy,
907  sigvert,cosPiDs,cosPiKPhi,normIP,normIPprong[0],normIPprong[1],normIPprong[2]};
908 
909  if(!fReadMC) {
910  fnSparse->Fill(var4nSparse);
911  }
912  else {
913  if(indexMCKKpi==GetSignalHistoIndex(iPtBin)) {
914  if(candType==1.5) fnSparseMC[2]->Fill(var4nSparse);
915  if(candType==2.5) fnSparseMC[3]->Fill(var4nSparse);
916  }
917  }
918  }
919  if(isPhipiKK) {
920  Double_t tmpNormIP[nProng];
921 
922  invMass = d->InvMassDspiKK();
923  massKK = d->InvMass2Prongs(1,2,321,321);
924  deltaMassKK = TMath::Abs(massKK-massPhi);
925  cosPiDs = d->CosPiDsLabFramepiKK();
926  cosPiKPhi = d->CosPiKPhiRFramepiKK();
927  cosPiKPhi = TMath::Abs(cosPiKPhi*cosPiKPhi*cosPiKPhi);
928 
929  for(Int_t ip=0; ip<nProng; ip++) {
930  Double_t diffIP, errdiffIP;
931  d->Getd0MeasMinusExpProng(ip,aod->GetMagneticField(),diffIP,errdiffIP);
932  tmpNormIP[ip] = diffIP/errdiffIP;
933  if(ip==0) normIP = tmpNormIP[ip];
934  else if(TMath::Abs(tmpNormIP[ip])>TMath::Abs(normIP)) normIP = tmpNormIP[ip];
935  }
936 
937  normIPprong[0] = tmpNormIP[2];
938  normIPprong[1] = tmpNormIP[1];
939  normIPprong[2] = tmpNormIP[0];
940 
941  Double_t var4nSparse[16] = {invMass,ptCand,deltaMassKK,dlen,dlenxy,normdl,normdlxy,cosp,cospxy,
942  sigvert,cosPiDs,cosPiKPhi,normIP,normIPprong[0],normIPprong[1],normIPprong[2]};
943 
944  if(!fReadMC) {
945  fnSparse->Fill(var4nSparse);
946  }
947  else {
948  if(indexMCpiKK==GetSignalHistoIndex(iPtBin)) {
949  if(candType==1.5) fnSparseMC[2]->Fill(var4nSparse);
950  if(candType==2.5) fnSparseMC[3]->Fill(var4nSparse);
951  }
952  }
953  }
954 
955  if(fReadMC && (isPhiKKpi || isPhiKKpi)) {
956  Double_t var[6] = {ptCand,normIP,normIPprong[0],normIPprong[1],normIPprong[2],candType};
957  fnSparseIP->Fill(var);
958  }
959  }
961 
962  if(fDoCutVarHistos){
963  Double_t dlen=d->DecayLength();
964  Double_t cosp=d->CosPointingAngle();
965  Double_t pt0=d->PtProng(0);
966  Double_t pt1=d->PtProng(1);
967  Double_t pt2=d->PtProng(2);
968  Double_t sigvert=d->GetSigmaVert();
969  Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
970  Double_t dca=d->GetDCA();
971  Double_t ptmax=0;
972  for(Int_t i=0;i<3;i++){
973  if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
974  }
975  fCosPHist[index]->Fill(cosp);
976  fDLenHist[index]->Fill(dlen);
977  fSigVertHist[index]->Fill(sigvert);
978  fSumd02Hist[index]->Fill(sumD02);
979  fPtMaxHist[index]->Fill(ptmax);
980  fDCAHist[index]->Fill(dca);
981  fPtProng0Hist[index]->Fill(pt0);
982  fPtProng1Hist[index]->Fill(pt1);
983  fPtProng2Hist[index]->Fill(pt2);
984  if(isKKpi){
985  Double_t massKK=d->InvMass2Prongs(0,1,321,321);
986  Double_t massKp=d->InvMass2Prongs(1,2,321,211);
987  fDalitz[index]->Fill(massKK,massKp);
988  if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);
989  if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);
990  if(fReadMC && indexMCKKpi!=-1){
991  fDalitz[indexMCKKpi]->Fill(massKK,massKp);
992  if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);
993  if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);
994  fCosPHist[indexMCKKpi]->Fill(cosp);
995  fDLenHist[indexMCKKpi]->Fill(dlen);
996  fSigVertHist[indexMCKKpi]->Fill(sigvert);
997  fSumd02Hist[indexMCKKpi]->Fill(sumD02);
998  fPtMaxHist[indexMCKKpi]->Fill(ptmax);
999  fPtCandHist[indexMCKKpi]->Fill(ptCand);
1000  fDCAHist[indexMCKKpi]->Fill(dca);
1001  fPtProng0Hist[indexMCKKpi]->Fill(pt0);
1002  fPtProng1Hist[indexMCKKpi]->Fill(pt1);
1003  fPtProng2Hist[indexMCKKpi]->Fill(pt2);
1004  }
1005  }
1006  if(ispiKK){
1007  Double_t massKK=d->InvMass2Prongs(1,2,321,321);
1008  Double_t massKp=d->InvMass2Prongs(0,1,211,321);
1009  fDalitz[index]->Fill(massKK,massKp);
1010  if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);
1011  if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);
1012 
1013 
1014  if(fReadMC && indexMCpiKK!=-1){
1015  fDalitz[indexMCpiKK]->Fill(massKK,massKp);
1016  if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);
1017  if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);
1018  fCosPHist[indexMCpiKK]->Fill(cosp);
1019  fDLenHist[indexMCpiKK]->Fill(dlen);
1020  fSigVertHist[indexMCpiKK]->Fill(sigvert);
1021  fSumd02Hist[indexMCpiKK]->Fill(sumD02);
1022  fPtMaxHist[indexMCpiKK]->Fill(ptmax);
1023  fPtCandHist[indexMCpiKK]->Fill(ptCand);
1024  fDCAHist[indexMCpiKK]->Fill(dca);
1025  fPtProng0Hist[indexMCpiKK]->Fill(pt0);
1026  fPtProng1Hist[indexMCpiKK]->Fill(pt1);
1027  fPtProng2Hist[indexMCpiKK]->Fill(pt2);
1028  }
1029  }
1030 
1031 
1032  }
1033 
1034  Float_t tmp[37];
1035  if(fFillNtuple>0){
1036 
1037  if ((fFillNtuple==1 && (isPhiKKpi || isPhipiKK)) || (fFillNtuple==2 && (isK0starKKpi || isK0starpiKK)) || (fFillNtuple==3 && (isKKpi || ispiKK))){
1038 
1039  AliAODTrack *track0=(AliAODTrack*)d->GetDaughter(0);
1040  AliAODTrack *track1=(AliAODTrack*)d->GetDaughter(1);
1041  AliAODTrack *track2=(AliAODTrack*)d->GetDaughter(2);
1042 
1043  UInt_t bitMapPIDTrack0=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track0);
1044  UInt_t bitMapPIDTrack1=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track1);
1045  UInt_t bitMapPIDTrack2=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track2);
1046 
1047  tmp[0]=Float_t(labDs);
1048  if(fReadMC && fWriteOnlySignal) tmp[0]=Float_t(isMCSignal);
1049  tmp[1]=Float_t(retCodeAnalysisCuts);
1050  tmp[2]=Float_t(pdgCode0);
1051  tmp[3]=d->PtProng(0);
1052  tmp[4]=d->PtProng(1);
1053  tmp[5]=d->PtProng(2);
1054  tmp[6]=d->Pt();
1055  tmp[7]=d->PProng(0);
1056  tmp[8]=d->PProng(1);
1057  tmp[9]=d->PProng(2);
1058  tmp[10]=Int_t(bitMapPIDTrack0);
1059  tmp[11]=Int_t(bitMapPIDTrack1);
1060  tmp[12]=Int_t(bitMapPIDTrack2);
1061  tmp[13]=d->CosPointingAngle();
1062  tmp[14]=d->CosPointingAngleXY();
1063  tmp[15]=d->DecayLength();
1064  tmp[16]=d->DecayLengthXY();
1065  tmp[17]=d->NormalizedDecayLength();
1066  tmp[18]=d->NormalizedDecayLengthXY();
1067  tmp[19]=d->InvMassDsKKpi();
1068  tmp[20]=d->InvMassDspiKK();
1069  tmp[21]=d->GetSigmaVert();
1070  tmp[22]=d->Getd0Prong(0);
1071  tmp[23]=d->Getd0Prong(1);
1072  tmp[24]=d->Getd0Prong(2);
1073  tmp[25]=d->GetDCA();
1074  tmp[26]=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
1075  tmp[27]=d->InvMass2Prongs(0,1,321,321);
1076  tmp[28]=d->InvMass2Prongs(1,2,321,321);
1077  tmp[29]=d->InvMass2Prongs(1,2,321,211);
1078  tmp[30]=d->InvMass2Prongs(0,1,211,321);
1079  tmp[31]=d->CosPiDsLabFrameKKpi();
1080  tmp[32]=d->CosPiDsLabFramepiKK();
1081  tmp[33]=d->CosPiKPhiRFrameKKpi();
1082  tmp[34]=d->CosPiKPhiRFramepiKK();
1083  tmp[35]=(Float_t)(centrality);
1084  tmp[36]=(Float_t)(runNumber);
1085 
1086  if(fReadMC && fWriteOnlySignal){
1087  if(isMCSignal>=0) fNtupleDs->Fill(tmp);
1088  }else{
1089  fNtupleDs->Fill(tmp);
1090  }
1091  PostData(4,fNtupleDs);
1092  }
1093  }
1094 
1095  if(unsetvtx) d->UnsetOwnPrimaryVtx();
1096  if(recVtx)fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
1097  }
1098 
1099  fCounter->StoreCandidates(aod,nFiltered,kTRUE);
1100  fCounter->StoreCandidates(aod,nSelected,kFALSE);
1101 
1102  PostData(1,fOutput);
1103  PostData(3,fCounter);
1104 
1105  return;
1106 }
1107 
1108 //_________________________________________________________________
1109 
1110 void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
1111 {
1113  //
1114  if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
1115  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1116  if (!fOutput) {
1117  printf("ERROR: fOutput not available\n");
1118  return;
1119  }
1120  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
1121  if(fHistNEvents){
1122  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1123  }else{
1124  printf("ERROR: fHistNEvents not available\n");
1125  return;
1126  }
1127  return;
1128 }
1129 
1130 //_________________________________________________________________
1131 void AliAnalysisTaskSEDs::FillMCGenAccHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader){
1133 
1134  Int_t nProng = 3;
1135  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
1136  if(TMath::Abs(zMCVertex) <= fAnalysisCuts->GetMaxVtxZ()) {
1137  for(Int_t iPart=0; iPart<arrayMC->GetEntriesFast(); iPart++){
1138 
1139  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
1140 
1141  if(TMath::Abs(mcPart->GetPdgCode()) == 431) {
1142  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,kTRUE);//Prompt = 4, FeedDown = 5
1143 
1144  Int_t deca = 0;
1145  Bool_t isGoodDecay = kFALSE;
1146  Int_t labDau[3] = {-1,-1,-1};
1147  Bool_t isFidAcc = kFALSE;
1148  Bool_t isDaugInAcc = kFALSE;
1149 
1150  deca = AliVertexingHFUtils::CheckDsDecay(arrayMC,mcPart,labDau);
1151  if(deca == 1) isGoodDecay=kTRUE; // == 1 -> Phi pi -> kkpi
1152 
1153  if(labDau[0]==-1) continue; //protection against unfilled array of labels
1154 
1155  if(isGoodDecay) {
1156  Double_t pt = mcPart->Pt();
1157  Double_t rapid = mcPart->Y();
1158  isFidAcc = fAnalysisCuts->IsInFiducialAcceptance(pt,rapid);
1159  isDaugInAcc = CheckDaugAcc(arrayMC,nProng,labDau);
1160 
1161  if(isFidAcc) {
1162  Double_t var4nSparseAcc[2] = {pt,rapid};
1163  if(isDaugInAcc) {
1164  if(orig==4) fnSparseMC[0]->Fill(var4nSparseAcc);
1165  if(orig==5) fnSparseMC[1]->Fill(var4nSparseAcc);
1166  }
1167  }
1168  }
1169  }
1170  }
1171  }
1172 }
1173 //_________________________________________________________________
1174 Bool_t AliAnalysisTaskSEDs::CheckDaugAcc(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1176 
1177  for (Int_t iProng = 0; iProng<nProng; iProng++){
1178  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1179  if(!mcPartDaughter) {
1180  return kFALSE;
1181  }
1182  Double_t eta = mcPartDaughter->Eta();
1183  Double_t pt = mcPartDaughter->Pt();
1184  if (TMath::Abs(eta) > 0.9 || pt < 0.1) {
1185  return kFALSE;
1186  }
1187  }
1188  return kTRUE;
1189 }
1190 
1191 
1192 
1193 
1194 
Double_t NormalizedDecayLengthXY() const
Bool_t IsEventRejectedDueToCentrality() const
Definition: AliRDHFCuts.h:310
TH1F * fPtMaxHist[4 *kMaxPtBins]
! hist. for Pt Max (Prod Cuts)
Bool_t IsEventRejectedDueToZVertexOutsideFiducialRegion() const
Definition: AliRDHFCuts.h:304
Double_t NormalizedDecayLength() const
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
THnSparseF * fnSparseMC[4]
!<!THnSparse for topomatic variable
Bool_t IsEventRejectedDueToNotRecoVertex() const
Definition: AliRDHFCuts.h:298
void StoreCandidates(AliVEvent *, Int_t nCand=0, Bool_t flagFilter=kTRUE)
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
TH1F * fPtProng1Hist[4 *kMaxPtBins]
! hist. for DCA (Prod Cuts)
void SetPtBins(Int_t n, Float_t *lim)
Int_t IsEventSelectedInCentrality(AliVEvent *event)
Bool_t HasSelectionBit(Int_t i) const
TH2F * fDalitzK0st[4 *kMaxPtBins]
! dalitz plot via K0* (sig,bkg,tot)
Double_t CosPiDsLabFrameKKpi() const
centrality
Int_t GetHistoIndex(Int_t iPtBin) const
Double_t fMassRange
limits for pt bins
Bool_t IsEventRejectedDueToVertexContributors() const
Definition: AliRDHFCuts.h:301
TH2F * fPtVsMass
! hist. of pt vs. mass (prod. cuts)
static Int_t CheckDsDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
Double_t CosPointingAngleXY() const
TH1F * fPtProng0Hist[4 *kMaxPtBins]
! hist. for Pt Max (Prod Cuts)
TH2F * fYVsPtSig
! hist. of Y vs. Pt (MC, only sig, prod. cuts)
Int_t GetPidOption() const
virtual void UserCreateOutputObjects()
Implementation of interface methods.
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Bool_t fFillSparse
flag for usage of HasSelectionBit
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:241
TH1F * fSumd02Hist[4 *kMaxPtBins]
! hist. for sum d02 (Prod Cuts)
UInt_t GetPIDTrackTPCTOFBitMap(AliAODTrack *track) const
TH1F * fMassHistPhi[4 *kMaxPtBins]
! hist. of mass spectra via phi (sig,bkg,tot)
THnSparseF * fnSparse
Cuts for Analysis.
TNtuple * fNtupleDs
! output ntuple
TH1F * fHistCentrality[3]
!hist. for cent distr (all,sel ev, )
AliNormalizationCounter * fCounter
bin size for inv. mass histo
TH2F * fHistCentralityMult[3]
!hist. for cent distr vs mult (all,sel ev, )
const Double_t ptmax
TH1F * fHistNEvents
! hist. for No. of events
AliRDHFCutsDstoKKpi * fAnalysisCuts
TList * fOutput
! list send on output slot 0
Double_t CosPiKPhiRFrameKKpi() const
AliAODVertex * GetOwnPrimaryVtx() const
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Double_t GetSigmaVert(const AliAODEvent *aod=0x0)
THnSparseF * fnSparseIP
!<!THnSparse for candidates on data
Bool_t CheckDaugAcc(TClonesArray *arrayMC, Int_t nProng, Int_t *labDau)
TH2F * fDalitzPhi[4 *kMaxPtBins]
! dalitz plot via phi (sig,bkg,tot)
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:242
Bool_t IsEventRejectedDueToPileup() const
Definition: AliRDHFCuts.h:307
TH1F * fMassHistKpi[kMaxPtBins]
! hist. of mass spectra of Kpi
Bool_t IsCutOnResonancesApplied() const
TH1F * fPtCandHist[4 *kMaxPtBins]
! hist. for Pt Max (Prod Cuts)
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
void FillMCGenAccHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader)
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
Double_t CosPiDsLabFramepiKK() const
Double_t DecayLengthXY() const
Bool_t GetIsPrimaryWithoutDaughters() const
Definition: AliRDHFCuts.h:261
Int_t GetSignalHistoIndex(Int_t iPtBin) const
Bool_t IsEventSelected(AliVEvent *event)
virtual void UserExec(Option_t *option)
Int_t GetReflSignalHistoIndex(Int_t iPtBin) const
virtual void Terminate(Option_t *option)
TH1F * fCosPHist[4 *kMaxPtBins]
! hist. of cos pointing angle (sig,bkg,tot)
Double_t minMass
void CleanOwnPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod, AliAODVertex *origownvtx) const
Float_t * GetPtBinLimits() const
Definition: AliRDHFCuts.h:232
TH1F * fDLenHist[4 *kMaxPtBins]
! hist. of decay length (sig,bkg,tot)
Bool_t fUseSelectionBit
flag to create and fill histos with distributions of cut variables
Double_t GetWeightForpiKK() const
Double_t GetWeightForKKpi() const
TH1F * fChanHist[4]
! hist. with KKpi and piKK candidates (sig,bkg,tot)
TH1F * fMassHist[4 *kMaxPtBins]
! hist. of mass spectra (sig,bkg,tot)
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:233
TH1F * fDCAHist[4 *kMaxPtBins]
! hist. for DCA (Prod Cuts)
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:295
Bool_t fDoCutVarHistos
flag to control ntuple writing in MC
Double_t maxMass
TH1F * fPtProng2Hist[4 *kMaxPtBins]
! hist. for DCA (Prod Cuts)
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999)
Double_t CosPointingAngle() const
TH2F * fDalitz[4 *kMaxPtBins]
! dalitz plot (sig,bkg,tot)
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:263
TH1F * fSigVertHist[4 *kMaxPtBins]
! hist. for sigVert (Prod Cuts)
Bool_t RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod) const
Double_t DecayLength() const
TH2F * fPtVsMassPhi
! hist. of pt vs. mass (phi selection)
TH1F * fMassHistKK[kMaxPtBins]
! hist. of mass spectra of KK
Double_t CosPiKPhiRFramepiKK() const
Int_t nptbins
TH1F * fMassHistK0st[4 *kMaxPtBins]
! hist. of mass spectra via K0* (sig,bkg,tot)
TH2F * fYVsPt
! hist. of Y vs. Pt (prod. cuts)
TH2F * fPtVsMassK0st
! hist. of pt vs. mass (K0* selection)
TList * fListCuts
number of Pt bins
void ApplyCutOnResonances(Bool_t opt=kTRUE)
UChar_t fNPtBins
flag for usage of THnSparse
Bool_t fWriteOnlySignal
flag for access to MC
Double_t fMassBinSize
range for mass histogram
Int_t GetBackgroundHistoIndex(Int_t iPtBin) const
Float_t fPtLimits[kMaxPtBins+1]