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