AliPhysics  v5-07-15-01 (b3d7633)
 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 ",13,-0.5,12.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  fHistNEvents->GetXaxis()->SetBinLabel(12,"no. of not on-the-fly rec Ds");
324 
325  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
326 
327  fHistNEvents->Sumw2();
328  fHistNEvents->SetMinimum(0);
329  fOutput->Add(fHistNEvents);
330 
331  fHistCentrality[0]=new TH1F("hCentr","centrality",10000,0.,100.);
332  fHistCentrality[1]=new TH1F("hCentr(selectedCent)","centrality(selectedCent)",10000,0.,100.);
333  fHistCentrality[2]=new TH1F("hCentr(OutofCent)","centrality(OutofCent)",10000,0.,100.);
334  fHistCentralityMult[0]=new TH2F("hCentrMult","centrality vs mult",100,0.5,30000.5,40,0.,100.);
335  fHistCentralityMult[1]=new TH2F("hCentrMult(selectedCent)","centrality vs mult(selectedCent)",100,0.5,30000.5,40,0.,100.);
336  fHistCentralityMult[2]=new TH2F("hCentrMult(OutofCent)","centrality vs mult(OutofCent)",100,0.5,30000.5,40,0.,100.);
337  for(Int_t i=0;i<3;i++){
338  fHistCentrality[i]->Sumw2();
339  fOutput->Add(fHistCentrality[i]);
340  fHistCentralityMult[i]->Sumw2();
341  fOutput->Add(fHistCentralityMult[i]);
342  }
343 
344  Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
345 
346  Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);
347  if(nInvMassBins%2==1) nInvMassBins++;
348  // Double_t minMass=massDs-1.0*nInvMassBins*fMassBinSize;
349  Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
350  // Double_t maxMass=massDs+1.0*nInvMassBins*fMassBinSize;
351  Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
352 
353  TString hisname;
354  TString htype;
355  Int_t index;
356  for(Int_t iType=0; iType<4; iType++){
357  for(Int_t i=0;i<fNPtBins;i++){
358  if(iType==0){
359  htype="All";
360  index=GetHistoIndex(i);
361  }else if(iType==1){
362  htype="Sig";
363  index=GetSignalHistoIndex(i);
364  }else if(iType==2){
365  htype="Bkg";
366  index=GetBackgroundHistoIndex(i);
367  }else{
368  htype="ReflSig";
369  index=GetReflSignalHistoIndex(i);
370  }
371  hisname.Form("hMass%sPt%d",htype.Data(),i);
372  fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
373  fMassHist[index]->Sumw2();
374  hisname.Form("hMass%sPt%dphi",htype.Data(),i);
375  fMassHistPhi[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
376  fMassHistPhi[index]->Sumw2();
377  hisname.Form("hMass%sPt%dk0st",htype.Data(),i);
378  fMassHistK0st[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
379  fMassHistK0st[index]->Sumw2();
380  hisname.Form("hCosP%sPt%d",htype.Data(),i);
381  fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
382  fCosPHist[index]->Sumw2();
383  hisname.Form("hDLen%sPt%d",htype.Data(),i);
384  fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
385  fDLenHist[index]->Sumw2();
386  hisname.Form("hSumd02%sPt%d",htype.Data(),i);
387  fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
388  fSumd02Hist[index]->Sumw2();
389  hisname.Form("hSigVert%sPt%d",htype.Data(),i);
390  fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
391  fSigVertHist[index]->Sumw2();
392  hisname.Form("hPtMax%sPt%d",htype.Data(),i);
393  fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
394  fPtMaxHist[index]->Sumw2();
395  hisname.Form("hPtCand%sPt%d",htype.Data(),i);
396  fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
397  fPtCandHist[index]->Sumw2();
398  hisname.Form("hDCA%sPt%d",htype.Data(),i);
399  fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
400  fDCAHist[index]->Sumw2();
401  hisname.Form("hPtProng0%sPt%d",htype.Data(),i);
402  fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
403  fPtProng0Hist[index]->Sumw2();
404  hisname.Form("hPtProng1%sPt%d",htype.Data(),i);
405  fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
406  fPtProng1Hist[index]->Sumw2();
407  hisname.Form("hPtProng2%sPt%d",htype.Data(),i);
408  fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
409  fPtProng2Hist[index]->Sumw2();
410  hisname.Form("hDalitz%sPt%d",htype.Data(),i);
411  fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
412  fDalitz[index]->Sumw2();
413  hisname.Form("hDalitz%sPt%dphi",htype.Data(),i);
414  fDalitzPhi[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
415  fDalitzPhi[index]->Sumw2();
416  hisname.Form("hDalitz%sPt%dk0st",htype.Data(),i);
417  fDalitzK0st[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
418  fDalitzK0st[index]->Sumw2();
419  }
420  }
421 
422  for(Int_t i=0; i<4*fNPtBins; i++){
423  fOutput->Add(fMassHist[i]);
424  fOutput->Add(fMassHistPhi[i]);
425  fOutput->Add(fMassHistK0st[i]);
426  fOutput->Add(fPtCandHist[i]);
427  if(fDoCutVarHistos){
428  fOutput->Add(fCosPHist[i]);
429  fOutput->Add(fDLenHist[i]);
430  fOutput->Add(fSumd02Hist[i]);
431  fOutput->Add(fSigVertHist[i]);
432  fOutput->Add(fPtMaxHist[i]);
433  fOutput->Add(fDCAHist[i]);
434  fOutput->Add(fPtProng0Hist[i]);
435  fOutput->Add(fPtProng1Hist[i]);
436  fOutput->Add(fPtProng2Hist[i]);
437  fOutput->Add(fDalitz[i]);
438  fOutput->Add(fDalitzPhi[i]);
439  fOutput->Add(fDalitzK0st[i]);
440  }
441  }
442 
443  fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",64,-0.5,63.5);
444  fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",64,-0.5,63.5);
445  fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",64,-0.5,63.5);
446  fChanHist[3] = new TH1F("hChanReflSig", "KKpi and piKK candidates",64,-0.5,63.5);
447  for(Int_t i=0;i<4;i++){
448  fChanHist[i]->Sumw2();
449  fChanHist[i]->SetMinimum(0);
450  fOutput->Add(fChanHist[i]);
451  }
452 
453 
454  fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nInvMassBins,minMass,maxMass,40,0.,20.);
455  fPtVsMassPhi=new TH2F("hPtVsMassPhi","PtVsMass (phi selection)",nInvMassBins,minMass,maxMass,200,0.,20.);
456  fPtVsMassK0st=new TH2F("hPtVsMassK0st","PtVsMass (K0* selection)",nInvMassBins,minMass,maxMass,200,0.,20.);
457  fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
458  fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
459 
460  for(Int_t i=0;i<fNPtBins;i++){
461  hisname.Form("hMassKKPt%d",i);
462  fMassHistKK[i]=new TH1F(hisname.Data(),hisname.Data(),200,0.95,1.35);
463  fMassHistKK[i]->Sumw2();
464  fOutput->Add(fMassHistKK[i]);
465  hisname.Form("hMassKpiPt%d",i);
466  fMassHistKpi[i]=new TH1F(hisname.Data(),hisname.Data(),200,0.7,1.1);
467  fMassHistKpi[i]->Sumw2();
468  fOutput->Add(fMassHistKpi[i]);
469  }
470 
471  fOutput->Add(fPtVsMass);
472  fOutput->Add(fPtVsMassPhi);
473  fOutput->Add(fPtVsMassK0st);
474  fOutput->Add(fYVsPt);
475  fOutput->Add(fYVsPtSig);
476 
477  const Int_t nvarsReco = 16;
478  const Int_t nvarsAcc = 2;
479  Int_t nBinsReco[nvarsReco] = {500,20,60,1000,1000,500,500,50,50,100,100,100,400,400,400,400};
480  Int_t nBinsAcc[nvarsAcc] = {100,20};
481  Double_t xminAcc[nvarsAcc] = {0.,-1.};
482  Double_t xmaxAcc[nvarsAcc] = {100.,1.};
483  Double_t xminReco[nvarsReco] = {1.5,0.,0.,0.,0.,0.,0.,0.5,0.5,0.,0.,0.,-10.,-10.,-10.,-10.};
484  Double_t xmaxReco[nvarsReco] = {2.5,20.,0.03,5.,5.,500.,500.,1.,1.,0.1,1.,1.,10.,10.,10.,10.};
485  TString axis[nvarsReco] = {"invMassAllPhi","p_{T}","#Delta Mass(KK)","dlen","dlen_{xy}","normdl","normdl_{xy}","cosP","cosP_{xy}",
486  "sigVert","cosPiDs","|cosPiKPhi^{3}|","normIP","normIP_kaon1","normIP_kaon2","normIP_pion"};
487 
488  const Int_t nvarsIP = 6;
489  Int_t nBinsIP[nvarsIP] = {20,400,400,400,400,3};
490  Double_t xminIP[nvarsIP] = {0.,-10.,-10.,-10.,-10.,0.};
491  Double_t xmaxIP[nvarsIP] = {20.,10.,10.,10.,10.,3.};
492  TString axisIP[nvarsIP] = {"motherPt","maxNormImp","IP0","IP1","IP2","candType"};
493 
494  if(fFillSparse) {
495 
496  if(fReadMC) {
497  TString label[nvarsAcc] = {"fromC","fromB"};
498  for (Int_t i=0; i<2; i++) {
499  fnSparseMC[i] = new THnSparseF(Form("fnSparseAcc_%s",label[i].Data()),Form("MC nSparse (Acc.Step)- %s",label[i].Data()),
500  nvarsAcc, nBinsAcc, xminAcc, xmaxAcc);
501  fnSparseMC[i]->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
502  fnSparseMC[i]->GetAxis(1)->SetTitle("y");
503  fOutput->Add(fnSparseMC[i]);
504  }
505  for (Int_t i=2; i<4; i++) {
506  fnSparseMC[i] = new THnSparseF(Form("fnSparseReco_%s",label[i-2].Data()),Form("MC nSparse (Reco Step)- %s",label[i-2].Data()),
507  nvarsReco, nBinsReco, xminReco, xmaxReco);
508  for (Int_t j=0; j<nvarsReco; j++) {
509  fnSparseMC[i]->GetAxis(j)->SetTitle(Form("%s",axis[j].Data()));
510  }
511  fOutput->Add(fnSparseMC[i]);
512  }
513 
514  fnSparseIP = new THnSparseF("fnSparseIP","nSparseIP", nvarsIP, nBinsIP, xminIP, xmaxIP);
515  for (Int_t j=0; j<nvarsIP; j++) {
516  fnSparseIP->GetAxis(j)->SetTitle(Form("%s",axisIP[j].Data()));
517  }
518  fnSparseIP->GetAxis(5)->SetTitle("candType (0.5=bkg; 1.5=prompt; 2.5=FD)");
519  fOutput->Add(fnSparseIP);
520  }
521  else {
522  fnSparse = new THnSparseF("fnSparse","nSparse", nvarsReco, nBinsReco, xminReco, xmaxReco);
523  for (Int_t j=0; j<nvarsReco; j++) {
524  fnSparse->GetAxis(j)->SetTitle(Form("%s",axis[j].Data()));
525  }
526  fOutput->Add(fnSparse);
527  }
528  }
529 
530  //Counter for Normalization
531  fCounter = new AliNormalizationCounter("NormalizationCounter");
532  fCounter->Init();
533 
534  PostData(1,fOutput);
535  PostData(3,fCounter);
536 
537  if(fFillNtuple>0){
538  OpenFile(4); // 4 is the slot number of the ntuple
539 
540  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");
541 
542  }
543 
544 
545  return;
546 }
547 
548 //________________________________________________________________________
549 void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
550 {
553 
554  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
555 
556  TClonesArray *array3Prong = 0;
557  if(!aod && AODEvent() && IsStandardAOD()) {
558  // In case there is an AOD handler writing a standard AOD, use the AOD
559  // event in memory rather than the input (ESD) event.
560  aod = dynamic_cast<AliAODEvent*> (AODEvent());
561  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
562  // have to taken from the AOD event hold by the AliAODExtension
563  AliAODHandler* aodHandler = (AliAODHandler*)
564  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
565  if(aodHandler->GetExtensions()) {
566  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
567  AliAODEvent *aodFromExt = ext->GetAOD();
568  array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
569  }
570  } else if(aod) {
571  array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
572  }
573 
574  if(!aod || !array3Prong) {
575  printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
576  return;
577  }
578 
579 
580  // fix for temporary bug in ESDfilter
581  // the AODs with null vertex pointer didn't pass the PhysSel
582  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
583 
584 
585  fHistNEvents->Fill(0); // count event
586  // Post the data already here
587  PostData(1,fOutput);
588 
590 
591 
592  Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
593  Float_t ntracks=aod->GetNumberOfTracks();
594  Float_t evCentr=fAnalysisCuts->GetCentrality(aod);
595 
596  fHistCentrality[0]->Fill(evCentr);
597  fHistCentralityMult[0]->Fill(ntracks,evCentr);
604  fHistNEvents->Fill(7);
605  fHistCentrality[2]->Fill(evCentr);
606  fHistCentralityMult[2]->Fill(ntracks,evCentr);
607  }
608 
609  Float_t centrality=fAnalysisCuts->GetCentrality(aod);
610  Int_t runNumber=aod->GetRunNumber();
611 
612 
613 
614  TClonesArray *arrayMC=0;
615  AliAODMCHeader *mcHeader=0;
616 
617  // AOD primary vertex
618  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
619  // vtx1->Print();
620 
621  // load MC particles
622  if(fReadMC){
623 
624  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
625  if(!arrayMC) {
626  printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
627  return;
628  }
629 
630  // load MC header
631  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
632  if(!mcHeader) {
633  printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
634  return;
635  }
636  }
637 
638 
639  if(fReadMC && fFillSparse){
640  if(aod->GetTriggerMask()==0 && (runNumber>=195344 && runNumber<=195677))
641  // protection for events with empty trigger mask in p-Pb
642  return;
644  // events not passing the centrality selection can be removed immediately.
645  return;
646  Double_t zMCVertex = mcHeader->GetVtxZ();
647  if (TMath::Abs(zMCVertex) > fAnalysisCuts->GetMaxVtxZ())
648  return;
649  FillMCGenAccHistos(arrayMC, mcHeader);
650  }
651 
652  if(!isEvSel)return;
653  fHistNEvents->Fill(1);
654  fHistCentrality[1]->Fill(evCentr);
655  fHistCentralityMult[1]->Fill(ntracks,evCentr);
656 
657  Int_t n3Prong = array3Prong->GetEntriesFast();
658  if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
659 
660 
661  Int_t pdgDstoKKpi[3]={321,321,211};
662  Int_t nSelected=0;
663  Int_t nFiltered=0;
664  Double_t massPhi=TDatabasePDG::Instance()->GetParticle(333)->Mass();
665 
666  // vHF object is needed to call the method that refills the missing info of the candidates
667  // if they have been deleted in dAOD reconstruction phase
668  // in order to reduce the size of the file
670 
671  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
672 
673  AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
674  fHistNEvents->Fill(8);
675 
677  continue;
678  }
679  nFiltered++;
680  fHistNEvents->Fill(9);
681 
682  if(!(vHF->FillRecoCand(aod,d))) {
683  fHistNEvents->Fill(18); //monitor how often this fails
684  continue;
685  }
686 
687  Bool_t unsetvtx=kFALSE;
688  if(!d->GetOwnPrimaryVtx()){
689  d->SetOwnPrimaryVtx(vtx1);
690  unsetvtx=kTRUE;
691  }
692 
693  Bool_t recVtx=kFALSE;
694  AliAODVertex *origownvtx=0x0;
695 
696  Double_t ptCand = d->Pt();
697  Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
698  Double_t rapid=d->YDs();
699  fYVsPt->Fill(ptCand,rapid);
700  Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);
701 
702  if(!isFidAcc) continue;
703 
704  Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
705  Int_t retCodeNoRes=retCodeAnalysisCuts;
706  Bool_t origRes=fAnalysisCuts->IsCutOnResonancesApplied();
707  if(origRes){
709  retCodeNoRes=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
711  }
712  Double_t massKK = 0.;
713 
714  if(retCodeNoRes&1){ //KKpi
715  massKK=d->InvMass2Prongs(0,1,321,321);
716  Double_t massKp=d->InvMass2Prongs(1,2,321,211);
717  fMassHistKK[iPtBin]->Fill(massKK);
718  fMassHistKpi[iPtBin]->Fill(massKp);
719  }
720  if(retCodeNoRes&2){ //piKK
721  massKK=d->InvMass2Prongs(1,2,321,321);
722  Double_t massKp=d->InvMass2Prongs(0,1,211,321);
723  fMassHistKK[iPtBin]->Fill(massKK);
724  fMassHistKpi[iPtBin]->Fill(massKp);
725  }
726 
727  Int_t isKKpi=retCodeAnalysisCuts&1;
728  Int_t ispiKK=retCodeAnalysisCuts&2;
729  Int_t isPhiKKpi=retCodeAnalysisCuts&4;
730  Int_t isPhipiKK=retCodeAnalysisCuts&8;
731  Int_t isK0starKKpi=retCodeAnalysisCuts&16;
732  Int_t isK0starpiKK=retCodeAnalysisCuts&32;
733 
734  if(retCodeAnalysisCuts<=0) continue;
736  if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
737  if(fAnalysisCuts->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
738  else fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
739  }
740 
741 
742  fHistNEvents->Fill(10);
743  nSelected++;
744 
745  Int_t index=GetHistoIndex(iPtBin);
746  fPtCandHist[index]->Fill(ptCand);
747 
748 
749  Double_t weightKKpi=1.;
750  Double_t weightpiKK=1.;
752  weightKKpi=fAnalysisCuts->GetWeightForKKpi();
753  weightpiKK=fAnalysisCuts->GetWeightForpiKK();
754  if(weightKKpi>1. || weightKKpi<0.) weightKKpi=0.;
755  if(weightpiKK>1. || weightpiKK<0.) weightpiKK=0.;
756  }
757 
758  fChanHist[0]->Fill(retCodeAnalysisCuts);
759 
760 
761  Double_t invMass = 0.;
762  Int_t indexMCKKpi=-1;
763  Int_t indexMCpiKK=-1;
764  Int_t labDs=-1;
765  Int_t pdgCode0=-999;
766  Int_t isMCSignal=-1;
767 
768 
769  if(fReadMC){
770 
771  labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
772  if(labDs>=0){
773  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
774  AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
775  pdgCode0=TMath::Abs(p->GetPdgCode());
776 
777  if(isKKpi){
778  if(pdgCode0==321) {
779  indexMCKKpi=GetSignalHistoIndex(iPtBin);
780  fYVsPtSig->Fill(ptCand,rapid);
781  fChanHist[1]->Fill(retCodeAnalysisCuts);
782  isMCSignal=1;
783  }else{
784  indexMCKKpi=GetReflSignalHistoIndex(iPtBin);
785  fChanHist[3]->Fill(retCodeAnalysisCuts);
786  isMCSignal=0;
787  }
788  }
789  if(ispiKK){
790  if(pdgCode0==211) {
791  indexMCpiKK=GetSignalHistoIndex(iPtBin);
792  fYVsPtSig->Fill(ptCand,rapid);
793  fChanHist[1]->Fill(retCodeAnalysisCuts);
794  isMCSignal=1;
795  }else{
796  indexMCpiKK=GetReflSignalHistoIndex(iPtBin);
797  fChanHist[3]->Fill(retCodeAnalysisCuts);
798  isMCSignal=0;
799  }
800  }
801  }else{
802  indexMCpiKK=GetBackgroundHistoIndex(iPtBin);
803  indexMCKKpi=GetBackgroundHistoIndex(iPtBin);
804  fChanHist[2]->Fill(retCodeAnalysisCuts);
805  }
806  }
807 
808  Double_t candType = 0.5; //for bkg
809 
810  if(isKKpi){
811  invMass=d->InvMassDsKKpi();
812  fMassHist[index]->Fill(invMass,weightKKpi);
813  fPtVsMass->Fill(invMass,ptCand,weightKKpi);
814  if(isPhiKKpi){
815  fMassHistPhi[index]->Fill(invMass,weightKKpi);
816  fPtVsMassPhi->Fill(invMass,ptCand,weightKKpi);
817  }
818  if(isK0starKKpi){
819  fMassHistK0st[index]->Fill(invMass,weightKKpi);
820  fPtVsMassK0st->Fill(invMass,ptCand,weightKKpi);
821  }
822  if(fReadMC && indexMCKKpi!=-1){
823  fMassHist[indexMCKKpi]->Fill(invMass,weightKKpi);
824  if(isPhiKKpi) {
825  fMassHistPhi[indexMCKKpi]->Fill(invMass,weightKKpi);
826  if(fFillSparse) {
827  if(indexMCKKpi==GetSignalHistoIndex(iPtBin)) {
828  AliAODMCParticle *partDs = (AliAODMCParticle*)arrayMC->At(labDs);
829  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partDs,kTRUE);
830  if(orig==4) {
831  candType = 1.5;
832  }
833  if(orig==5) {
834  candType = 2.5;
835  }
836  }
837  }
838  }
839  if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass,weightKKpi);
840  }
841  }
842  if(ispiKK){
843  invMass=d->InvMassDspiKK();
844  fMassHist[index]->Fill(invMass,weightpiKK);
845  fPtVsMass->Fill(invMass,ptCand,weightpiKK);
846  if(isPhipiKK){
847  fMassHistPhi[index]->Fill(invMass,weightpiKK);
848  fPtVsMassPhi->Fill(invMass,ptCand,weightpiKK);
849  }
850  if(isK0starpiKK){
851  fMassHistK0st[index]->Fill(invMass,weightpiKK);
852  fPtVsMassK0st->Fill(invMass,ptCand,weightpiKK);
853  }
854  if(fReadMC && indexMCpiKK!=-1){
855  fMassHist[indexMCpiKK]->Fill(invMass,weightpiKK);
856  if(isPhipiKK) {
857  fMassHistPhi[indexMCpiKK]->Fill(invMass,weightpiKK);
858  if(fFillSparse) {
859  if(indexMCpiKK==GetSignalHistoIndex(iPtBin)) {
860  AliAODMCParticle *partDs = (AliAODMCParticle*)arrayMC->At(labDs);
861  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partDs,kTRUE);
862  if(orig==4) {
863  candType = 1.5;
864  }
865  if(orig==5) {
866  candType = 2.5;
867  }
868  }
869  }
870  }
871  if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass,weightpiKK);
872  }
873  }
874 
875 
877 
878  if(fFillSparse) {
879 
880  const Int_t nProng = 3;
881  Double_t deltaMassKK=999.;
882  Double_t dlen=d->DecayLength();
883  Double_t dlenxy=d->DecayLengthXY();
884  Double_t normdl=d->NormalizedDecayLength();
885  Double_t normdlxy=d->NormalizedDecayLengthXY();
886  Double_t cosp=d->CosPointingAngle();
887  Double_t cospxy=d->CosPointingAngleXY();
888  Double_t pt0=d->PtProng(0);
889  Double_t pt1=d->PtProng(1);
890  Double_t pt2=d->PtProng(2);
891  Double_t sigvert=d->GetSigmaVert();
892  Double_t cosPiDs=-99.;
893  Double_t cosPiKPhi=-99.;
894  Double_t normIP; //to store the maximum topomatic var. among the 3 prongs
895  Double_t normIPprong[nProng]; //to store IP of k,k,pi
896  if(isPhiKKpi) {
897  Double_t tmpNormIP[nProng];
898 
899  invMass = d->InvMassDsKKpi();
900  massKK = d->InvMass2Prongs(0,1,321,321);
901  deltaMassKK = TMath::Abs(massKK-massPhi);
902  cosPiDs = d->CosPiDsLabFrameKKpi();
903  cosPiKPhi = d->CosPiKPhiRFrameKKpi();
904  cosPiKPhi = TMath::Abs(cosPiKPhi*cosPiKPhi*cosPiKPhi);
905 
906  for(Int_t ip=0; ip<nProng; ip++) {
907  Double_t diffIP, errdiffIP;
908  d->Getd0MeasMinusExpProng(ip,aod->GetMagneticField(),diffIP,errdiffIP);
909  tmpNormIP[ip] = diffIP/errdiffIP;
910  if(ip==0) normIP = tmpNormIP[ip];
911  else if(TMath::Abs(tmpNormIP[ip])>TMath::Abs(normIP)) normIP = tmpNormIP[ip];
912  }
913  normIPprong[0] = tmpNormIP[0];
914  normIPprong[1] = tmpNormIP[1];
915  normIPprong[2] = tmpNormIP[2];
916 
917  Double_t var4nSparse[16] = {invMass,ptCand,deltaMassKK,dlen,dlenxy,normdl,normdlxy,cosp,cospxy,
918  sigvert,cosPiDs,cosPiKPhi,normIP,normIPprong[0],normIPprong[1],normIPprong[2]};
919 
920  if(!fReadMC) {
921  fnSparse->Fill(var4nSparse);
922  }
923  else {
924  if(indexMCKKpi==GetSignalHistoIndex(iPtBin)) {
925  if(candType==1.5) fnSparseMC[2]->Fill(var4nSparse);
926  if(candType==2.5) fnSparseMC[3]->Fill(var4nSparse);
927  }
928  }
929  }
930  if(isPhipiKK) {
931  Double_t tmpNormIP[nProng];
932 
933  invMass = d->InvMassDspiKK();
934  massKK = d->InvMass2Prongs(1,2,321,321);
935  deltaMassKK = TMath::Abs(massKK-massPhi);
936  cosPiDs = d->CosPiDsLabFramepiKK();
937  cosPiKPhi = d->CosPiKPhiRFramepiKK();
938  cosPiKPhi = TMath::Abs(cosPiKPhi*cosPiKPhi*cosPiKPhi);
939 
940  for(Int_t ip=0; ip<nProng; ip++) {
941  Double_t diffIP, errdiffIP;
942  d->Getd0MeasMinusExpProng(ip,aod->GetMagneticField(),diffIP,errdiffIP);
943  tmpNormIP[ip] = diffIP/errdiffIP;
944  if(ip==0) normIP = tmpNormIP[ip];
945  else if(TMath::Abs(tmpNormIP[ip])>TMath::Abs(normIP)) normIP = tmpNormIP[ip];
946  }
947 
948  normIPprong[0] = tmpNormIP[2];
949  normIPprong[1] = tmpNormIP[1];
950  normIPprong[2] = tmpNormIP[0];
951 
952  Double_t var4nSparse[16] = {invMass,ptCand,deltaMassKK,dlen,dlenxy,normdl,normdlxy,cosp,cospxy,
953  sigvert,cosPiDs,cosPiKPhi,normIP,normIPprong[0],normIPprong[1],normIPprong[2]};
954 
955  if(!fReadMC) {
956  fnSparse->Fill(var4nSparse);
957  }
958  else {
959  if(indexMCpiKK==GetSignalHistoIndex(iPtBin)) {
960  if(candType==1.5) fnSparseMC[2]->Fill(var4nSparse);
961  if(candType==2.5) fnSparseMC[3]->Fill(var4nSparse);
962  }
963  }
964  }
965 
966  if(fReadMC && (isPhiKKpi || isPhiKKpi)) {
967  Double_t var[6] = {ptCand,normIP,normIPprong[0],normIPprong[1],normIPprong[2],candType};
968  fnSparseIP->Fill(var);
969  }
970  }
972 
973  if(fDoCutVarHistos){
974  Double_t dlen=d->DecayLength();
975  Double_t cosp=d->CosPointingAngle();
976  Double_t pt0=d->PtProng(0);
977  Double_t pt1=d->PtProng(1);
978  Double_t pt2=d->PtProng(2);
979  Double_t sigvert=d->GetSigmaVert();
980  Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
981  Double_t dca=d->GetDCA();
982  Double_t ptmax=0;
983  for(Int_t i=0;i<3;i++){
984  if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
985  }
986  fCosPHist[index]->Fill(cosp);
987  fDLenHist[index]->Fill(dlen);
988  fSigVertHist[index]->Fill(sigvert);
989  fSumd02Hist[index]->Fill(sumD02);
990  fPtMaxHist[index]->Fill(ptmax);
991  fDCAHist[index]->Fill(dca);
992  fPtProng0Hist[index]->Fill(pt0);
993  fPtProng1Hist[index]->Fill(pt1);
994  fPtProng2Hist[index]->Fill(pt2);
995  if(isKKpi){
996  Double_t massKK=d->InvMass2Prongs(0,1,321,321);
997  Double_t massKp=d->InvMass2Prongs(1,2,321,211);
998  fDalitz[index]->Fill(massKK,massKp);
999  if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);
1000  if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);
1001  if(fReadMC && indexMCKKpi!=-1){
1002  fDalitz[indexMCKKpi]->Fill(massKK,massKp);
1003  if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);
1004  if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);
1005  fCosPHist[indexMCKKpi]->Fill(cosp);
1006  fDLenHist[indexMCKKpi]->Fill(dlen);
1007  fSigVertHist[indexMCKKpi]->Fill(sigvert);
1008  fSumd02Hist[indexMCKKpi]->Fill(sumD02);
1009  fPtMaxHist[indexMCKKpi]->Fill(ptmax);
1010  fPtCandHist[indexMCKKpi]->Fill(ptCand);
1011  fDCAHist[indexMCKKpi]->Fill(dca);
1012  fPtProng0Hist[indexMCKKpi]->Fill(pt0);
1013  fPtProng1Hist[indexMCKKpi]->Fill(pt1);
1014  fPtProng2Hist[indexMCKKpi]->Fill(pt2);
1015  }
1016  }
1017  if(ispiKK){
1018  Double_t massKK=d->InvMass2Prongs(1,2,321,321);
1019  Double_t massKp=d->InvMass2Prongs(0,1,211,321);
1020  fDalitz[index]->Fill(massKK,massKp);
1021  if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);
1022  if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);
1023 
1024 
1025  if(fReadMC && indexMCpiKK!=-1){
1026  fDalitz[indexMCpiKK]->Fill(massKK,massKp);
1027  if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);
1028  if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);
1029  fCosPHist[indexMCpiKK]->Fill(cosp);
1030  fDLenHist[indexMCpiKK]->Fill(dlen);
1031  fSigVertHist[indexMCpiKK]->Fill(sigvert);
1032  fSumd02Hist[indexMCpiKK]->Fill(sumD02);
1033  fPtMaxHist[indexMCpiKK]->Fill(ptmax);
1034  fPtCandHist[indexMCpiKK]->Fill(ptCand);
1035  fDCAHist[indexMCpiKK]->Fill(dca);
1036  fPtProng0Hist[indexMCpiKK]->Fill(pt0);
1037  fPtProng1Hist[indexMCpiKK]->Fill(pt1);
1038  fPtProng2Hist[indexMCpiKK]->Fill(pt2);
1039  }
1040  }
1041 
1042 
1043  }
1044 
1045  Float_t tmp[37];
1046  if(fFillNtuple>0){
1047 
1048  if ((fFillNtuple==1 && (isPhiKKpi || isPhipiKK)) || (fFillNtuple==2 && (isK0starKKpi || isK0starpiKK)) || (fFillNtuple==3 && (isKKpi || ispiKK))){
1049 
1050  AliAODTrack *track0=(AliAODTrack*)d->GetDaughter(0);
1051  AliAODTrack *track1=(AliAODTrack*)d->GetDaughter(1);
1052  AliAODTrack *track2=(AliAODTrack*)d->GetDaughter(2);
1053 
1054  UInt_t bitMapPIDTrack0=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track0);
1055  UInt_t bitMapPIDTrack1=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track1);
1056  UInt_t bitMapPIDTrack2=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track2);
1057 
1058  tmp[0]=Float_t(labDs);
1059  if(fReadMC && fWriteOnlySignal) tmp[0]=Float_t(isMCSignal);
1060  tmp[1]=Float_t(retCodeAnalysisCuts);
1061  tmp[2]=Float_t(pdgCode0);
1062  tmp[3]=d->PtProng(0);
1063  tmp[4]=d->PtProng(1);
1064  tmp[5]=d->PtProng(2);
1065  tmp[6]=d->Pt();
1066  tmp[7]=d->PProng(0);
1067  tmp[8]=d->PProng(1);
1068  tmp[9]=d->PProng(2);
1069  tmp[10]=Int_t(bitMapPIDTrack0);
1070  tmp[11]=Int_t(bitMapPIDTrack1);
1071  tmp[12]=Int_t(bitMapPIDTrack2);
1072  tmp[13]=d->CosPointingAngle();
1073  tmp[14]=d->CosPointingAngleXY();
1074  tmp[15]=d->DecayLength();
1075  tmp[16]=d->DecayLengthXY();
1076  tmp[17]=d->NormalizedDecayLength();
1077  tmp[18]=d->NormalizedDecayLengthXY();
1078  tmp[19]=d->InvMassDsKKpi();
1079  tmp[20]=d->InvMassDspiKK();
1080  tmp[21]=d->GetSigmaVert();
1081  tmp[22]=d->Getd0Prong(0);
1082  tmp[23]=d->Getd0Prong(1);
1083  tmp[24]=d->Getd0Prong(2);
1084  tmp[25]=d->GetDCA();
1085  tmp[26]=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
1086  tmp[27]=d->InvMass2Prongs(0,1,321,321);
1087  tmp[28]=d->InvMass2Prongs(1,2,321,321);
1088  tmp[29]=d->InvMass2Prongs(1,2,321,211);
1089  tmp[30]=d->InvMass2Prongs(0,1,211,321);
1090  tmp[31]=d->CosPiDsLabFrameKKpi();
1091  tmp[32]=d->CosPiDsLabFramepiKK();
1092  tmp[33]=d->CosPiKPhiRFrameKKpi();
1093  tmp[34]=d->CosPiKPhiRFramepiKK();
1094  tmp[35]=(Float_t)(centrality);
1095  tmp[36]=(Float_t)(runNumber);
1096 
1097  if(fReadMC && fWriteOnlySignal){
1098  if(isMCSignal>=0) fNtupleDs->Fill(tmp);
1099  }else{
1100  fNtupleDs->Fill(tmp);
1101  }
1102  PostData(4,fNtupleDs);
1103  }
1104  }
1105 
1106  if(unsetvtx) d->UnsetOwnPrimaryVtx();
1107  if(recVtx)fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
1108  }
1109 
1110  fCounter->StoreCandidates(aod,nFiltered,kTRUE);
1111  fCounter->StoreCandidates(aod,nSelected,kFALSE);
1112 
1113  delete vHF;
1114 
1115  PostData(1,fOutput);
1116  PostData(3,fCounter);
1117 
1118  return;
1119 }
1120 
1121 //_________________________________________________________________
1122 
1123 void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
1124 {
1126  //
1127  if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
1128  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1129  if (!fOutput) {
1130  printf("ERROR: fOutput not available\n");
1131  return;
1132  }
1133  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
1134  if(fHistNEvents){
1135  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1136  }else{
1137  printf("ERROR: fHistNEvents not available\n");
1138  return;
1139  }
1140  return;
1141 }
1142 
1143 //_________________________________________________________________
1144 void AliAnalysisTaskSEDs::FillMCGenAccHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader){
1146 
1147  Int_t nProng = 3;
1148  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
1149  if(TMath::Abs(zMCVertex) <= fAnalysisCuts->GetMaxVtxZ()) {
1150  for(Int_t iPart=0; iPart<arrayMC->GetEntriesFast(); iPart++){
1151 
1152  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
1153 
1154  if(TMath::Abs(mcPart->GetPdgCode()) == 431) {
1155  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,kTRUE);//Prompt = 4, FeedDown = 5
1156 
1157  Int_t deca = 0;
1158  Bool_t isGoodDecay = kFALSE;
1159  Int_t labDau[3] = {-1,-1,-1};
1160  Bool_t isFidAcc = kFALSE;
1161  Bool_t isDaugInAcc = kFALSE;
1162 
1163  deca = AliVertexingHFUtils::CheckDsDecay(arrayMC,mcPart,labDau);
1164  if(deca == 1) isGoodDecay=kTRUE; // == 1 -> Phi pi -> kkpi
1165 
1166  if(labDau[0]==-1) continue; //protection against unfilled array of labels
1167 
1168  if(isGoodDecay) {
1169  Double_t pt = mcPart->Pt();
1170  Double_t rapid = mcPart->Y();
1171  isFidAcc = fAnalysisCuts->IsInFiducialAcceptance(pt,rapid);
1172  isDaugInAcc = CheckDaugAcc(arrayMC,nProng,labDau);
1173 
1174  if(isFidAcc) {
1175  Double_t var4nSparseAcc[2] = {pt,rapid};
1176  if(isDaugInAcc) {
1177  if(orig==4) fnSparseMC[0]->Fill(var4nSparseAcc);
1178  if(orig==5) fnSparseMC[1]->Fill(var4nSparseAcc);
1179  }
1180  }
1181  }
1182  }
1183  }
1184  }
1185 }
1186 //_________________________________________________________________
1187 Bool_t AliAnalysisTaskSEDs::CheckDaugAcc(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1189 
1190  for (Int_t iProng = 0; iProng<nProng; iProng++){
1191  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1192  if(!mcPartDaughter) {
1193  return kFALSE;
1194  }
1195  Double_t eta = mcPartDaughter->Eta();
1196  Double_t pt = mcPartDaughter->Pt();
1197  if (TMath::Abs(eta) > 0.9 || pt < 0.1) {
1198  return kFALSE;
1199  }
1200  }
1201  return kTRUE;
1202 }
1203 
1204 
1205 
1206 
1207 
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)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
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)
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
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)
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]