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