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