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