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