AliPhysics  ec7afe5 (ec7afe5)
 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  fHistAllV0multNTPCout(0),
66  fHistSelV0multNTPCout(0),
67  fCosPHist3D(0x0),
68  fCosPxyHist3D(0x0),
69  fDLenHist3D(0x0),
70  fDLenxyHist3D(0x0),
71  fNDLenxyHist3D(0x0),
72  fSigVertHist3D(0x0),
73  fDCAHist3D(0x0),
74  fNormIPHist3D(0x0),
75  fCosPiDsHist3D(0x0),
76  fCosPiKPhiHist3D(0x0),
77  fPtProng0Hist3D(0x0),
78  fPtProng1Hist3D(0x0),
79  fPtProng2Hist3D(0x0),
80  fNtupleDs(0),
81  fFillNtuple(0),
82  fSystem(0),
83  fReadMC(kFALSE),
84  fWriteOnlySignal(kFALSE),
85  fDoCutVarHistos(kTRUE),
86  fUseSelectionBit(kFALSE),
87  fFillSparse(kTRUE),
88  fFillSparseDplus(kFALSE),
89  fFillImpParSparse(kFALSE),
90  fDoRotBkg(kFALSE),
91  fDoBkgPhiSB(kFALSE),
92  fDoCutV0multTPCout(kFALSE),
93  fAODProtection(1),
94  fNPtBins(0),
95  fListCuts(0),
96  fMassRange(0.8),
97  fMassBinSize(0.002),
98  fminMass(1.6),
99  fmaxMass(2.5),
100  fMaxDeltaPhiMass4Rot(0.010),
101  fCounter(0),
102  fAnalysisCuts(0),
103  fnSparse(0),
104  fnSparseIP(0),
105  fImpParSparse(0x0),
106  fImpParSparseMC(0x0)
107 {
109 
110  for(Int_t i=0;i<3;i++){
111  fHistCentrality[i]=0;
112  fHistCentralityMult[i]=0;
113  }
114  for(Int_t i=0;i<4;i++) {
115  fChanHist[i]=0;
116  }
117  for(Int_t i=0;i<4*kMaxPtBins;i++){
118  fPtCandHist[i]=0;
119  fMassHist[i]=0;
120  fMassHistPhi[i]=0;
121  fMassHistK0st[i]=0;
122  fCosPHist[i]=0;
123  fDLenHist[i]=0;
124  fSumd02Hist[i]=0;
125  fSigVertHist[i]=0;
126  fPtMaxHist[i]=0;
127  fDCAHist[i]=0;
128  fPtProng0Hist[i]=0;
129  fPtProng1Hist[i]=0;
130  fPtProng2Hist[i]=0;
131  fDalitz[i]=0;
132  fDalitzPhi[i]=0;
133  fDalitzK0st[i]=0;
134  }
135  for(Int_t i=0;i<kMaxPtBins;i++){
136  fMassHistKK[i]=0;
137  fMassHistKpi[i]=0;
138  fMassRotBkgHistPhi[i]=0;
139  fMassLSBkgHistPhi[i]=0;
140  fMassRSBkgHistPhi[i]=0;
141  }
142  for(Int_t i=0;i<kMaxPtBins+1;i++){
143  fPtLimits[i]=0;
144  }
145  for (Int_t i=0; i<4; i++) {
146  fnSparseMC[i]=0;
147  fnSparseMCDplus[i]=0;
148  }
149 }
150 
151 //________________________________________________________________________
152 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name,AliRDHFCutsDstoKKpi* analysiscuts,Int_t fillNtuple):
153  AliAnalysisTaskSE(name),
154  fOutput(0),
155  fHistNEvents(0),
156  fPtVsMass(0),
157  fPtVsMassPhi(0),
158  fPtVsMassK0st(0),
159  fYVsPt(0),
160  fYVsPtSig(0),
161  fHistAllV0multNTPCout(0),
162  fHistSelV0multNTPCout(0),
163  fCosPHist3D(0x0),
164  fCosPxyHist3D(0x0),
165  fDLenHist3D(0x0),
166  fDLenxyHist3D(0x0),
167  fNDLenxyHist3D(0x0),
168  fSigVertHist3D(0x0),
169  fDCAHist3D(0x0),
170  fNormIPHist3D(0x0),
171  fCosPiDsHist3D(0x0),
172  fCosPiKPhiHist3D(0x0),
173  fPtProng0Hist3D(0x0),
174  fPtProng1Hist3D(0x0),
175  fPtProng2Hist3D(0x0),
176  fNtupleDs(0),
177  fFillNtuple(fillNtuple),
178  fSystem(0),
179  fReadMC(kFALSE),
180  fWriteOnlySignal(kFALSE),
181  fDoCutVarHistos(kTRUE),
182  fUseSelectionBit(kFALSE),
183  fFillSparse(kTRUE),
184  fFillSparseDplus(kFALSE),
185  fFillImpParSparse(kFALSE),
186  fDoRotBkg(kTRUE),
187  fDoBkgPhiSB(kTRUE),
188  fDoCutV0multTPCout(kFALSE),
189  fAODProtection(1),
190  fNPtBins(0),
191  fListCuts(0),
192  fMassRange(0.8),
193  fMassBinSize(0.002),
194  fminMass(1.6),
195  fmaxMass(2.5),
196  fMaxDeltaPhiMass4Rot(0.010),
197  fCounter(0),
198  fAnalysisCuts(analysiscuts),
199  fnSparse(0),
200  fnSparseIP(0),
201  fImpParSparse(0x0),
202  fImpParSparseMC(0x0)
203 {
206 
207  for(Int_t i=0;i<3;i++){
208  fHistCentrality[i]=0;
209  fHistCentralityMult[i]=0;
210  }
211  for(Int_t i=0;i<4;i++) {
212  fChanHist[i]=0;
213  }
214  for(Int_t i=0;i<4*kMaxPtBins;i++){
215  fPtCandHist[i]=0;
216  fMassHist[i]=0;
217  fMassHistPhi[i]=0;
218  fMassHistK0st[i]=0;
219  fCosPHist[i]=0;
220  fDLenHist[i]=0;
221  fSumd02Hist[i]=0;
222  fSigVertHist[i]=0;
223  fPtMaxHist[i]=0;
224  fDCAHist[i]=0;
225  fPtProng0Hist[i]=0;
226  fPtProng1Hist[i]=0;
227  fPtProng2Hist[i]=0;
228  fDalitz[i]=0;
229  fDalitzPhi[i]=0;
230  fDalitzK0st[i]=0;
231  }
232  for(Int_t i=0;i<kMaxPtBins;i++){
233  fMassHistKK[i]=0;
234  fMassHistKpi[i]=0;
235  fMassRotBkgHistPhi[i]=0;
236  fMassLSBkgHistPhi[i]=0;
237  fMassRSBkgHistPhi[i]=0;
238  }
239  for(Int_t i=0;i<kMaxPtBins+1;i++){
240  fPtLimits[i]=0;
241  }
242 
243  for (Int_t i=0; i<4; i++) {
244  fnSparseMC[i]=0;
245  fnSparseMCDplus[i]=0;
246  }
247 
250  SetPtBins(nptbins,ptlim);
251 
252  DefineOutput(1,TList::Class()); //My private output
253 
254  DefineOutput(2,TList::Class());
255 
256  DefineOutput(3,AliNormalizationCounter::Class());
257 
258  if(fFillNtuple>0){
259  // Output slot #4 writes into a TNtuple container
260  DefineOutput(4,TNtuple::Class()); //My private output
261  }
262 
263 }
264 
265 //________________________________________________________________________
268  if(n>kMaxPtBins){
269  printf("Max. number of Pt bins = %d\n",kMaxPtBins);
271  fPtLimits[0]=0.;
272  fPtLimits[1]=1.;
273  fPtLimits[2]=3.;
274  fPtLimits[3]=5.;
275  fPtLimits[4]=10.;
276  for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
277  }else{
278  fNPtBins=n;
279  for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
280  for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
281  }
282  if(fDebug > 1){
283  printf("Number of Pt bins = %d\n",fNPtBins);
284  for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);
285  }
286 }
287 //________________________________________________________________________
289 {
290  // Destructor
291  if(fOutput && !fOutput->IsOwner()){
292  delete fHistNEvents;
293  delete fHistAllV0multNTPCout;
294  delete fHistSelV0multNTPCout;
295  delete fCosPHist3D;
296  delete fCosPxyHist3D;
297  delete fDLenHist3D;
298  delete fDLenxyHist3D;
299  delete fNDLenxyHist3D;
300  delete fSigVertHist3D;
301  delete fDCAHist3D;
302  delete fNormIPHist3D;
303  delete fCosPiDsHist3D;
304  delete fCosPiKPhiHist3D;
305  delete fPtProng0Hist3D;
306  delete fPtProng1Hist3D;
307  delete fPtProng2Hist3D;
308 
309  for(Int_t i=0;i<4;i++){
310  delete fChanHist[i];
311  }
312  for(Int_t i=0;i<4*fNPtBins;i++){
313  delete fMassHist[i];
314  delete fMassHistPhi[i];
315  delete fMassHistK0st[i];
316  delete fCosPHist[i];
317  delete fDLenHist[i];
318  delete fSumd02Hist[i];
319  delete fSigVertHist[i];
320  delete fPtMaxHist[i];
321  delete fPtCandHist[i];
322  delete fDCAHist[i];
323  delete fPtProng0Hist[i];
324  delete fPtProng1Hist[i];
325  delete fPtProng2Hist[i];
326  delete fDalitz[i];
327  delete fDalitzPhi[i];
328  delete fDalitzK0st[i];
329  }
330  for(Int_t i=0;i<fNPtBins;i++){
331  delete fMassHistKK[i];
332  delete fMassHistKpi[i];
333  delete fMassRotBkgHistPhi[i];
334  delete fMassLSBkgHistPhi[i];
335  delete fMassRSBkgHistPhi[i];
336  }
337  delete fPtVsMass;
338  delete fPtVsMassPhi;
339  delete fPtVsMassK0st;
340  delete fYVsPt;
341  delete fYVsPtSig;
342  for(Int_t i=0;i<3;i++){
343  delete fHistCentrality[i];
344  delete fHistCentralityMult[i];
345  }
346 
347  delete fnSparse;
348  delete fnSparseIP;
349  if(fFillImpParSparse) {
350  delete fImpParSparse;
351  delete fImpParSparseMC;
352  }
353  for (Int_t i=0; i<4; i++) {
354  delete fnSparseMC[i];
355  if(fFillSparseDplus)delete fnSparseMCDplus[i];
356  }
357  }
358  delete fOutput;
359  delete fListCuts;
360  delete fNtupleDs;
361  delete fCounter;
362  delete fAnalysisCuts;
363 
364 }
365 
366 //________________________________________________________________________
368 {
370 
371  if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
372 
373  fListCuts=new TList();
374  fListCuts->SetOwner();
375  fListCuts->SetName("CutObjects");
376 
378  analysis->SetName("AnalysisCuts");
379 
380  fListCuts->Add(analysis);
381  PostData(2,fListCuts);
382  return;
383 }
384 
385 //________________________________________________________________________
387 {
389  //
390  if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
391 
392  // Several histograms are more conveniently managed in a TList
393  fOutput = new TList();
394  fOutput->SetOwner();
395  fOutput->SetName("OutputHistos");
396 
397  fHistNEvents = new TH1F("hNEvents", "number of events ",15,-0.5,14.5);
398  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsRead");
399  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents Matched dAOD");
400  fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents Mismatched dAOD");
401  fHistNEvents->GetXaxis()->SetBinLabel(4,"nEventsAnal");
402  fHistNEvents->GetXaxis()->SetBinLabel(5,"n. passing IsEvSelected");
403  fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected due to trigger");
404  fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected due to not reco vertex");
405  fHistNEvents->GetXaxis()->SetBinLabel(8,"n. rejected for contr vertex");
406  fHistNEvents->GetXaxis()->SetBinLabel(9,"n. rejected for vertex out of accept");
407  fHistNEvents->GetXaxis()->SetBinLabel(10,"n. rejected for pileup events");
408  fHistNEvents->GetXaxis()->SetBinLabel(11,"no. of out centrality events");
409  fHistNEvents->GetXaxis()->SetBinLabel(12,"no. of 3 prong candidates");
410  fHistNEvents->GetXaxis()->SetBinLabel(13,"no. of Ds after filtering cuts");
411  fHistNEvents->GetXaxis()->SetBinLabel(14,"no. of Ds after selection cuts");
412  fHistNEvents->GetXaxis()->SetBinLabel(15,"no. of not on-the-fly rec Ds");
413 
414  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
415 
416  fHistNEvents->Sumw2();
417  fHistNEvents->SetMinimum(0);
418  fOutput->Add(fHistNEvents);
419 
420  fHistCentrality[0]=new TH1F("hCentr","centrality",10000,0.,100.);
421  fHistCentrality[1]=new TH1F("hCentr(selectedCent)","centrality(selectedCent)",10000,0.,100.);
422  fHistCentrality[2]=new TH1F("hCentr(OutofCent)","centrality(OutofCent)",10000,0.,100.);
423  fHistCentralityMult[0]=new TH2F("hCentrMult","centrality vs mult",100,0.5,30000.5,40,0.,100.);
424  fHistCentralityMult[1]=new TH2F("hCentrMult(selectedCent)","centrality vs mult(selectedCent)",100,0.5,30000.5,40,0.,100.);
425  fHistCentralityMult[2]=new TH2F("hCentrMult(OutofCent)","centrality vs mult(OutofCent)",100,0.5,30000.5,40,0.,100.);
426  for(Int_t i=0;i<3;i++){
427  fHistCentrality[i]->Sumw2();
428  fOutput->Add(fHistCentrality[i]);
429  fHistCentralityMult[i]->Sumw2();
430  fOutput->Add(fHistCentralityMult[i]);
431  }
432  if(fDoCutV0multTPCout) {
433  fHistAllV0multNTPCout = new TH2F("HistAllV0multNTPCout", "V0mult vs # TPCout (all) ;V0mult ;# TPCout", 1000, 0., 40000, 1000, 0, 30000);
434  fHistSelV0multNTPCout = new TH2F("HistSelV0multNTPCout", "V0mult vs # TPCout (sel) ;V0mult ;# TPCout", 1000, 0., 40000, 1000, 0, 30000);
435  fHistAllV0multNTPCout->Sumw2();
436  fHistSelV0multNTPCout->Sumw2();
439  }
440 
441  Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
442 
443  Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);
444  if(nInvMassBins%2==1) nInvMassBins++;
445  // Double_t minMass=massDs-1.0*nInvMassBins*fMassBinSize;
446  Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
447  // Double_t maxMass=massDs+1.0*nInvMassBins*fMassBinSize;
448  Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
449  fminMass = minMass;
450  fmaxMass = maxMass;
451 
452  TString hisname;
453  TString htype;
454  Int_t index;
455  for(Int_t iType=0; iType<4; iType++){
456  for(Int_t i=0;i<fNPtBins;i++){
457  if(iType==0){
458  htype="All";
459  index=GetHistoIndex(i);
460  }else if(iType==1){
461  htype="Sig";
462  index=GetSignalHistoIndex(i);
463  }else if(iType==2){
464  htype="Bkg";
465  index=GetBackgroundHistoIndex(i);
466  }else{
467  htype="ReflSig";
468  index=GetReflSignalHistoIndex(i);
469  }
470  hisname.Form("hMass%sPt%d",htype.Data(),i);
471  fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
472  fMassHist[index]->Sumw2();
473  hisname.Form("hMass%sPt%dphi",htype.Data(),i);
474  fMassHistPhi[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
475  fMassHistPhi[index]->Sumw2();
476  hisname.Form("hMass%sPt%dk0st",htype.Data(),i);
477  fMassHistK0st[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
478  fMassHistK0st[index]->Sumw2();
479  hisname.Form("hCosP%sPt%d",htype.Data(),i);
480  fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
481  fCosPHist[index]->Sumw2();
482  hisname.Form("hCosPxy%sPt%d",htype.Data(),i);
483  fCosPxyHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
484  fCosPxyHist[index]->Sumw2();
485  hisname.Form("hDLen%sPt%d",htype.Data(),i);
486  fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
487  fDLenHist[index]->Sumw2();
488  hisname.Form("hDLenxy%sPt%d",htype.Data(),i);
489  fDLenxyHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
490  fDLenxyHist[index]->Sumw2();
491  hisname.Form("hNDLenxy%sPt%d",htype.Data(),i);
492  fNDLenxyHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,11.);
493  fNDLenxyHist[index]->Sumw2();
494  hisname.Form("hSumd02%sPt%d",htype.Data(),i);
495  fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
496  fSumd02Hist[index]->Sumw2();
497  hisname.Form("hSigVert%sPt%d",htype.Data(),i);
498  fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
499  fSigVertHist[index]->Sumw2();
500  hisname.Form("hPtMax%sPt%d",htype.Data(),i);
501  fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
502  fPtMaxHist[index]->Sumw2();
503  hisname.Form("hPtCand%sPt%d",htype.Data(),i);
504  fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
505  fPtCandHist[index]->Sumw2();
506  hisname.Form("hDCA%sPt%d",htype.Data(),i);
507  fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
508  fDCAHist[index]->Sumw2();
509  hisname.Form("hNormIP%sPt%d",htype.Data(),i);
510  fNormIPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,6.);
511  fNormIPHist[index]->Sumw2();
512  hisname.Form("hCosPiDs%sPt%d",htype.Data(),i);
513  fCosPiDsHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
514  fCosPiDsHist[index]->Sumw2();
515  hisname.Form("hCosPiKPhi%sPt%d",htype.Data(),i);
516  fCosPiKPhiHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
517  fCosPiKPhiHist[index]->Sumw2();
518  hisname.Form("hPtProng0%sPt%d",htype.Data(),i);
519  fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
520  fPtProng0Hist[index]->Sumw2();
521  hisname.Form("hPtProng1%sPt%d",htype.Data(),i);
522  fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
523  fPtProng1Hist[index]->Sumw2();
524  hisname.Form("hPtProng2%sPt%d",htype.Data(),i);
525  fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
526  fPtProng2Hist[index]->Sumw2();
527  hisname.Form("hDalitz%sPt%d",htype.Data(),i);
528  fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
529  fDalitz[index]->Sumw2();
530  hisname.Form("hDalitz%sPt%dphi",htype.Data(),i);
531  fDalitzPhi[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
532  fDalitzPhi[index]->Sumw2();
533  hisname.Form("hDalitz%sPt%dk0st",htype.Data(),i);
534  fDalitzK0st[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
535  fDalitzK0st[index]->Sumw2();
536  }
537  }
538 
539  for(Int_t i=0; i<4*fNPtBins; i++){
540  fOutput->Add(fMassHist[i]);
541  fOutput->Add(fMassHistPhi[i]);
542  fOutput->Add(fMassHistK0st[i]);
543  fOutput->Add(fPtCandHist[i]);
544  if(fDoCutVarHistos){
545  fOutput->Add(fCosPHist[i]);
546  fOutput->Add(fCosPxyHist[i]);
547  fOutput->Add(fDLenHist[i]);
548  fOutput->Add(fDLenxyHist[i]);
549  fOutput->Add(fNDLenxyHist[i]);
550  fOutput->Add(fSumd02Hist[i]);
551  fOutput->Add(fSigVertHist[i]);
552  fOutput->Add(fPtMaxHist[i]);
553  fOutput->Add(fDCAHist[i]);
554  fOutput->Add(fNormIPHist[i]);
555  fOutput->Add(fCosPiDsHist[i]);
556  fOutput->Add(fCosPiKPhiHist[i]);
557  fOutput->Add(fPtProng0Hist[i]);
558  fOutput->Add(fPtProng1Hist[i]);
559  fOutput->Add(fPtProng2Hist[i]);
560  fOutput->Add(fDalitz[i]);
561  fOutput->Add(fDalitzPhi[i]);
562  fOutput->Add(fDalitzK0st[i]);
563  }
564  }
565 
566  fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",64,-0.5,63.5);
567  fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",64,-0.5,63.5);
568  fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",64,-0.5,63.5);
569  fChanHist[3] = new TH1F("hChanReflSig", "KKpi and piKK candidates",64,-0.5,63.5);
570  for(Int_t i=0;i<4;i++){
571  fChanHist[i]->Sumw2();
572  fChanHist[i]->SetMinimum(0);
573  fOutput->Add(fChanHist[i]);
574  }
575 
576  fCosPHist3D = new TH3F("fCosPHist3D","CosP vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.5,1.);
577  fCosPxyHist3D = new TH3F("fCosPxyHist3D","CosPxy vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.5,1.);
578  fDLenHist3D = new TH3F("fDLenHist3D","DLen vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,0.5);
579  fDLenxyHist3D = new TH3F("fDLenxyHist3D","DLenxy vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,0.5);
580  fNDLenxyHist3D = new TH3F("fNDLenxyHist3D","NDLenxy vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,11.);
581  fSigVertHist3D = new TH3F("fSigVertHist3D","SigVert vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,0.1);
582  fDCAHist3D = new TH3F("fDCAHist3D","DCA vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,0.1);
583  fNormIPHist3D = new TH3F("fNormIPHist3D","nIP vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,6.);
584  fCosPiDsHist3D = new TH3F("fCosPiDsHist3D","CosPiDs vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.5,1.);
585  fCosPiKPhiHist3D = new TH3F("fCosPiKPhiHist3D","CosPiKPhi vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,0.5);
586  fPtProng0Hist3D = new TH3F("fPtProng0Hist3D","Pt prong0 vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.0,20.);
587  fPtProng1Hist3D = new TH3F("fPtProng1Hist3D","Pt prong1 vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.0,20.);
588  fPtProng2Hist3D = new TH3F("fPtProng2Hist3D","Pt prong2 vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.0,20.);
589 
590  if(!fReadMC && fDoCutVarHistos) {
591  fOutput->Add(fCosPHist3D);
592  fOutput->Add(fCosPxyHist3D);
593  fOutput->Add(fDLenHist3D);
594  fOutput->Add(fDLenxyHist3D);
595  fOutput->Add(fNDLenxyHist3D);
596  fOutput->Add(fSigVertHist3D);
597  fOutput->Add(fDCAHist3D);
598  fOutput->Add(fNormIPHist3D);
599  fOutput->Add(fCosPiDsHist3D);
601  fOutput->Add(fPtProng0Hist3D);
602  fOutput->Add(fPtProng1Hist3D);
603  fOutput->Add(fPtProng2Hist3D);
604  }
605 
606  fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nInvMassBins,minMass,maxMass,40,0.,20.);
607  fPtVsMassPhi=new TH2F("hPtVsMassPhi","PtVsMass (phi selection)",nInvMassBins,minMass,maxMass,200,0.,20.);
608  fPtVsMassK0st=new TH2F("hPtVsMassK0st","PtVsMass (K0* selection)",nInvMassBins,minMass,maxMass,200,0.,20.);
609  fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
610  fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
611 
612  for(Int_t i=0;i<fNPtBins;i++){
613  hisname.Form("hMassKKPt%d",i);
614  fMassHistKK[i]=new TH1F(hisname.Data(),hisname.Data(),200,0.95,1.35);
615  fMassHistKK[i]->Sumw2();
616  fOutput->Add(fMassHistKK[i]);
617  hisname.Form("hMassKpiPt%d",i);
618  fMassHistKpi[i]=new TH1F(hisname.Data(),hisname.Data(),200,0.7,1.1);
619  fMassHistKpi[i]->Sumw2();
620  fOutput->Add(fMassHistKpi[i]);
621  if(fDoRotBkg) {
622  hisname.Form("hMassAllPt%dphi_RotBkg",i);
623  fMassRotBkgHistPhi[i]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
624  fMassRotBkgHistPhi[i]->Sumw2();
625  fOutput->Add(fMassRotBkgHistPhi[i]);
626  }
627  if(fDoBkgPhiSB) {
628  hisname.Form("fMassLSBkgHistPhiPt%d",i);
629  fMassLSBkgHistPhi[i]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
630  fMassLSBkgHistPhi[i]->Sumw2();
631  fOutput->Add(fMassLSBkgHistPhi[i]);
632  hisname.Form("fMassRSBkgHistPhiPt%d",i);
633  fMassRSBkgHistPhi[i]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
634  fMassRSBkgHistPhi[i]->Sumw2();
635  fOutput->Add(fMassRSBkgHistPhi[i]);
636  }
637  }
638 
639  fOutput->Add(fPtVsMass);
640  fOutput->Add(fPtVsMassPhi);
641  fOutput->Add(fPtVsMassK0st);
642  fOutput->Add(fYVsPt);
643  fOutput->Add(fYVsPtSig);
644 
645  nInvMassBins=(Int_t)(0.7/fMassBinSize+0.5);
646  minMass=massDs-0.5*nInvMassBins*fMassBinSize;
647  maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
648 
649  Int_t nBinsReco[knVarForSparse] = {nInvMassBins, 20, 30, 14, 14, 20, 10, 10, 14, 6, 6, 12};
650  Double_t xminReco[knVarForSparse] = {minMass, 0., 0., 0., 0., 0., 90., 90., 0., 7., 0., 0.};
651  Double_t xmaxReco[knVarForSparse] = {maxMass, 20., 15, 70., 70., 10., 100., 100., 70., 10., 3., 6.};
652  TString axis[knVarForSparse] = {"invMassDsAllPhi","p_{T}","#Delta Mass(KK)","dlen","dlen_{xy}","normdl_{xy}","cosP","cosP_{xy}","sigVert","cosPiDs","|cosPiKPhi^{3}|","normIP"};
653  if(fSystem == 1) { //pPb,PbPb
654  nInvMassBins=(Int_t)(0.45/fMassBinSize+0.5);
655  minMass=massDs-0.5*nInvMassBins*fMassBinSize;
656  maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
657  nBinsReco[0] = nInvMassBins; //Ds mass
658  xminReco[0] = minMass;
659  xmaxReco[0] = maxMass;
660 
661  nBinsReco[1] = 16; //pt
662  xminReco[1] = 0.;
663  xmaxReco[1] = 16.;
664 
665  nBinsReco[2] = 12; //#Delta Mass(KK)
666  xmaxReco[2] = 12.;
667 
668  nBinsReco[3] = 7; //dlen
669  nBinsReco[4] = 7; //dlenxy
670  nBinsReco[5] = 10; //ndlenxy
671 
672  nBinsReco[6] = 6; //cosP
673  xminReco[6] = 97.;
674  xmaxReco[6] = 100.;
675 
676  nBinsReco[7] = 6; //cosPxy
677  xminReco[7] = 97.;
678  xmaxReco[7] = 100.;
679  }
680 
681  Int_t nBinsAcc[knVarForSparseAcc] = {20, 20};
682  Double_t xminAcc[knVarForSparseAcc] = {0., -10.};
683  Double_t xmaxAcc[knVarForSparseAcc] = {20, 10.};
684 
685  Int_t nBinsIP[knVarForSparseIP] = { 20, 400, 400, 400, 400, 3};
686  Double_t xminIP[knVarForSparseIP] = { 0., -10., -10., -10., -10., 0.};
687  Double_t xmaxIP[knVarForSparseIP] = {20., 10., 10., 10., 10., 3.};
688  TString axisIP[knVarForSparseIP] = {"motherPt","maxNormImp","IP0","IP1","IP2","candType"};
689 
690  if(fFillSparse) {
691 
692  if(fReadMC) {
693  TString label[knVarForSparseAcc] = {"fromC","fromB"};
694  for (Int_t i=0; i<2; i++) {
695  fnSparseMC[i] = new THnSparseF(Form("fnSparseAcc_%s",label[i].Data()),Form("MC nSparse (Acc.Step)- %s",label[i].Data()),
696  knVarForSparseAcc, nBinsAcc, xminAcc, xmaxAcc);
697  fnSparseMC[i]->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
698  fnSparseMC[i]->GetAxis(1)->SetTitle("y");
699  fOutput->Add(fnSparseMC[i]);
700 
701  //Dplus
702  if(fFillSparseDplus) {
703  fnSparseMCDplus[i] = new THnSparseF(Form("fnSparseAccDplus_%s",label[i].Data()),Form("MC nSparse D^{+} (Acc.Step)- %s",label[i].Data()),
704  knVarForSparseAcc, nBinsAcc, xminAcc, xmaxAcc);
705  fnSparseMCDplus[i]->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
706  fnSparseMCDplus[i]->GetAxis(1)->SetTitle("y");
707  fOutput->Add(fnSparseMCDplus[i]);
708  }
709  }
710  for (Int_t i=2; i<4; i++) {
711  fnSparseMC[i] = new THnSparseF(Form("fnSparseReco_%s",label[i-2].Data()),Form("MC nSparse (Reco Step)- %s",label[i-2].Data()),
712  knVarForSparse, nBinsReco, xminReco, xmaxReco);
713  for (Int_t j=0; j<knVarForSparse; j++) {
714  fnSparseMC[i]->GetAxis(j)->SetTitle(Form("%s",axis[j].Data()));
715  }
716  fOutput->Add(fnSparseMC[i]);
717 
718  //Dplus
719  if(fFillSparseDplus) {
720  fnSparseMCDplus[i] = new THnSparseF(Form("fnSparseRecoDplus_%s",label[i-2].Data()),Form("MC nSparse D^{+} (Reco Step)- %s",label[i-2].Data()),
721  knVarForSparse, nBinsReco, xminReco, xmaxReco);
722  for (Int_t j=0; j<knVarForSparse; j++) {
723  fnSparseMCDplus[i]->GetAxis(j)->SetTitle(Form("%s",axis[j].Data()));
724  }
725  fOutput->Add(fnSparseMCDplus[i]);
726  }
727  }
728 
729  fnSparseIP = new THnSparseF("fnSparseIP","nSparseIP", knVarForSparseIP, nBinsIP, xminIP, xmaxIP);
730  for (Int_t j=0; j<knVarForSparseIP; j++) {
731  fnSparseIP->GetAxis(j)->SetTitle(Form("%s",axisIP[j].Data()));
732  }
733  fnSparseIP->GetAxis(5)->SetTitle("candType (0.5=bkg; 1.5=prompt; 2.5=FD)");
734  fOutput->Add(fnSparseIP);
735  }
736  else {
737  fnSparse = new THnSparseF("fnSparse","nSparse", knVarForSparse, nBinsReco, xminReco, xmaxReco);
738  for (Int_t j=0; j<knVarForSparse; j++) {
739  fnSparse->GetAxis(j)->SetTitle(Form("%s",axis[j].Data()));
740  }
741  fOutput->Add(fnSparse);
742  }
743  }
744 
745  if(fFillImpParSparse) {
746  Int_t nBinsImpPar[3] = { 20, 200, 350};
747  Double_t xminImpPar[3] = { 0., 0., 1.6};
748  Double_t xmaxImpPar[3] = {20., 1000., 2.3};
749  TString axisImpPar[3] = {"Pt","imp.par. (#mum)","invMassDsAllPhi"};
750  if(!fReadMC) {
751  fImpParSparse = new THnSparseF("fImpParSparse","ImpParSparse", 3, nBinsImpPar, xminImpPar, xmaxImpPar);
752  for (Int_t j=0; j<3; j++) {
753  fImpParSparse->GetAxis(j)->SetTitle(Form("%s",axisImpPar[j].Data()));
754  }
755  fOutput->Add(fImpParSparse);
756  }
757  else {
758  nBinsImpPar[2] = 3;
759  xminImpPar[2] = 0.;
760  xmaxImpPar[2] = 3.;
761  axisImpPar[2] = "candType (0.5=bkg; 1.5=prompt; 2.5=FD)";
762  fImpParSparseMC = new THnSparseF("fImpParSparseMC","ImpParSparseMC", 3, nBinsImpPar, xminImpPar, xmaxImpPar);
763  for (Int_t j=0; j<3; j++) {
764  fImpParSparseMC->GetAxis(j)->SetTitle(Form("%s",axisImpPar[j].Data()));
765  }
766  fOutput->Add(fImpParSparseMC);
767  }
768  }
769 
770  //Counter for Normalization
771  fCounter = new AliNormalizationCounter("NormalizationCounter");
772  fCounter->Init();
773 
774  PostData(1,fOutput);
775  PostData(3,fCounter);
776 
777  if(fFillNtuple>0){
778  OpenFile(4); // 4 is the slot number of the ntuple
779 
780  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");
781 
782  }
783 
784 
785  return;
786 }
787 
788 //________________________________________________________________________
790 {
793 
794  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
795 
796  fHistNEvents->Fill(0); // all events
797  if(fAODProtection>=0){
798  // Protection against different number of events in the AOD and deltaAOD
799  // In case of discrepancy the event is rejected.
800  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
801  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
802  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
803  fHistNEvents->Fill(2);
804  return;
805  }
806  fHistNEvents->Fill(1);
807  }
808 
809  TClonesArray *array3Prong = 0;
810  if(!aod && AODEvent() && IsStandardAOD()) {
811  // In case there is an AOD handler writing a standard AOD, use the AOD
812  // event in memory rather than the input (ESD) event.
813  aod = dynamic_cast<AliAODEvent*> (AODEvent());
814  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
815  // have to taken from the AOD event hold by the AliAODExtension
816  AliAODHandler* aodHandler = (AliAODHandler*)
817  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
818  if(aodHandler->GetExtensions()) {
819  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
820  AliAODEvent *aodFromExt = ext->GetAOD();
821  array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
822  }
823  } else if(aod) {
824  array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
825  }
826 
827  if(!aod || !array3Prong) {
828  printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
829  return;
830  }
831 
832 
833  // fix for temporary bug in ESDfilter
834  // the AODs with null vertex pointer didn't pass the PhysSel
835  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
836 
837 
838  fHistNEvents->Fill(3); // count event
839  // Post the data already here
840  PostData(1,fOutput);
841 
843 
844 
845  Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
846  Float_t ntracks=aod->GetNumberOfTracks();
847  Float_t evCentr=fAnalysisCuts->GetCentrality(aod);
848 
849  fHistCentrality[0]->Fill(evCentr);
850  fHistCentralityMult[0]->Fill(ntracks,evCentr);
857  fHistNEvents->Fill(10);
858  fHistCentrality[2]->Fill(evCentr);
859  fHistCentralityMult[2]->Fill(ntracks,evCentr);
860  }
861 
863  Int_t runNumber=aod->GetRunNumber();
864 
865 
866 
867  TClonesArray *arrayMC=0;
868  AliAODMCHeader *mcHeader=0;
869 
870  // AOD primary vertex
871  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
872  // vtx1->Print();
873 
874  // load MC particles
875  if(fReadMC){
876 
877  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
878  if(!arrayMC) {
879  printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
880  return;
881  }
882 
883  // load MC header
884  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
885  if(!mcHeader) {
886  printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
887  return;
888  }
889  }
890 
891 
892  if(fReadMC && fFillSparse){
893  if(aod->GetTriggerMask()==0 && (runNumber>=195344 && runNumber<=195677))
894  // protection for events with empty trigger mask in p-Pb
895  return;
897  // events not passing the centrality selection can be removed immediately.
898  return;
899  Double_t zMCVertex = mcHeader->GetVtxZ();
900  if (TMath::Abs(zMCVertex) > fAnalysisCuts->GetMaxVtxZ())
901  return;
902  FillMCGenAccHistos(arrayMC, mcHeader);
903  }
904 
905  if(!isEvSel)return;
906  fHistNEvents->Fill(4);
907  fHistCentrality[1]->Fill(evCentr);
908  fHistCentralityMult[1]->Fill(ntracks,evCentr);
909 
910  if(fDoCutV0multTPCout) {
911  //cut on V0mult. vs #tracks kTPCout
912  Float_t V0mult=0.;
913  AliAODVZERO *aodVZERO = (AliAODVZERO*)aod->GetVZEROData();
914  if(aodVZERO) {
915  for(int ich=0; ich<64;ich++) V0mult += aodVZERO->GetMultiplicity(ich);
916  }
917  Int_t nTPCout=0;
918  for (Int_t i=0;i<aod->GetNumberOfTracks();++i) {
919  AliVTrack *track = aod->GetTrack(i);
920  if (!track) continue;
921  if((track->GetStatus() & AliVTrack::kTPCout)) nTPCout++;
922  }
923  fHistAllV0multNTPCout->Fill(V0mult,nTPCout);
924  if(nTPCout > (0.32*V0mult+750)) return;
925  else fHistSelV0multNTPCout->Fill(V0mult,nTPCout);
926  }
927 
928  Int_t n3Prong = array3Prong->GetEntriesFast();
929  if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
930 
931 
932  Int_t pdgDstoKKpi[3]={321,321,211};
933  Int_t nSelected=0;
934  Int_t nFiltered=0;
935  Double_t massPhi=TDatabasePDG::Instance()->GetParticle(333)->Mass();
936 
937  // vHF object is needed to call the method that refills the missing info of the candidates
938  // if they have been deleted in dAOD reconstruction phase
939  // in order to reduce the size of the file
941 
942  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
943 
944 
945  AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
946  fHistNEvents->Fill(11);
947 
949  continue;
950  }
951  nFiltered++;
952  fHistNEvents->Fill(12);
953 
954  if(!(vHF->FillRecoCand(aod,d))) {
955  fHistNEvents->Fill(14); //monitor how often this fails
956  continue;
957  }
958 
959  Bool_t unsetvtx=kFALSE;
960  if(!d->GetOwnPrimaryVtx()){
961  d->SetOwnPrimaryVtx(vtx1);
962  unsetvtx=kTRUE;
963  // NOTE: the ow primary vertex should be unset, otherwise there is a memory leak
964  // Pay attention if you use continue inside this loop!!!
965  }
966 
967  Bool_t recVtx=kFALSE;
968  AliAODVertex *origownvtx=0x0;
969 
970  Double_t ptCand = d->Pt();
971  Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
972  Double_t rapid=d->YDs();
973  fYVsPt->Fill(ptCand,rapid);
974  Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);
975 
976  if(isFidAcc){
977 
978  Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
979  Int_t retCodeNoRes=retCodeAnalysisCuts;
981  if(origRes){
983  retCodeNoRes=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
985  }
986  Double_t massKK = 0.;
987 
988  if(retCodeNoRes&1){ //KKpi
989  massKK=d->InvMass2Prongs(0,1,321,321);
990  Double_t massKp=d->InvMass2Prongs(1,2,321,211);
991  fMassHistKK[iPtBin]->Fill(massKK);
992  fMassHistKpi[iPtBin]->Fill(massKp);
993  }
994  if(retCodeNoRes&2){ //piKK
995  massKK=d->InvMass2Prongs(1,2,321,321);
996  Double_t massKp=d->InvMass2Prongs(0,1,211,321);
997  fMassHistKK[iPtBin]->Fill(massKK);
998  fMassHistKpi[iPtBin]->Fill(massKp);
999  }
1000 
1001  Int_t isKKpi=retCodeAnalysisCuts&1;
1002  Int_t ispiKK=retCodeAnalysisCuts&2;
1003  Int_t isPhiKKpi=retCodeAnalysisCuts&4;
1004  Int_t isPhipiKK=retCodeAnalysisCuts&8;
1005  Int_t isK0starKKpi=retCodeAnalysisCuts&16;
1006  Int_t isK0starpiKK=retCodeAnalysisCuts&32;
1007 
1008  if(retCodeAnalysisCuts>0){
1010  if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
1011  if(fAnalysisCuts->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
1012  else fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
1013  }
1014 
1015 
1016  fHistNEvents->Fill(13);
1017  nSelected++;
1018 
1019  Int_t index=GetHistoIndex(iPtBin);
1020  fPtCandHist[index]->Fill(ptCand);
1021 
1022 
1023  Double_t weightKKpi=1.;
1024  Double_t weightpiKK=1.;
1026  weightKKpi=fAnalysisCuts->GetWeightForKKpi();
1027  weightpiKK=fAnalysisCuts->GetWeightForpiKK();
1028  if(weightKKpi>1. || weightKKpi<0.) weightKKpi=0.;
1029  if(weightpiKK>1. || weightpiKK<0.) weightpiKK=0.;
1030  }
1031 
1032  fChanHist[0]->Fill(retCodeAnalysisCuts);
1033 
1034 
1035  Double_t invMass = 0.;
1036  Int_t indexMCKKpi=-1;
1037  Int_t indexMCpiKK=-1;
1038  Int_t labDs=-1;
1039  Int_t labDplus=-1;
1040  Int_t pdgCode0=-999;
1041  Int_t isMCSignal=-1;
1042 
1043 
1044  if(fReadMC){
1045 
1046  labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
1047  labDplus = d->MatchToMC(411,arrayMC,3,pdgDstoKKpi);
1048  if(labDs>=0){
1049  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
1050  AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
1051  pdgCode0=TMath::Abs(p->GetPdgCode());
1052 
1053  if(isKKpi){
1054  if(pdgCode0==321) {
1055  indexMCKKpi=GetSignalHistoIndex(iPtBin);
1056  fYVsPtSig->Fill(ptCand,rapid);
1057  fChanHist[1]->Fill(retCodeAnalysisCuts);
1058  isMCSignal=1;
1059  }else{
1060  indexMCKKpi=GetReflSignalHistoIndex(iPtBin);
1061  fChanHist[3]->Fill(retCodeAnalysisCuts);
1062  isMCSignal=0;
1063  }
1064  }
1065  if(ispiKK){
1066  if(pdgCode0==211) {
1067  indexMCpiKK=GetSignalHistoIndex(iPtBin);
1068  fYVsPtSig->Fill(ptCand,rapid);
1069  fChanHist[1]->Fill(retCodeAnalysisCuts);
1070  isMCSignal=1;
1071  }else{
1072  indexMCpiKK=GetReflSignalHistoIndex(iPtBin);
1073  fChanHist[3]->Fill(retCodeAnalysisCuts);
1074  isMCSignal=0;
1075  }
1076  }
1077  }else{
1078  indexMCpiKK=GetBackgroundHistoIndex(iPtBin);
1079  indexMCKKpi=GetBackgroundHistoIndex(iPtBin);
1080  fChanHist[2]->Fill(retCodeAnalysisCuts);
1081  if(labDplus>=0) {
1082  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
1083  AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
1084  pdgCode0=TMath::Abs(p->GetPdgCode());
1085  }
1086  }
1087  }
1088 
1089  Double_t candType = 0.5; //for bkg
1090  Double_t massPhi=TDatabasePDG::Instance()->GetParticle(333)->Mass();
1091 
1092  if(isKKpi){
1093  if(fDoRotBkg && TMath::Abs(massKK-massPhi)<=fMaxDeltaPhiMass4Rot)GenerateRotBkg(d,1,iPtBin);
1094 
1095  invMass=d->InvMassDsKKpi();
1096  fMassHist[index]->Fill(invMass,weightKKpi);
1097  fPtVsMass->Fill(invMass,ptCand,weightKKpi);
1098 
1099  if(fDoBkgPhiSB && 0.010<TMath::Abs(massKK-massPhi)<0.030) {
1100  if(massKK<massPhi)fMassLSBkgHistPhi[iPtBin]->Fill(invMass);
1101  else fMassRSBkgHistPhi[iPtBin]->Fill(invMass);
1102  }
1103 
1104  if(isPhiKKpi){
1105  fMassHistPhi[index]->Fill(invMass,weightKKpi);
1106  fPtVsMassPhi->Fill(invMass,ptCand,weightKKpi);
1107  }
1108  if(isK0starKKpi){
1109  fMassHistK0st[index]->Fill(invMass,weightKKpi);
1110  fPtVsMassK0st->Fill(invMass,ptCand,weightKKpi);
1111  }
1112  if(fReadMC && indexMCKKpi!=-1){
1113  fMassHist[indexMCKKpi]->Fill(invMass,weightKKpi);
1114  if(isPhiKKpi) {
1115  fMassHistPhi[indexMCKKpi]->Fill(invMass,weightKKpi);
1116  if(fFillSparse) {
1117  if(indexMCKKpi==GetSignalHistoIndex(iPtBin) || labDplus >= 0) {
1118  AliAODMCParticle *partDs;
1119  if(indexMCKKpi==GetSignalHistoIndex(iPtBin)) partDs = (AliAODMCParticle*)arrayMC->At(labDs);
1120  if(labDplus >= 0) partDs = (AliAODMCParticle*)arrayMC->At(labDplus);
1121  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partDs,kTRUE);
1122  if(orig==4) {
1123  candType = 1.5;
1124  }
1125  if(orig==5) {
1126  candType = 2.5;
1127  }
1128  }
1129  }
1130  }
1131  if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass,weightKKpi);
1132  }
1133  if(isPhiKKpi && fFillImpParSparse) {
1134  Double_t impParxy = d->ImpParXY()*10000.;
1135  if(!fReadMC) {
1136  Double_t array4ImpPar[3] = {ptCand,impParxy,invMass};
1137  fImpParSparse->Fill(array4ImpPar);
1138  }
1139  else {
1140  Double_t array4ImpPar[3] = {ptCand,impParxy,candType};
1141  fImpParSparseMC->Fill(array4ImpPar);
1142  }
1143  }
1144  }
1145  if(ispiKK){
1146  if(fDoRotBkg && TMath::Abs(massKK-massPhi)<=fMaxDeltaPhiMass4Rot)GenerateRotBkg(d,2,iPtBin);
1147 
1148  invMass=d->InvMassDspiKK();
1149  fMassHist[index]->Fill(invMass,weightpiKK);
1150  fPtVsMass->Fill(invMass,ptCand,weightpiKK);
1151 
1152  if(fDoBkgPhiSB && 0.010<TMath::Abs(massKK-massPhi)<0.030) {
1153  if(massKK<massPhi)fMassLSBkgHistPhi[iPtBin]->Fill(invMass);
1154  else fMassRSBkgHistPhi[iPtBin]->Fill(invMass);
1155  }
1156 
1157  if(isPhipiKK){
1158  fMassHistPhi[index]->Fill(invMass,weightpiKK);
1159  fPtVsMassPhi->Fill(invMass,ptCand,weightpiKK);
1160  }
1161  if(isK0starpiKK){
1162  fMassHistK0st[index]->Fill(invMass,weightpiKK);
1163  fPtVsMassK0st->Fill(invMass,ptCand,weightpiKK);
1164  }
1165  if(fReadMC && indexMCpiKK!=-1){
1166  fMassHist[indexMCpiKK]->Fill(invMass,weightpiKK);
1167  if(isPhipiKK) {
1168  fMassHistPhi[indexMCpiKK]->Fill(invMass,weightpiKK);
1169  if(fFillSparse) {
1170  if(indexMCpiKK==GetSignalHistoIndex(iPtBin) || labDplus >= 0) {
1171  AliAODMCParticle *partDs;
1172  if(indexMCpiKK==GetSignalHistoIndex(iPtBin)) partDs = (AliAODMCParticle*)arrayMC->At(labDs);
1173  if(labDplus >= 0) partDs = (AliAODMCParticle*)arrayMC->At(labDplus);
1174  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partDs,kTRUE);
1175  if(orig==4) {
1176  candType = 1.5;
1177  }
1178  if(orig==5) {
1179  candType = 2.5;
1180  }
1181  }
1182  }
1183  }
1184  if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass,weightpiKK);
1185  }
1186  if(isPhipiKK && fFillImpParSparse) {
1187  Double_t impParxy = d->ImpParXY()*10000.;
1188  if(!fReadMC) {
1189  Double_t array4ImpPar[3] = {ptCand,impParxy,invMass};
1190  fImpParSparse->Fill(array4ImpPar);
1191  }
1192  else {
1193  Double_t array4ImpPar[3] = {ptCand,impParxy,candType};
1194  fImpParSparseMC->Fill(array4ImpPar);
1195  }
1196  }
1197  }
1198 
1199 
1201 
1202  const Int_t nProng = 3;
1203  Double_t deltaMassKK=999.;
1204  Double_t dlen=d->DecayLength();
1205  Double_t dlenxy=d->DecayLengthXY();
1206  Double_t normdl=d->NormalizedDecayLength();
1207  Double_t normdlxy=d->NormalizedDecayLengthXY();
1208  Double_t cosp=d->CosPointingAngle();
1209  Double_t cospxy=d->CosPointingAngleXY();
1210  Double_t pt0=d->PtProng(0);
1211  Double_t pt1=d->PtProng(1);
1212  Double_t pt2=d->PtProng(2);
1213  Double_t sigvert=d->GetSigmaVert();
1214  Double_t cosPiDs=-99.;
1215  Double_t cosPiKPhi=-99.;
1216  Double_t normIP; //to store the maximum topomatic var. among the 3 prongs
1217  Double_t normIPprong[nProng]; //to store IP of k,k,pi
1218 
1219  if(fFillSparse) {
1220 
1221  if(isPhiKKpi) {
1222  Double_t tmpNormIP[nProng];
1223 
1224  invMass = d->InvMassDsKKpi();
1225  massKK = d->InvMass2Prongs(0,1,321,321);
1226  deltaMassKK = TMath::Abs(massKK-massPhi);
1227  cosPiDs = d->CosPiDsLabFrameKKpi();
1228  cosPiKPhi = d->CosPiKPhiRFrameKKpi();
1229  cosPiKPhi = TMath::Abs(cosPiKPhi*cosPiKPhi*cosPiKPhi);
1230 
1231  for(Int_t ip=0; ip<nProng; ip++) {
1232  Double_t diffIP, errdiffIP;
1233  d->Getd0MeasMinusExpProng(ip,aod->GetMagneticField(),diffIP,errdiffIP);
1234  tmpNormIP[ip] = diffIP/errdiffIP;
1235  if(ip==0) normIP = tmpNormIP[ip];
1236  else if(TMath::Abs(tmpNormIP[ip])>TMath::Abs(normIP)) normIP = tmpNormIP[ip];
1237  }
1238  normIPprong[0] = tmpNormIP[0];
1239  normIPprong[1] = tmpNormIP[1];
1240  normIPprong[2] = tmpNormIP[2];
1241 
1242  Double_t var4nSparse[knVarForSparse] = {invMass,ptCand,deltaMassKK*1000,dlen*1000,dlenxy*1000,normdlxy,cosp*100,cospxy*100,
1243  sigvert*1000,cosPiDs*10,cosPiKPhi*10,TMath::Abs(normIP)};
1244 
1245  if(!fReadMC) {
1246  fnSparse->Fill(var4nSparse);
1247  }
1248  else {
1249  if(indexMCKKpi==GetSignalHistoIndex(iPtBin)) {
1250  if(candType==1.5) fnSparseMC[2]->Fill(var4nSparse);
1251  if(candType==2.5) fnSparseMC[3]->Fill(var4nSparse);
1252  }
1253  else if(fFillSparseDplus && labDplus>=0 && pdgCode0==321) {
1254  if(candType==1.5) fnSparseMCDplus[2]->Fill(var4nSparse);
1255  if(candType==2.5) fnSparseMCDplus[3]->Fill(var4nSparse);
1256  }
1257  }
1258  }
1259  if(isPhipiKK) {
1260  Double_t tmpNormIP[nProng];
1261 
1262  invMass = d->InvMassDspiKK();
1263  massKK = d->InvMass2Prongs(1,2,321,321);
1264  deltaMassKK = TMath::Abs(massKK-massPhi);
1265  cosPiDs = d->CosPiDsLabFramepiKK();
1266  cosPiKPhi = d->CosPiKPhiRFramepiKK();
1267  cosPiKPhi = TMath::Abs(cosPiKPhi*cosPiKPhi*cosPiKPhi);
1268 
1269  for(Int_t ip=0; ip<nProng; ip++) {
1270  Double_t diffIP, errdiffIP;
1271  d->Getd0MeasMinusExpProng(ip,aod->GetMagneticField(),diffIP,errdiffIP);
1272  tmpNormIP[ip] = diffIP/errdiffIP;
1273  if(ip==0) normIP = tmpNormIP[ip];
1274  else if(TMath::Abs(tmpNormIP[ip])>TMath::Abs(normIP)) normIP = tmpNormIP[ip];
1275  }
1276 
1277  normIPprong[0] = tmpNormIP[2];
1278  normIPprong[1] = tmpNormIP[1];
1279  normIPprong[2] = tmpNormIP[0];
1280 
1281  Double_t var4nSparse[knVarForSparse] = {invMass,ptCand,deltaMassKK*1000,dlen*1000,dlenxy*1000,normdlxy,cosp*100,cospxy*100,
1282  sigvert*1000,cosPiDs*10,cosPiKPhi*10,TMath::Abs(normIP)};
1283 
1284  if(!fReadMC) {
1285  fnSparse->Fill(var4nSparse);
1286  }
1287  else {
1288  if(indexMCpiKK==GetSignalHistoIndex(iPtBin)) {
1289  if(candType==1.5) fnSparseMC[2]->Fill(var4nSparse);
1290  if(candType==2.5) fnSparseMC[3]->Fill(var4nSparse);
1291  }
1292  else if(fFillSparseDplus && labDplus>=0 && pdgCode0==211) {
1293  if(candType==1.5) fnSparseMCDplus[2]->Fill(var4nSparse);
1294  if(candType==2.5) fnSparseMCDplus[3]->Fill(var4nSparse);
1295  }
1296  }
1297  }
1298 
1299  if(fReadMC && (isPhiKKpi || isPhiKKpi)) {
1300  Double_t var[6] = {ptCand,normIP,normIPprong[0],normIPprong[1],normIPprong[2],candType};
1301  fnSparseIP->Fill(var);
1302  }
1303  }
1305 
1306  if(fDoCutVarHistos){
1307  Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
1308  Double_t dca=d->GetDCA();
1309  Double_t ptmax=0;
1310 
1311  for(Int_t i=0;i<3;i++){
1312  if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
1313  }
1314  fCosPHist[index]->Fill(cosp);
1315  fCosPxyHist[index]->Fill(cospxy);
1316  fDLenHist[index]->Fill(dlen);
1317  fDLenxyHist[index]->Fill(dlenxy);
1318  fNDLenxyHist[index]->Fill(normdlxy);
1319  fSigVertHist[index]->Fill(sigvert);
1320  fSumd02Hist[index]->Fill(sumD02);
1321  fPtMaxHist[index]->Fill(ptmax);
1322  fDCAHist[index]->Fill(dca);
1323  fNormIPHist[index]->Fill(normIP);
1324  fCosPiDsHist[index]->Fill(cosPiDs);
1325  fCosPiKPhiHist[index]->Fill(cosPiKPhi);
1326  fPtProng0Hist[index]->Fill(pt0);
1327  fPtProng1Hist[index]->Fill(pt1);
1328  fPtProng2Hist[index]->Fill(pt2);
1329  if(!fReadMC) {
1330  fCosPHist3D->Fill(invMass,ptCand,cosp);
1331  fCosPxyHist3D->Fill(invMass,ptCand,cospxy);
1332  fDLenHist3D->Fill(invMass,ptCand,dlen);
1333  fDLenxyHist3D->Fill(invMass,ptCand,dlenxy);
1334  fNDLenxyHist3D->Fill(invMass,ptCand,normdlxy);
1335  fSigVertHist3D->Fill(invMass,ptCand,sigvert);
1336  fDCAHist3D->Fill(invMass,ptCand,dca);
1337  fNormIPHist3D->Fill(invMass,ptCand,normIP);
1338  fCosPiDsHist3D->Fill(invMass,ptCand,cosPiDs);
1339  fCosPiKPhiHist3D->Fill(invMass,ptCand,cosPiKPhi);
1340  fPtProng0Hist3D->Fill(invMass,ptCand,pt0);
1341  fPtProng1Hist3D->Fill(invMass,ptCand,pt1);
1342  fPtProng2Hist3D->Fill(invMass,ptCand,pt2);
1343  }
1344  if(isKKpi){
1345  Double_t massKK=d->InvMass2Prongs(0,1,321,321);
1346  Double_t massKp=d->InvMass2Prongs(1,2,321,211);
1347  fDalitz[index]->Fill(massKK,massKp);
1348  if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);
1349  if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);
1350  if(fReadMC && indexMCKKpi!=-1){
1351  fDalitz[indexMCKKpi]->Fill(massKK,massKp);
1352  if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);
1353  if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);
1354  fCosPHist[indexMCKKpi]->Fill(cosp);
1355  fCosPxyHist[indexMCKKpi]->Fill(cospxy);
1356  fDLenHist[indexMCKKpi]->Fill(dlen);
1357  fDLenxyHist[indexMCKKpi]->Fill(dlenxy);
1358  fNDLenxyHist[indexMCKKpi]->Fill(normdlxy);
1359  fSigVertHist[indexMCKKpi]->Fill(sigvert);
1360  fSumd02Hist[indexMCKKpi]->Fill(sumD02);
1361  fPtMaxHist[indexMCKKpi]->Fill(ptmax);
1362  fPtCandHist[indexMCKKpi]->Fill(ptCand);
1363  fDCAHist[indexMCKKpi]->Fill(dca);
1364  fNormIPHist[indexMCKKpi]->Fill(normIP);
1365  fCosPiDsHist[indexMCKKpi]->Fill(cosPiDs);
1366  fCosPiKPhiHist[indexMCKKpi]->Fill(cosPiKPhi);
1367  fPtProng0Hist[indexMCKKpi]->Fill(pt0);
1368  fPtProng1Hist[indexMCKKpi]->Fill(pt1);
1369  fPtProng2Hist[indexMCKKpi]->Fill(pt2);
1370  }
1371  }
1372  if(ispiKK){
1373  Double_t massKK=d->InvMass2Prongs(1,2,321,321);
1374  Double_t massKp=d->InvMass2Prongs(0,1,211,321);
1375  fDalitz[index]->Fill(massKK,massKp);
1376  if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);
1377  if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);
1378 
1379 
1380  if(fReadMC && indexMCpiKK!=-1){
1381  fDalitz[indexMCpiKK]->Fill(massKK,massKp);
1382  if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);
1383  if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);
1384  fCosPHist[indexMCpiKK]->Fill(cosp);
1385  fCosPxyHist[indexMCpiKK]->Fill(cospxy);
1386  fDLenHist[indexMCpiKK]->Fill(dlen);
1387  fDLenxyHist[indexMCpiKK]->Fill(dlenxy);
1388  fNDLenxyHist[indexMCpiKK]->Fill(normdlxy);
1389  fSigVertHist[indexMCpiKK]->Fill(sigvert);
1390  fSumd02Hist[indexMCpiKK]->Fill(sumD02);
1391  fPtMaxHist[indexMCpiKK]->Fill(ptmax);
1392  fPtCandHist[indexMCpiKK]->Fill(ptCand);
1393  fDCAHist[indexMCpiKK]->Fill(dca);
1394  fNormIPHist[indexMCpiKK]->Fill(normIP);
1395  fCosPiDsHist[indexMCpiKK]->Fill(cosPiDs);
1396  fCosPiKPhiHist[indexMCpiKK]->Fill(cosPiKPhi);
1397  fPtProng0Hist[indexMCpiKK]->Fill(pt0);
1398  fPtProng1Hist[indexMCpiKK]->Fill(pt1);
1399  fPtProng2Hist[indexMCpiKK]->Fill(pt2);
1400  }
1401  }
1402  }
1403 
1404  Float_t tmp[37];
1405  if(fFillNtuple>0){
1406 
1407  if ((fFillNtuple==1 && (isPhiKKpi || isPhipiKK)) || (fFillNtuple==2 && (isK0starKKpi || isK0starpiKK)) || (fFillNtuple==3 && (isKKpi || ispiKK))){
1408 
1409  AliAODTrack *track0=(AliAODTrack*)d->GetDaughter(0);
1410  AliAODTrack *track1=(AliAODTrack*)d->GetDaughter(1);
1411  AliAODTrack *track2=(AliAODTrack*)d->GetDaughter(2);
1412 
1413  UInt_t bitMapPIDTrack0=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track0);
1414  UInt_t bitMapPIDTrack1=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track1);
1415  UInt_t bitMapPIDTrack2=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track2);
1416 
1417  tmp[0]=Float_t(labDs);
1418  if(fReadMC && fWriteOnlySignal) tmp[0]=Float_t(isMCSignal);
1419  tmp[1]=Float_t(retCodeAnalysisCuts);
1420  tmp[2]=Float_t(pdgCode0);
1421  tmp[3]=d->PtProng(0);
1422  tmp[4]=d->PtProng(1);
1423  tmp[5]=d->PtProng(2);
1424  tmp[6]=d->Pt();
1425  tmp[7]=d->PProng(0);
1426  tmp[8]=d->PProng(1);
1427  tmp[9]=d->PProng(2);
1428  tmp[10]=Int_t(bitMapPIDTrack0);
1429  tmp[11]=Int_t(bitMapPIDTrack1);
1430  tmp[12]=Int_t(bitMapPIDTrack2);
1431  tmp[13]=d->CosPointingAngle();
1432  tmp[14]=d->CosPointingAngleXY();
1433  tmp[15]=d->DecayLength();
1434  tmp[16]=d->DecayLengthXY();
1435  tmp[17]=d->NormalizedDecayLength();
1436  tmp[18]=d->NormalizedDecayLengthXY();
1437  tmp[19]=d->InvMassDsKKpi();
1438  tmp[20]=d->InvMassDspiKK();
1439  tmp[21]=d->GetSigmaVert();
1440  tmp[22]=d->Getd0Prong(0);
1441  tmp[23]=d->Getd0Prong(1);
1442  tmp[24]=d->Getd0Prong(2);
1443  tmp[25]=d->GetDCA();
1444  tmp[26]=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
1445  tmp[27]=d->InvMass2Prongs(0,1,321,321);
1446  tmp[28]=d->InvMass2Prongs(1,2,321,321);
1447  tmp[29]=d->InvMass2Prongs(1,2,321,211);
1448  tmp[30]=d->InvMass2Prongs(0,1,211,321);
1449  tmp[31]=d->CosPiDsLabFrameKKpi();
1450  tmp[32]=d->CosPiDsLabFramepiKK();
1451  tmp[33]=d->CosPiKPhiRFrameKKpi();
1452  tmp[34]=d->CosPiKPhiRFramepiKK();
1453  tmp[35]=(Float_t)(centrality);
1454  tmp[36]=(Float_t)(runNumber);
1455 
1456  if(fReadMC && fWriteOnlySignal){
1457  if(isMCSignal>=0) fNtupleDs->Fill(tmp);
1458  }else{
1459  fNtupleDs->Fill(tmp);
1460  }
1461  PostData(4,fNtupleDs);
1462  }
1463  }
1464  } //if(retCodeAnalysisCuts
1465  } // if(isFidAcc)
1466 
1467  if(unsetvtx) d->UnsetOwnPrimaryVtx();
1468  if(recVtx)fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
1469  }
1470 
1471  fCounter->StoreCandidates(aod,nFiltered,kTRUE);
1472  fCounter->StoreCandidates(aod,nSelected,kFALSE);
1473 
1474  delete vHF;
1475 
1476  PostData(1,fOutput);
1477  PostData(3,fCounter);
1478 
1479  return;
1480 }
1481 
1482 //_________________________________________________________________
1483 
1485 {
1487  //
1488  if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
1489  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1490  if (!fOutput) {
1491  printf("ERROR: fOutput not available\n");
1492  return;
1493  }
1494  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
1495  if(fHistNEvents){
1496  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1497  }else{
1498  printf("ERROR: fHistNEvents not available\n");
1499  return;
1500  }
1501  return;
1502 }
1503 
1504 //_________________________________________________________________
1505 void AliAnalysisTaskSEDs::FillMCGenAccHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader){
1507 
1508  Int_t nProng = 3;
1509  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
1510  if(TMath::Abs(zMCVertex) <= fAnalysisCuts->GetMaxVtxZ()) {
1511  for(Int_t iPart=0; iPart<arrayMC->GetEntriesFast(); iPart++){
1512 
1513  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
1514 
1515  if(TMath::Abs(mcPart->GetPdgCode()) == 431) {
1516  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,kTRUE);//Prompt = 4, FeedDown = 5
1517 
1518  Int_t deca = 0;
1519  Bool_t isGoodDecay = kFALSE;
1520  Int_t labDau[3] = {-1,-1,-1};
1521  Bool_t isFidAcc = kFALSE;
1522  Bool_t isDaugInAcc = kFALSE;
1523 
1524  deca = AliVertexingHFUtils::CheckDsDecay(arrayMC,mcPart,labDau);
1525  if(deca == 1) isGoodDecay=kTRUE; // == 1 -> Phi pi -> kkpi
1526 
1527  if(labDau[0]==-1) continue; //protection against unfilled array of labels
1528 
1529  if(isGoodDecay) {
1530  Double_t pt = mcPart->Pt();
1531  Double_t rapid = mcPart->Y();
1532  isFidAcc = fAnalysisCuts->IsInFiducialAcceptance(pt,rapid);
1533  isDaugInAcc = CheckDaugAcc(arrayMC,nProng,labDau);
1534 
1535  if(isFidAcc) {
1536  Double_t var4nSparseAcc[2] = {pt,rapid*10};
1537  if(isDaugInAcc) {
1538  if(orig==4) fnSparseMC[0]->Fill(var4nSparseAcc);
1539  if(orig==5) fnSparseMC[1]->Fill(var4nSparseAcc);
1540  }
1541  }
1542  }
1543  }
1544 
1545  if(fFillSparseDplus && TMath::Abs(mcPart->GetPdgCode()) == 411) {
1546  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,kTRUE);//Prompt = 4, FeedDown = 5
1547 
1548  Int_t deca = 0;
1549  Bool_t isGoodDecay = kFALSE;
1550  Int_t labDau[3] = {-1,-1,-1};
1551  Bool_t isFidAcc = kFALSE;
1552  Bool_t isDaugInAcc = kFALSE;
1553 
1554  deca = AliVertexingHFUtils::CheckDplusKKpiDecay(arrayMC,mcPart,labDau);
1555  if(deca == 1) isGoodDecay=kTRUE; // == 1 -> Phi pi -> kkpi
1556 
1557  if(labDau[0]==-1) continue; //protection against unfilled array of labels
1558 
1559  if(isGoodDecay) {
1560  Double_t pt = mcPart->Pt();
1561  Double_t rapid = mcPart->Y();
1562  isFidAcc = fAnalysisCuts->IsInFiducialAcceptance(pt,rapid);
1563  isDaugInAcc = CheckDaugAcc(arrayMC,nProng,labDau);
1564 
1565  if(isFidAcc) {
1566  Double_t var4nSparseAcc[2] = {pt,rapid*10};
1567  if(isDaugInAcc) {
1568  if(orig==4) fnSparseMCDplus[0]->Fill(var4nSparseAcc);
1569  if(orig==5) fnSparseMCDplus[1]->Fill(var4nSparseAcc);
1570  }
1571  }
1572  }
1573  }
1574  }
1575  }
1576 }
1577 //_________________________________________________________________
1578 Bool_t AliAnalysisTaskSEDs::CheckDaugAcc(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1580 
1581  for (Int_t iProng = 0; iProng<nProng; iProng++){
1582  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1583  if(!mcPartDaughter) {
1584  return kFALSE;
1585  }
1586  Double_t eta = mcPartDaughter->Eta();
1587  Double_t pt = mcPartDaughter->Pt();
1588  if (TMath::Abs(eta) > 0.9 || pt < 0.1) {
1589  return kFALSE;
1590  }
1591  }
1592  return kTRUE;
1593 }
1594 
1595 //_________________________________________________________________
1596 
1598 
1599  const Int_t nprongs = 3;
1600  Double_t PxProng[nprongs], PyProng[nprongs], PtProng[nprongs], PzProng[nprongs], P2Prong[nprongs], mProng[nprongs];
1601  Double_t Px, Py, Pz, P2;
1602  UInt_t pdg[3]={321,321,211};
1603  int idPion = 2;
1604  if(dec==2) {
1605  pdg[0]=211;
1606  pdg[2]=321;
1607  idPion = 0;
1608  }
1609 
1610  for (Int_t ip=0; ip<nprongs; ip++) {
1611  PxProng[ip] = d->PxProng(ip);
1612  PyProng[ip] = d->PxProng(ip);
1613  PtProng[ip] = d->PtProng(ip);
1614  PzProng[ip] = d->PzProng(ip);
1615  P2Prong[ip] = d->P2Prong(ip);
1616  mProng[ip] = TDatabasePDG::Instance()->GetParticle(pdg[ip])->Mass();
1617  }
1618 
1619  for(Int_t i=0; i<9; i++) { //9 rotations implemented for the pion track around Pi
1620  Px = 0.;
1621  Py = 0.;
1622  Pz = 0.;
1623 
1624  Double_t phirot=TMath::Pi()*(5/6.+1/27.*i);
1625 
1626  PxProng[idPion] = PxProng[idPion]*TMath::Cos(phirot)-PyProng[idPion]*TMath::Sin(phirot);
1627  PyProng[idPion] = PxProng[idPion]*TMath::Sin(phirot)+PyProng[idPion]*TMath::Cos(phirot);
1628 
1629  for (Int_t j=0; j<nprongs; j++) {
1630  Px += PxProng[j];
1631  Py += PyProng[j];
1632  Pz += PzProng[j];
1633  }
1634  P2 = Px*Px + Py*Py + Pz*Pz;
1635 
1636  Double_t energysum = 0.;
1637  for(Int_t j=0; j<nprongs; j++) {
1638  energysum += TMath::Sqrt(mProng[j]*mProng[j]+P2Prong[j]);
1639  }
1640  Double_t mass = TMath::Sqrt(energysum*energysum-P2);
1641  if(fminMass<=mass<fmaxMass) fMassRotBkgHistPhi[iPtBin]->Fill(mass);
1642  }
1643 }
1644 
1645 
1646 
1647 
1648 
1649 
1650 
1651 
1652 
Double_t NormalizedDecayLengthXY() const
Bool_t IsEventRejectedDueToCentrality() const
Definition: AliRDHFCuts.h:321
Int_t pdg
THnSparseF * fnSparseMCDplus[4]
TH1F * fPtMaxHist[4 *kMaxPtBins]
! hist. for Pt Max (Prod Cuts)
Bool_t IsEventRejectedDueToZVertexOutsideFiducialRegion() const
Definition: AliRDHFCuts.h:315
Double_t NormalizedDecayLength() const
THnSparseF * fnSparseMC[4]
!<!THnSparse for topomatic variable
Bool_t IsEventRejectedDueToNotRecoVertex() const
Definition: AliRDHFCuts.h:309
double Double_t
Definition: External.C:58
Definition: External.C:260
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)
TH3F * fCosPxyHist3D
! cosPxy vs Ds mass vs pt
TH1F * fDLenxyHist[4 *kMaxPtBins]
! hist. of norm decay length XY (sig,bkg,tot)
TH1F * fNormIPHist[4 *kMaxPtBins]
! hist. for topomatic variable
TH3F * fPtProng1Hist3D
! Pt prong1 vs Ds mass vs pt
void SetPtBins(Int_t n, Float_t *lim)
Int_t IsEventSelectedInCentrality(AliVEvent *event)
Bool_t HasSelectionBit(Int_t i) const
TH1F * fNDLenxyHist[4 *kMaxPtBins]
! hist. of decay length XY (sig,bkg,tot)
TH2F * fDalitzK0st[4 *kMaxPtBins]
! dalitz plot via K0* (sig,bkg,tot)
static Int_t CheckDsDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
Bool_t fFillSparseDplus
flag for usage of THnSparse
Double_t mass
Double_t CosPiDsLabFrameKKpi() const
Double_t ImpParXY() const
centrality
Int_t GetHistoIndex(Int_t iPtBin) const
static Int_t CheckMatchingAODdeltaAODevents()
THnSparseF * fImpParSparse
!<!THnSparse for MC for D+->kkpi
TH3F * fSigVertHist3D
! SigVert vs Ds mass vs pt
Double_t fMassRange
limits for pt bins
Bool_t IsEventRejectedDueToVertexContributors() const
Definition: AliRDHFCuts.h:312
TH2F * fPtVsMass
! hist. of pt vs. mass (prod. cuts)
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:247
TH1F * fSumd02Hist[4 *kMaxPtBins]
! hist. for sum d02 (Prod Cuts)
TH3F * fCosPiKPhiHist3D
! cosPiKPhi vs Ds mass vs pt
UInt_t GetPIDTrackTPCTOFBitMap(AliAODTrack *track) const
TH1F * fMassHistPhi[4 *kMaxPtBins]
! hist. of mass spectra via phi (sig,bkg,tot)
TH1F * fCosPiKPhiHist[4 *kMaxPtBins]
! hist. for CosPiKPhi
Bool_t fDoBkgPhiSB
flag to create rotational bkg (rotating pi track)
THnSparseF * fnSparse
Cuts for Analysis.
TNtuple * fNtupleDs
! output ntuple
THnSparseF * fImpParSparseMC
!<!THnSparse for imp. par. on data
TH1F * fMassRotBkgHistPhi[kMaxPtBins]
! hist. of bkg generated from rot. of the pion
void GenerateRotBkg(AliAODRecoDecayHF3Prong *d, Int_t dec, Int_t iPtBin)
Bool_t fDoRotBkg
flag for usage of sparse for imp. parameter
TH3F * fPtProng0Hist3D
! Pt prong0 vs Ds mass vs pt
TH1F * fHistCentrality[3]
!hist. for cent distr (all,sel ev, )
TH3F * fDCAHist3D
! DCA vs Ds mass vs pt
int Int_t
Definition: External.C:63
TH3F * fCosPHist3D
! cosP vs Ds mass vs pt
AliNormalizationCounter * fCounter
flag to set mass window of phi meson (when using pion rotation to create bkg)
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
static Int_t CheckDplusKKpiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
AliRDHFCutsDstoKKpi * fAnalysisCuts
TList * fOutput
! list send on output slot 0
TH2F * fHistAllV0multNTPCout
! histo for V0mult vs #tracks TPCout (all)
Double_t CosPiKPhiRFrameKKpi() const
TH1F * fCosPxyHist[4 *kMaxPtBins]
! hist. of cosXY pointing angle (sig,bkg,tot)
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
TH3F * fPtProng2Hist3D
! Pt prong2 vs Ds mass vs pt
Bool_t CheckDaugAcc(TClonesArray *arrayMC, Int_t nProng, Int_t *labDau)
TH2F * fDalitzPhi[4 *kMaxPtBins]
! dalitz plot via phi (sig,bkg,tot)
TH2F * fHistSelV0multNTPCout
! histo for V0mult vs #tracks TPCout (sel)
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:248
Bool_t IsEventRejectedDueToPileup() const
Definition: AliRDHFCuts.h:318
TH3F * fDLenHist3D
! Dlen vs Ds mass vs pt
TH1F * fMassHistKpi[kMaxPtBins]
! hist. of mass spectra of Kpi
TH1F * fMassRSBkgHistPhi[kMaxPtBins]
! hist. of bkg generated from right phi sideband + pion
TH3F * fNormIPHist3D
! nIP vs Ds mass vs pt
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 to activate cut on V0mult vs #tracks TPCout
Bool_t GetIsPrimaryWithoutDaughters() const
Definition: AliRDHFCuts.h:268
TH1F * fCosPiDsHist[4 *kMaxPtBins]
! hist. for CosPiDs
Int_t GetSignalHistoIndex(Int_t iPtBin) const
Bool_t IsEventSelected(AliVEvent *event)
TH1F * fMassLSBkgHistPhi[kMaxPtBins]
! hist. of bkg generated from left phi sideband + pion
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
TH3F * fDLenxyHist3D
! Dlenxy vs Ds mass vs pt
Bool_t fFillImpParSparse
flag for usage of THnSparse
Double_t minMass
Bool_t fDoCutV0multTPCout
flag to create bkg from phi sidebands
void CleanOwnPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod, AliAODVertex *origownvtx) const
Float_t * GetPtBinLimits() const
Definition: AliRDHFCuts.h:238
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:239
TH1F * fDCAHist[4 *kMaxPtBins]
! hist. for DCA (Prod Cuts)
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:306
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:270
TH1F * fSigVertHist[4 *kMaxPtBins]
! hist. for sigVert (Prod Cuts)
Bool_t RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod) const
Double_t DecayLength() const
TH3F * fNDLenxyHist3D
! NDlenxy vs Ds mass vs pt
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
TH3F * fCosPiDsHist3D
! cosPiDs vs Ds mass vs pt
Double_t fMassBinSize
range for mass histogram
Double_t fminMass
bin size for inv. mass histo
Int_t GetBackgroundHistoIndex(Int_t iPtBin) const
Float_t fPtLimits[kMaxPtBins+1]