AliPhysics  32b88a8 (32b88a8)
 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 different number of events in the AOD and deltaAOD
564  // In case of discrepancy the event is rejected.
566  fHistNEvents->Fill(2);
567  return;
568  }
569  fHistNEvents->Fill(1);
570  }
571 
572  TClonesArray *array3Prong = 0;
573  if(!aod && AODEvent() && IsStandardAOD()) {
574  // In case there is an AOD handler writing a standard AOD, use the AOD
575  // event in memory rather than the input (ESD) event.
576  aod = dynamic_cast<AliAODEvent*> (AODEvent());
577  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
578  // have to taken from the AOD event hold by the AliAODExtension
579  AliAODHandler* aodHandler = (AliAODHandler*)
580  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
581  if(aodHandler->GetExtensions()) {
582  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
583  AliAODEvent *aodFromExt = ext->GetAOD();
584  array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
585  }
586  } else if(aod) {
587  array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
588  }
589 
590  if(!aod || !array3Prong) {
591  printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
592  return;
593  }
594 
595 
596  // fix for temporary bug in ESDfilter
597  // the AODs with null vertex pointer didn't pass the PhysSel
598  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
599 
600 
601  fHistNEvents->Fill(3); // count event
602  // Post the data already here
603  PostData(1,fOutput);
604 
606 
607 
608  Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
609  Float_t ntracks=aod->GetNumberOfTracks();
610  Float_t evCentr=fAnalysisCuts->GetCentrality(aod);
611 
612  fHistCentrality[0]->Fill(evCentr);
613  fHistCentralityMult[0]->Fill(ntracks,evCentr);
620  fHistNEvents->Fill(10);
621  fHistCentrality[2]->Fill(evCentr);
622  fHistCentralityMult[2]->Fill(ntracks,evCentr);
623  }
624 
626  Int_t runNumber=aod->GetRunNumber();
627 
628 
629 
630  TClonesArray *arrayMC=0;
631  AliAODMCHeader *mcHeader=0;
632 
633  // AOD primary vertex
634  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
635  // vtx1->Print();
636 
637  // load MC particles
638  if(fReadMC){
639 
640  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
641  if(!arrayMC) {
642  printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
643  return;
644  }
645 
646  // load MC header
647  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
648  if(!mcHeader) {
649  printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
650  return;
651  }
652  }
653 
654 
655  if(fReadMC && fFillSparse){
656  if(aod->GetTriggerMask()==0 && (runNumber>=195344 && runNumber<=195677))
657  // protection for events with empty trigger mask in p-Pb
658  return;
660  // events not passing the centrality selection can be removed immediately.
661  return;
662  Double_t zMCVertex = mcHeader->GetVtxZ();
663  if (TMath::Abs(zMCVertex) > fAnalysisCuts->GetMaxVtxZ())
664  return;
665  FillMCGenAccHistos(arrayMC, mcHeader);
666  }
667 
668  if(!isEvSel)return;
669  fHistNEvents->Fill(4);
670  fHistCentrality[1]->Fill(evCentr);
671  fHistCentralityMult[1]->Fill(ntracks,evCentr);
672 
673  Int_t n3Prong = array3Prong->GetEntriesFast();
674  if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
675 
676 
677  Int_t pdgDstoKKpi[3]={321,321,211};
678  Int_t nSelected=0;
679  Int_t nFiltered=0;
680  Double_t massPhi=TDatabasePDG::Instance()->GetParticle(333)->Mass();
681 
682  // vHF object is needed to call the method that refills the missing info of the candidates
683  // if they have been deleted in dAOD reconstruction phase
684  // in order to reduce the size of the file
686 
687  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
688 
689 
690  AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
691  fHistNEvents->Fill(11);
692 
694  continue;
695  }
696  nFiltered++;
697  fHistNEvents->Fill(12);
698 
699  if(!(vHF->FillRecoCand(aod,d))) {
700  fHistNEvents->Fill(14); //monitor how often this fails
701  continue;
702  }
703 
704  Bool_t unsetvtx=kFALSE;
705  if(!d->GetOwnPrimaryVtx()){
706  d->SetOwnPrimaryVtx(vtx1);
707  unsetvtx=kTRUE;
708  // NOTE: the ow primary vertex should be unset, otherwise there is a memory leak
709  // Pay attention if you use continue inside this loop!!!
710  }
711 
712  Bool_t recVtx=kFALSE;
713  AliAODVertex *origownvtx=0x0;
714 
715  Double_t ptCand = d->Pt();
716  Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
717  Double_t rapid=d->YDs();
718  fYVsPt->Fill(ptCand,rapid);
719  Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);
720 
721  if(isFidAcc){
722 
723  Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
724  Int_t retCodeNoRes=retCodeAnalysisCuts;
726  if(origRes){
728  retCodeNoRes=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
730  }
731  Double_t massKK = 0.;
732 
733  if(retCodeNoRes&1){ //KKpi
734  massKK=d->InvMass2Prongs(0,1,321,321);
735  Double_t massKp=d->InvMass2Prongs(1,2,321,211);
736  fMassHistKK[iPtBin]->Fill(massKK);
737  fMassHistKpi[iPtBin]->Fill(massKp);
738  }
739  if(retCodeNoRes&2){ //piKK
740  massKK=d->InvMass2Prongs(1,2,321,321);
741  Double_t massKp=d->InvMass2Prongs(0,1,211,321);
742  fMassHistKK[iPtBin]->Fill(massKK);
743  fMassHistKpi[iPtBin]->Fill(massKp);
744  }
745 
746  Int_t isKKpi=retCodeAnalysisCuts&1;
747  Int_t ispiKK=retCodeAnalysisCuts&2;
748  Int_t isPhiKKpi=retCodeAnalysisCuts&4;
749  Int_t isPhipiKK=retCodeAnalysisCuts&8;
750  Int_t isK0starKKpi=retCodeAnalysisCuts&16;
751  Int_t isK0starpiKK=retCodeAnalysisCuts&32;
752 
753  if(retCodeAnalysisCuts>0){
755  if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
756  if(fAnalysisCuts->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
757  else fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
758  }
759 
760 
761  fHistNEvents->Fill(13);
762  nSelected++;
763 
764  Int_t index=GetHistoIndex(iPtBin);
765  fPtCandHist[index]->Fill(ptCand);
766 
767 
768  Double_t weightKKpi=1.;
769  Double_t weightpiKK=1.;
771  weightKKpi=fAnalysisCuts->GetWeightForKKpi();
772  weightpiKK=fAnalysisCuts->GetWeightForpiKK();
773  if(weightKKpi>1. || weightKKpi<0.) weightKKpi=0.;
774  if(weightpiKK>1. || weightpiKK<0.) weightpiKK=0.;
775  }
776 
777  fChanHist[0]->Fill(retCodeAnalysisCuts);
778 
779 
780  Double_t invMass = 0.;
781  Int_t indexMCKKpi=-1;
782  Int_t indexMCpiKK=-1;
783  Int_t labDs=-1;
784  Int_t pdgCode0=-999;
785  Int_t isMCSignal=-1;
786 
787 
788  if(fReadMC){
789 
790  labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
791  if(labDs>=0){
792  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
793  AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
794  pdgCode0=TMath::Abs(p->GetPdgCode());
795 
796  if(isKKpi){
797  if(pdgCode0==321) {
798  indexMCKKpi=GetSignalHistoIndex(iPtBin);
799  fYVsPtSig->Fill(ptCand,rapid);
800  fChanHist[1]->Fill(retCodeAnalysisCuts);
801  isMCSignal=1;
802  }else{
803  indexMCKKpi=GetReflSignalHistoIndex(iPtBin);
804  fChanHist[3]->Fill(retCodeAnalysisCuts);
805  isMCSignal=0;
806  }
807  }
808  if(ispiKK){
809  if(pdgCode0==211) {
810  indexMCpiKK=GetSignalHistoIndex(iPtBin);
811  fYVsPtSig->Fill(ptCand,rapid);
812  fChanHist[1]->Fill(retCodeAnalysisCuts);
813  isMCSignal=1;
814  }else{
815  indexMCpiKK=GetReflSignalHistoIndex(iPtBin);
816  fChanHist[3]->Fill(retCodeAnalysisCuts);
817  isMCSignal=0;
818  }
819  }
820  }else{
821  indexMCpiKK=GetBackgroundHistoIndex(iPtBin);
822  indexMCKKpi=GetBackgroundHistoIndex(iPtBin);
823  fChanHist[2]->Fill(retCodeAnalysisCuts);
824  }
825  }
826 
827  Double_t candType = 0.5; //for bkg
828 
829  if(isKKpi){
830  invMass=d->InvMassDsKKpi();
831  fMassHist[index]->Fill(invMass,weightKKpi);
832  fPtVsMass->Fill(invMass,ptCand,weightKKpi);
833  if(isPhiKKpi){
834  fMassHistPhi[index]->Fill(invMass,weightKKpi);
835  fPtVsMassPhi->Fill(invMass,ptCand,weightKKpi);
836  }
837  if(isK0starKKpi){
838  fMassHistK0st[index]->Fill(invMass,weightKKpi);
839  fPtVsMassK0st->Fill(invMass,ptCand,weightKKpi);
840  }
841  if(fReadMC && indexMCKKpi!=-1){
842  fMassHist[indexMCKKpi]->Fill(invMass,weightKKpi);
843  if(isPhiKKpi) {
844  fMassHistPhi[indexMCKKpi]->Fill(invMass,weightKKpi);
845  if(fFillSparse) {
846  if(indexMCKKpi==GetSignalHistoIndex(iPtBin)) {
847  AliAODMCParticle *partDs = (AliAODMCParticle*)arrayMC->At(labDs);
848  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partDs,kTRUE);
849  if(orig==4) {
850  candType = 1.5;
851  }
852  if(orig==5) {
853  candType = 2.5;
854  }
855  }
856  }
857  }
858  if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass,weightKKpi);
859  }
860  }
861  if(ispiKK){
862  invMass=d->InvMassDspiKK();
863  fMassHist[index]->Fill(invMass,weightpiKK);
864  fPtVsMass->Fill(invMass,ptCand,weightpiKK);
865  if(isPhipiKK){
866  fMassHistPhi[index]->Fill(invMass,weightpiKK);
867  fPtVsMassPhi->Fill(invMass,ptCand,weightpiKK);
868  }
869  if(isK0starpiKK){
870  fMassHistK0st[index]->Fill(invMass,weightpiKK);
871  fPtVsMassK0st->Fill(invMass,ptCand,weightpiKK);
872  }
873  if(fReadMC && indexMCpiKK!=-1){
874  fMassHist[indexMCpiKK]->Fill(invMass,weightpiKK);
875  if(isPhipiKK) {
876  fMassHistPhi[indexMCpiKK]->Fill(invMass,weightpiKK);
877  if(fFillSparse) {
878  if(indexMCpiKK==GetSignalHistoIndex(iPtBin)) {
879  AliAODMCParticle *partDs = (AliAODMCParticle*)arrayMC->At(labDs);
880  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partDs,kTRUE);
881  if(orig==4) {
882  candType = 1.5;
883  }
884  if(orig==5) {
885  candType = 2.5;
886  }
887  }
888  }
889  }
890  if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass,weightpiKK);
891  }
892  }
893 
894 
896 
897  if(fFillSparse) {
898 
899  const Int_t nProng = 3;
900  Double_t deltaMassKK=999.;
901  Double_t dlen=d->DecayLength();
902  Double_t dlenxy=d->DecayLengthXY();
903  Double_t normdl=d->NormalizedDecayLength();
904  Double_t normdlxy=d->NormalizedDecayLengthXY();
905  Double_t cosp=d->CosPointingAngle();
906  Double_t cospxy=d->CosPointingAngleXY();
907  Double_t pt0=d->PtProng(0);
908  Double_t pt1=d->PtProng(1);
909  Double_t pt2=d->PtProng(2);
910  Double_t sigvert=d->GetSigmaVert();
911  Double_t cosPiDs=-99.;
912  Double_t cosPiKPhi=-99.;
913  Double_t normIP; //to store the maximum topomatic var. among the 3 prongs
914  Double_t normIPprong[nProng]; //to store IP of k,k,pi
915  if(isPhiKKpi) {
916  Double_t tmpNormIP[nProng];
917 
918  invMass = d->InvMassDsKKpi();
919  massKK = d->InvMass2Prongs(0,1,321,321);
920  deltaMassKK = TMath::Abs(massKK-massPhi);
921  cosPiDs = d->CosPiDsLabFrameKKpi();
922  cosPiKPhi = d->CosPiKPhiRFrameKKpi();
923  cosPiKPhi = TMath::Abs(cosPiKPhi*cosPiKPhi*cosPiKPhi);
924 
925  for(Int_t ip=0; ip<nProng; ip++) {
926  Double_t diffIP, errdiffIP;
927  d->Getd0MeasMinusExpProng(ip,aod->GetMagneticField(),diffIP,errdiffIP);
928  tmpNormIP[ip] = diffIP/errdiffIP;
929  if(ip==0) normIP = tmpNormIP[ip];
930  else if(TMath::Abs(tmpNormIP[ip])>TMath::Abs(normIP)) normIP = tmpNormIP[ip];
931  }
932  normIPprong[0] = tmpNormIP[0];
933  normIPprong[1] = tmpNormIP[1];
934  normIPprong[2] = tmpNormIP[2];
935 
936  Double_t var4nSparse[16] = {invMass,ptCand,deltaMassKK,dlen,dlenxy,normdl,normdlxy,cosp,cospxy,
937  sigvert,cosPiDs,cosPiKPhi,normIP,normIPprong[0],normIPprong[1],normIPprong[2]};
938 
939  if(!fReadMC) {
940  fnSparse->Fill(var4nSparse);
941  }
942  else {
943  if(indexMCKKpi==GetSignalHistoIndex(iPtBin)) {
944  if(candType==1.5) fnSparseMC[2]->Fill(var4nSparse);
945  if(candType==2.5) fnSparseMC[3]->Fill(var4nSparse);
946  }
947  }
948  }
949  if(isPhipiKK) {
950  Double_t tmpNormIP[nProng];
951 
952  invMass = d->InvMassDspiKK();
953  massKK = d->InvMass2Prongs(1,2,321,321);
954  deltaMassKK = TMath::Abs(massKK-massPhi);
955  cosPiDs = d->CosPiDsLabFramepiKK();
956  cosPiKPhi = d->CosPiKPhiRFramepiKK();
957  cosPiKPhi = TMath::Abs(cosPiKPhi*cosPiKPhi*cosPiKPhi);
958 
959  for(Int_t ip=0; ip<nProng; ip++) {
960  Double_t diffIP, errdiffIP;
961  d->Getd0MeasMinusExpProng(ip,aod->GetMagneticField(),diffIP,errdiffIP);
962  tmpNormIP[ip] = diffIP/errdiffIP;
963  if(ip==0) normIP = tmpNormIP[ip];
964  else if(TMath::Abs(tmpNormIP[ip])>TMath::Abs(normIP)) normIP = tmpNormIP[ip];
965  }
966 
967  normIPprong[0] = tmpNormIP[2];
968  normIPprong[1] = tmpNormIP[1];
969  normIPprong[2] = tmpNormIP[0];
970 
971  Double_t var4nSparse[16] = {invMass,ptCand,deltaMassKK,dlen,dlenxy,normdl,normdlxy,cosp,cospxy,
972  sigvert,cosPiDs,cosPiKPhi,normIP,normIPprong[0],normIPprong[1],normIPprong[2]};
973 
974  if(!fReadMC) {
975  fnSparse->Fill(var4nSparse);
976  }
977  else {
978  if(indexMCpiKK==GetSignalHistoIndex(iPtBin)) {
979  if(candType==1.5) fnSparseMC[2]->Fill(var4nSparse);
980  if(candType==2.5) fnSparseMC[3]->Fill(var4nSparse);
981  }
982  }
983  }
984 
985  if(fReadMC && (isPhiKKpi || isPhiKKpi)) {
986  Double_t var[6] = {ptCand,normIP,normIPprong[0],normIPprong[1],normIPprong[2],candType};
987  fnSparseIP->Fill(var);
988  }
989  }
991 
992  if(fDoCutVarHistos){
993  Double_t dlen=d->DecayLength();
994  Double_t cosp=d->CosPointingAngle();
995  Double_t pt0=d->PtProng(0);
996  Double_t pt1=d->PtProng(1);
997  Double_t pt2=d->PtProng(2);
998  Double_t sigvert=d->GetSigmaVert();
999  Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
1000  Double_t dca=d->GetDCA();
1001  Double_t ptmax=0;
1002  for(Int_t i=0;i<3;i++){
1003  if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
1004  }
1005  fCosPHist[index]->Fill(cosp);
1006  fDLenHist[index]->Fill(dlen);
1007  fSigVertHist[index]->Fill(sigvert);
1008  fSumd02Hist[index]->Fill(sumD02);
1009  fPtMaxHist[index]->Fill(ptmax);
1010  fDCAHist[index]->Fill(dca);
1011  fPtProng0Hist[index]->Fill(pt0);
1012  fPtProng1Hist[index]->Fill(pt1);
1013  fPtProng2Hist[index]->Fill(pt2);
1014  if(isKKpi){
1015  Double_t massKK=d->InvMass2Prongs(0,1,321,321);
1016  Double_t massKp=d->InvMass2Prongs(1,2,321,211);
1017  fDalitz[index]->Fill(massKK,massKp);
1018  if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);
1019  if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);
1020  if(fReadMC && indexMCKKpi!=-1){
1021  fDalitz[indexMCKKpi]->Fill(massKK,massKp);
1022  if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);
1023  if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);
1024  fCosPHist[indexMCKKpi]->Fill(cosp);
1025  fDLenHist[indexMCKKpi]->Fill(dlen);
1026  fSigVertHist[indexMCKKpi]->Fill(sigvert);
1027  fSumd02Hist[indexMCKKpi]->Fill(sumD02);
1028  fPtMaxHist[indexMCKKpi]->Fill(ptmax);
1029  fPtCandHist[indexMCKKpi]->Fill(ptCand);
1030  fDCAHist[indexMCKKpi]->Fill(dca);
1031  fPtProng0Hist[indexMCKKpi]->Fill(pt0);
1032  fPtProng1Hist[indexMCKKpi]->Fill(pt1);
1033  fPtProng2Hist[indexMCKKpi]->Fill(pt2);
1034  }
1035  }
1036  if(ispiKK){
1037  Double_t massKK=d->InvMass2Prongs(1,2,321,321);
1038  Double_t massKp=d->InvMass2Prongs(0,1,211,321);
1039  fDalitz[index]->Fill(massKK,massKp);
1040  if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);
1041  if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);
1042 
1043 
1044  if(fReadMC && indexMCpiKK!=-1){
1045  fDalitz[indexMCpiKK]->Fill(massKK,massKp);
1046  if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);
1047  if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);
1048  fCosPHist[indexMCpiKK]->Fill(cosp);
1049  fDLenHist[indexMCpiKK]->Fill(dlen);
1050  fSigVertHist[indexMCpiKK]->Fill(sigvert);
1051  fSumd02Hist[indexMCpiKK]->Fill(sumD02);
1052  fPtMaxHist[indexMCpiKK]->Fill(ptmax);
1053  fPtCandHist[indexMCpiKK]->Fill(ptCand);
1054  fDCAHist[indexMCpiKK]->Fill(dca);
1055  fPtProng0Hist[indexMCpiKK]->Fill(pt0);
1056  fPtProng1Hist[indexMCpiKK]->Fill(pt1);
1057  fPtProng2Hist[indexMCpiKK]->Fill(pt2);
1058  }
1059  }
1060 
1061 
1062  }
1063 
1064  Float_t tmp[37];
1065  if(fFillNtuple>0){
1066 
1067  if ((fFillNtuple==1 && (isPhiKKpi || isPhipiKK)) || (fFillNtuple==2 && (isK0starKKpi || isK0starpiKK)) || (fFillNtuple==3 && (isKKpi || ispiKK))){
1068 
1069  AliAODTrack *track0=(AliAODTrack*)d->GetDaughter(0);
1070  AliAODTrack *track1=(AliAODTrack*)d->GetDaughter(1);
1071  AliAODTrack *track2=(AliAODTrack*)d->GetDaughter(2);
1072 
1073  UInt_t bitMapPIDTrack0=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track0);
1074  UInt_t bitMapPIDTrack1=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track1);
1075  UInt_t bitMapPIDTrack2=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track2);
1076 
1077  tmp[0]=Float_t(labDs);
1078  if(fReadMC && fWriteOnlySignal) tmp[0]=Float_t(isMCSignal);
1079  tmp[1]=Float_t(retCodeAnalysisCuts);
1080  tmp[2]=Float_t(pdgCode0);
1081  tmp[3]=d->PtProng(0);
1082  tmp[4]=d->PtProng(1);
1083  tmp[5]=d->PtProng(2);
1084  tmp[6]=d->Pt();
1085  tmp[7]=d->PProng(0);
1086  tmp[8]=d->PProng(1);
1087  tmp[9]=d->PProng(2);
1088  tmp[10]=Int_t(bitMapPIDTrack0);
1089  tmp[11]=Int_t(bitMapPIDTrack1);
1090  tmp[12]=Int_t(bitMapPIDTrack2);
1091  tmp[13]=d->CosPointingAngle();
1092  tmp[14]=d->CosPointingAngleXY();
1093  tmp[15]=d->DecayLength();
1094  tmp[16]=d->DecayLengthXY();
1095  tmp[17]=d->NormalizedDecayLength();
1096  tmp[18]=d->NormalizedDecayLengthXY();
1097  tmp[19]=d->InvMassDsKKpi();
1098  tmp[20]=d->InvMassDspiKK();
1099  tmp[21]=d->GetSigmaVert();
1100  tmp[22]=d->Getd0Prong(0);
1101  tmp[23]=d->Getd0Prong(1);
1102  tmp[24]=d->Getd0Prong(2);
1103  tmp[25]=d->GetDCA();
1104  tmp[26]=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
1105  tmp[27]=d->InvMass2Prongs(0,1,321,321);
1106  tmp[28]=d->InvMass2Prongs(1,2,321,321);
1107  tmp[29]=d->InvMass2Prongs(1,2,321,211);
1108  tmp[30]=d->InvMass2Prongs(0,1,211,321);
1109  tmp[31]=d->CosPiDsLabFrameKKpi();
1110  tmp[32]=d->CosPiDsLabFramepiKK();
1111  tmp[33]=d->CosPiKPhiRFrameKKpi();
1112  tmp[34]=d->CosPiKPhiRFramepiKK();
1113  tmp[35]=(Float_t)(centrality);
1114  tmp[36]=(Float_t)(runNumber);
1115 
1116  if(fReadMC && fWriteOnlySignal){
1117  if(isMCSignal>=0) fNtupleDs->Fill(tmp);
1118  }else{
1119  fNtupleDs->Fill(tmp);
1120  }
1121  PostData(4,fNtupleDs);
1122  }
1123  }
1124  } //if(retCodeAnalysisCuts
1125  } // if(isFidAcc)
1126 
1127  if(unsetvtx) d->UnsetOwnPrimaryVtx();
1128  if(recVtx)fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
1129  }
1130 
1131  fCounter->StoreCandidates(aod,nFiltered,kTRUE);
1132  fCounter->StoreCandidates(aod,nSelected,kFALSE);
1133 
1134  delete vHF;
1135 
1136  PostData(1,fOutput);
1137  PostData(3,fCounter);
1138 
1139  return;
1140 }
1141 
1142 //_________________________________________________________________
1143 
1145 {
1147  //
1148  if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
1149  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1150  if (!fOutput) {
1151  printf("ERROR: fOutput not available\n");
1152  return;
1153  }
1154  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
1155  if(fHistNEvents){
1156  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1157  }else{
1158  printf("ERROR: fHistNEvents not available\n");
1159  return;
1160  }
1161  return;
1162 }
1163 
1164 //_________________________________________________________________
1165 void AliAnalysisTaskSEDs::FillMCGenAccHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader){
1167 
1168  Int_t nProng = 3;
1169  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
1170  if(TMath::Abs(zMCVertex) <= fAnalysisCuts->GetMaxVtxZ()) {
1171  for(Int_t iPart=0; iPart<arrayMC->GetEntriesFast(); iPart++){
1172 
1173  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
1174 
1175  if(TMath::Abs(mcPart->GetPdgCode()) == 431) {
1176  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,kTRUE);//Prompt = 4, FeedDown = 5
1177 
1178  Int_t deca = 0;
1179  Bool_t isGoodDecay = kFALSE;
1180  Int_t labDau[3] = {-1,-1,-1};
1181  Bool_t isFidAcc = kFALSE;
1182  Bool_t isDaugInAcc = kFALSE;
1183 
1184  deca = AliVertexingHFUtils::CheckDsDecay(arrayMC,mcPart,labDau);
1185  if(deca == 1) isGoodDecay=kTRUE; // == 1 -> Phi pi -> kkpi
1186 
1187  if(labDau[0]==-1) continue; //protection against unfilled array of labels
1188 
1189  if(isGoodDecay) {
1190  Double_t pt = mcPart->Pt();
1191  Double_t rapid = mcPart->Y();
1192  isFidAcc = fAnalysisCuts->IsInFiducialAcceptance(pt,rapid);
1193  isDaugInAcc = CheckDaugAcc(arrayMC,nProng,labDau);
1194 
1195  if(isFidAcc) {
1196  Double_t var4nSparseAcc[2] = {pt,rapid};
1197  if(isDaugInAcc) {
1198  if(orig==4) fnSparseMC[0]->Fill(var4nSparseAcc);
1199  if(orig==5) fnSparseMC[1]->Fill(var4nSparseAcc);
1200  }
1201  }
1202  }
1203  }
1204  }
1205  }
1206 }
1207 //_________________________________________________________________
1208 Bool_t AliAnalysisTaskSEDs::CheckDaugAcc(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1210 
1211  for (Int_t iProng = 0; iProng<nProng; iProng++){
1212  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1213  if(!mcPartDaughter) {
1214  return kFALSE;
1215  }
1216  Double_t eta = mcPartDaughter->Eta();
1217  Double_t pt = mcPartDaughter->Pt();
1218  if (TMath::Abs(eta) > 0.9 || pt < 0.1) {
1219  return kFALSE;
1220  }
1221  }
1222  return kTRUE;
1223 }
1224 
1225 
1226 
1227 
1228 
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
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)
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:311
TH1F * fMassHistKpi[kMaxPtBins]
! hist. of mass spectra of Kpi
static Bool_t CheckMatchingAODdeltaAODevents()
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: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)
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]