AliPhysics  d565ceb (d565ceb)
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 #include "AliMultSelection.h"
51 
53 ClassImp(AliAnalysisTaskSEDs);
55 
56 //________________________________________________________________________
59  fOutput(0),
60  fHistNEvents(0),
61  fHistoPtWeight(0x0),
62  fPtVsMass(0),
63  fPtVsMassPhi(0),
64  fPtVsMassK0st(0),
65  fYVsPt(0),
66  fYVsPtSig(0),
67  fHistAllV0multNTPCout(0),
68  fHistSelV0multNTPCout(0),
69  fCosPHist3D(0x0),
70  fCosPxyHist3D(0x0),
71  fDLenHist3D(0x0),
72  fDLenxyHist3D(0x0),
73  fNDLenxyHist3D(0x0),
74  fSigVertHist3D(0x0),
75  fDCAHist3D(0x0),
76  fNormIPHist3D(0x0),
77  fCosPiDsHist3D(0x0),
78  fCosPiKPhiHist3D(0x0),
79  fPtProng0Hist3D(0x0),
80  fPtProng1Hist3D(0x0),
81  fPtProng2Hist3D(0x0),
82  fNtupleDs(0),
83  fFillNtuple(0),
84  fUseCentrAxis(0),
85  fSystem(0),
86  fReadMC(kFALSE),
87  fWriteOnlySignal(kFALSE),
88  fDoCutVarHistos(kTRUE),
89  fUseSelectionBit(kFALSE),
90  fFillSparse(kTRUE),
91  fFillSparseDplus(kFALSE),
92  fFillImpParSparse(kFALSE),
93  fFillAcceptanceLevel(kFALSE),
94  fDoRotBkg(kFALSE),
95  fDoBkgPhiSB(kFALSE),
96  fDoCutV0multTPCout(kFALSE),
97  fUseWeight(kFALSE),
98  fUseTrkl(kFALSE),
99  fAODProtection(1),
100  fNPtBins(0),
101  fListCuts(0),
102  fMassRange(0.8),
103  fMassBinSize(0.002),
104  fminMass(1.6),
105  fmaxMass(2.5),
106  fMaxDeltaPhiMass4Rot(0.010),
107  fCounter(0),
108  fAnalysisCuts(0),
109  fnSparse(0),
110  fnSparseIP(0),
111  fImpParSparse(0x0),
112  fMultSelectionObjectName("MultSelection"),
113  fCentEstName("off")
114 {
116 
117  for(Int_t i=0;i<3;i++){
118  fHistCentrality[i]=0;
119  fHistCentralityMult[i]=0;
120  }
121  for(Int_t i=0;i<4;i++) {
122  fChanHist[i]=0;
123  }
124  for(Int_t i=0;i<4*kMaxPtBins;i++){
125  fPtCandHist[i]=0;
126  fMassHist[i]=0;
127  fMassHistPhi[i]=0;
128  fMassHistK0st[i]=0;
129  fCosPHist[i]=0;
130  fDLenHist[i]=0;
131  fSumd02Hist[i]=0;
132  fSigVertHist[i]=0;
133  fPtMaxHist[i]=0;
134  fDCAHist[i]=0;
135  fPtProng0Hist[i]=0;
136  fPtProng1Hist[i]=0;
137  fPtProng2Hist[i]=0;
138  fDalitz[i]=0;
139  fDalitzPhi[i]=0;
140  fDalitzK0st[i]=0;
141  }
142  for(Int_t i=0;i<kMaxPtBins;i++){
143  fMassHistKK[i]=0;
144  fMassHistKpi[i]=0;
145  fMassRotBkgHistPhi[i]=0;
146  fMassLSBkgHistPhi[i]=0;
147  fMassRSBkgHistPhi[i]=0;
148  }
149  for(Int_t i=0;i<kMaxPtBins+1;i++){
150  fPtLimits[i]=0;
151  }
152  for (Int_t i=0; i<4; i++) {
153  fnSparseMC[i]=0;
154  fnSparseMCDplus[i]=0;
155  fImpParSparseMC[i]=0;
156  }
157 }
158 
159 //________________________________________________________________________
160 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name,AliRDHFCutsDstoKKpi* analysiscuts,Int_t fillNtuple):
161  AliAnalysisTaskSE(name),
162  fOutput(0),
163  fHistNEvents(0),
164  fHistoPtWeight(0x0),
165  fPtVsMass(0),
166  fPtVsMassPhi(0),
167  fPtVsMassK0st(0),
168  fYVsPt(0),
169  fYVsPtSig(0),
172  fCosPHist3D(0x0),
173  fCosPxyHist3D(0x0),
174  fDLenHist3D(0x0),
175  fDLenxyHist3D(0x0),
176  fNDLenxyHist3D(0x0),
177  fSigVertHist3D(0x0),
178  fDCAHist3D(0x0),
179  fNormIPHist3D(0x0),
180  fCosPiDsHist3D(0x0),
181  fCosPiKPhiHist3D(0x0),
182  fPtProng0Hist3D(0x0),
183  fPtProng1Hist3D(0x0),
184  fPtProng2Hist3D(0x0),
185  fNtupleDs(0),
186  fFillNtuple(fillNtuple),
187  fUseCentrAxis(0),
188  fSystem(0),
189  fReadMC(kFALSE),
190  fWriteOnlySignal(kFALSE),
191  fDoCutVarHistos(kTRUE),
192  fUseSelectionBit(kFALSE),
193  fFillSparse(kTRUE),
194  fFillSparseDplus(kFALSE),
195  fFillImpParSparse(kFALSE),
196  fFillAcceptanceLevel(kFALSE),
197  fDoRotBkg(kTRUE),
198  fDoBkgPhiSB(kTRUE),
199  fDoCutV0multTPCout(kFALSE),
200  fUseWeight(kFALSE),
201  fUseTrkl(kFALSE),
202  fAODProtection(1),
203  fNPtBins(0),
204  fListCuts(0),
205  fMassRange(0.8),
206  fMassBinSize(0.002),
207  fminMass(1.6),
208  fmaxMass(2.5),
209  fMaxDeltaPhiMass4Rot(0.010),
210  fCounter(0),
211  fAnalysisCuts(analysiscuts),
212  fnSparse(0),
213  fnSparseIP(0),
214  fImpParSparse(0x0),
215  fMultSelectionObjectName("MultSelection"),
216  fCentEstName("off")
217 {
220 
221  for(Int_t i=0;i<3;i++){
222  fHistCentrality[i]=0;
223  fHistCentralityMult[i]=0;
224  }
225  for(Int_t i=0;i<4;i++) {
226  fChanHist[i]=0;
227  }
228  for(Int_t i=0;i<4*kMaxPtBins;i++){
229  fPtCandHist[i]=0;
230  fMassHist[i]=0;
231  fMassHistPhi[i]=0;
232  fMassHistK0st[i]=0;
233  fCosPHist[i]=0;
234  fDLenHist[i]=0;
235  fSumd02Hist[i]=0;
236  fSigVertHist[i]=0;
237  fPtMaxHist[i]=0;
238  fDCAHist[i]=0;
239  fPtProng0Hist[i]=0;
240  fPtProng1Hist[i]=0;
241  fPtProng2Hist[i]=0;
242  fDalitz[i]=0;
243  fDalitzPhi[i]=0;
244  fDalitzK0st[i]=0;
245  }
246  for(Int_t i=0;i<kMaxPtBins;i++){
247  fMassHistKK[i]=0;
248  fMassHistKpi[i]=0;
249  fMassRotBkgHistPhi[i]=0;
250  fMassLSBkgHistPhi[i]=0;
251  fMassRSBkgHistPhi[i]=0;
252  }
253  for(Int_t i=0;i<kMaxPtBins+1;i++){
254  fPtLimits[i]=0;
255  }
256 
257  for (Int_t i=0; i<4; i++) {
258  fnSparseMC[i]=0;
259  fnSparseMCDplus[i]=0;
260  fImpParSparseMC[i]=0;
261  }
262 
265  SetPtBins(nptbins,ptlim);
266 
267  DefineOutput(1,TList::Class()); //My private output
268 
269  DefineOutput(2,TList::Class());
270 
271  DefineOutput(3,AliNormalizationCounter::Class());
272 
273  if(fFillNtuple>0){
274  // Output slot #4 writes into a TNtuple container
275  DefineOutput(4,TNtuple::Class()); //My private output
276  }
277 
278 }
279 
280 //________________________________________________________________________
283  if(n>kMaxPtBins){
284  printf("Max. number of Pt bins = %d\n",kMaxPtBins);
286  fPtLimits[0]=0.;
287  fPtLimits[1]=1.;
288  fPtLimits[2]=3.;
289  fPtLimits[3]=5.;
290  fPtLimits[4]=10.;
291  for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
292  }else{
293  fNPtBins=n;
294  for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
295  for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
296  }
297  if(fDebug > 1){
298  printf("Number of Pt bins = %d\n",fNPtBins);
299  for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);
300  }
301 }
302 //________________________________________________________________________
304 {
305  // Destructor
306  if(fOutput && !fOutput->IsOwner()){
307  delete fHistNEvents;
308  delete fHistAllV0multNTPCout;
309  delete fHistSelV0multNTPCout;
310  delete fCosPHist3D;
311  delete fCosPxyHist3D;
312  delete fDLenHist3D;
313  delete fDLenxyHist3D;
314  delete fNDLenxyHist3D;
315  delete fSigVertHist3D;
316  delete fDCAHist3D;
317  delete fNormIPHist3D;
318  delete fCosPiDsHist3D;
319  delete fCosPiKPhiHist3D;
320  delete fPtProng0Hist3D;
321  delete fPtProng1Hist3D;
322  delete fPtProng2Hist3D;
323 
324  for(Int_t i=0;i<4;i++){
325  delete fChanHist[i];
326  }
327  for(Int_t i=0;i<4*fNPtBins;i++){
328  delete fMassHist[i];
329  delete fMassHistPhi[i];
330  delete fMassHistK0st[i];
331  delete fCosPHist[i];
332  delete fDLenHist[i];
333  delete fSumd02Hist[i];
334  delete fSigVertHist[i];
335  delete fPtMaxHist[i];
336  delete fPtCandHist[i];
337  delete fDCAHist[i];
338  delete fPtProng0Hist[i];
339  delete fPtProng1Hist[i];
340  delete fPtProng2Hist[i];
341  delete fDalitz[i];
342  delete fDalitzPhi[i];
343  delete fDalitzK0st[i];
344  }
345  for(Int_t i=0;i<fNPtBins;i++){
346  delete fMassHistKK[i];
347  delete fMassHistKpi[i];
348  delete fMassRotBkgHistPhi[i];
349  delete fMassLSBkgHistPhi[i];
350  delete fMassRSBkgHistPhi[i];
351  }
352  delete fPtVsMass;
353  delete fPtVsMassPhi;
354  delete fPtVsMassK0st;
355  delete fYVsPt;
356  delete fYVsPtSig;
357  for(Int_t i=0;i<3;i++){
358  delete fHistCentrality[i];
359  delete fHistCentralityMult[i];
360  }
361 
362  delete fnSparse;
363  delete fnSparseIP;
364  if(fFillImpParSparse) {
365  delete fImpParSparse;
366  for(Int_t i=0; i<4; i++) delete fImpParSparseMC[i];
367  }
368  for (Int_t i=0; i<4; i++) {
369  delete fnSparseMC[i];
370  if(fFillSparseDplus)delete fnSparseMCDplus[i];
371  }
372  }
373  if(fHistoPtWeight) delete fHistoPtWeight;
374  delete fOutput;
375  delete fListCuts;
376  delete fNtupleDs;
377  delete fCounter;
378  delete fAnalysisCuts;
379 
380 }
381 
382 //________________________________________________________________________
384  //
385  // set centrality estimator
386  //
387  fUseCentrAxis = flag;
388  if(fUseCentrAxis<kCentOff||fUseCentrAxis>=kCentInvalid) AliWarning("Centrality estimator not valid");
389  if(fUseCentrAxis) {
390  if(fUseCentrAxis==1) fCentEstName = "V0M";
391  if(fUseCentrAxis==2) fCentEstName = "V0A";
392  if(fUseCentrAxis==3) fCentEstName = "CL1";
393  if(fUseCentrAxis==4) fCentEstName = "ZNA";
394  }
395  return;
396 }
397 
398 //________________________________________________________________________
400 {
402 
403  if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
404 
405  fListCuts=new TList();
406  fListCuts->SetOwner();
407  fListCuts->SetName("CutObjects");
408 
410  analysis->SetName("AnalysisCuts");
411 
412  fListCuts->Add(analysis);
413  PostData(2,fListCuts);
414  return;
415 }
416 
417 //________________________________________________________________________
419 {
421  //
422  if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
423 
424  // Several histograms are more conveniently managed in a TList
425  fOutput = new TList();
426  fOutput->SetOwner();
427  fOutput->SetName("OutputHistos");
428 
429  fHistNEvents = new TH1F("hNEvents", "number of events ",15,-0.5,14.5);
430  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsRead");
431  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents Matched dAOD");
432  fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents Mismatched dAOD");
433  fHistNEvents->GetXaxis()->SetBinLabel(4,"nEventsAnal");
434  fHistNEvents->GetXaxis()->SetBinLabel(5,"n. passing IsEvSelected");
435  fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected due to trigger");
436  fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected due to not reco vertex");
437  fHistNEvents->GetXaxis()->SetBinLabel(8,"n. rejected for contr vertex");
438  fHistNEvents->GetXaxis()->SetBinLabel(9,"n. rejected for vertex out of accept");
439  fHistNEvents->GetXaxis()->SetBinLabel(10,"n. rejected for pileup events");
440  fHistNEvents->GetXaxis()->SetBinLabel(11,"no. of out centrality events");
441  fHistNEvents->GetXaxis()->SetBinLabel(12,"no. of 3 prong candidates");
442  fHistNEvents->GetXaxis()->SetBinLabel(13,"no. of Ds after filtering cuts");
443  fHistNEvents->GetXaxis()->SetBinLabel(14,"no. of Ds after selection cuts");
444  fHistNEvents->GetXaxis()->SetBinLabel(15,"no. of not on-the-fly rec Ds");
445 
446  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
447 
448  fHistNEvents->SetMinimum(0);
449  fOutput->Add(fHistNEvents);
450 
451  fHistCentrality[0]=new TH1F("hCentr","centrality",10000,0.,100.);
452  fHistCentrality[1]=new TH1F("hCentr(selectedCent)","centrality(selectedCent)",10000,0.,100.);
453  fHistCentrality[2]=new TH1F("hCentr(OutofCent)","centrality(OutofCent)",10000,0.,100.);
454  fHistCentralityMult[0]=new TH2F("hCentrMult","centrality vs mult",100,0.5,30000.5,40,0.,100.);
455  fHistCentralityMult[1]=new TH2F("hCentrMult(selectedCent)","centrality vs mult(selectedCent)",100,0.5,30000.5,40,0.,100.);
456  fHistCentralityMult[2]=new TH2F("hCentrMult(OutofCent)","centrality vs mult(OutofCent)",100,0.5,30000.5,40,0.,100.);
457  for(Int_t i=0;i<3;i++){
458  fOutput->Add(fHistCentrality[i]);
459  fOutput->Add(fHistCentralityMult[i]);
460  }
461  if(fDoCutV0multTPCout) {
462  fHistAllV0multNTPCout = new TH2F("HistAllV0multNTPCout", "V0mult vs # TPCout (all) ;V0mult ;# TPCout", 1000, 0., 40000, 1000, 0, 30000);
463  fHistSelV0multNTPCout = new TH2F("HistSelV0multNTPCout", "V0mult vs # TPCout (sel) ;V0mult ;# TPCout", 1000, 0., 40000, 1000, 0, 30000);
466  }
467 
468  Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
469 
470  Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);
471  if(nInvMassBins%2==1) nInvMassBins++;
472  // Double_t minMass=massDs-1.0*nInvMassBins*fMassBinSize;
473  Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
474  // Double_t maxMass=massDs+1.0*nInvMassBins*fMassBinSize;
475  Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
476  fminMass = minMass;
477  fmaxMass = maxMass;
478 
479  TString hisname;
480  TString htype;
481  Int_t index;
482  for(Int_t iType=0; iType<4; iType++){
483  for(Int_t i=0;i<fNPtBins;i++){
484  if(iType==0){
485  htype="All";
486  index=GetHistoIndex(i);
487  }else if(iType==1){
488  htype="Sig";
489  index=GetSignalHistoIndex(i);
490  }else if(iType==2){
491  htype="Bkg";
492  index=GetBackgroundHistoIndex(i);
493  }else{
494  htype="ReflSig";
495  index=GetReflSignalHistoIndex(i);
496  }
497  hisname.Form("hMass%sPt%d",htype.Data(),i);
498  fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
499  fMassHist[index]->Sumw2();
500  hisname.Form("hMass%sPt%dphi",htype.Data(),i);
501  fMassHistPhi[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
502  fMassHistPhi[index]->Sumw2();
503  hisname.Form("hMass%sPt%dk0st",htype.Data(),i);
504  fMassHistK0st[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
505  fMassHistK0st[index]->Sumw2();
506  hisname.Form("hCosP%sPt%d",htype.Data(),i);
507  fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
508  hisname.Form("hCosPxy%sPt%d",htype.Data(),i);
509  fCosPxyHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
510  hisname.Form("hDLen%sPt%d",htype.Data(),i);
511  fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
512  hisname.Form("hDLenxy%sPt%d",htype.Data(),i);
513  fDLenxyHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
514  hisname.Form("hNDLenxy%sPt%d",htype.Data(),i);
515  fNDLenxyHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,11.);
516  hisname.Form("hSumd02%sPt%d",htype.Data(),i);
517  fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
518  hisname.Form("hSigVert%sPt%d",htype.Data(),i);
519  fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
520  hisname.Form("hPtMax%sPt%d",htype.Data(),i);
521  fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
522  hisname.Form("hPtCand%sPt%d",htype.Data(),i);
523  fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
524  hisname.Form("hDCA%sPt%d",htype.Data(),i);
525  fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
526  hisname.Form("hNormIP%sPt%d",htype.Data(),i);
527  fNormIPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,6.);
528  hisname.Form("hCosPiDs%sPt%d",htype.Data(),i);
529  fCosPiDsHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
530  hisname.Form("hCosPiKPhi%sPt%d",htype.Data(),i);
531  fCosPiKPhiHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
532  hisname.Form("hPtProng0%sPt%d",htype.Data(),i);
533  fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
534  hisname.Form("hPtProng1%sPt%d",htype.Data(),i);
535  fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
536  hisname.Form("hPtProng2%sPt%d",htype.Data(),i);
537  fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
538  hisname.Form("hDalitz%sPt%d",htype.Data(),i);
539  fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
540  hisname.Form("hDalitz%sPt%dphi",htype.Data(),i);
541  fDalitzPhi[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
542  hisname.Form("hDalitz%sPt%dk0st",htype.Data(),i);
543  fDalitzK0st[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
544  }
545  }
546 
547  for(Int_t i=0; i<4*fNPtBins; i++){
548  fOutput->Add(fMassHist[i]);
549  fOutput->Add(fMassHistPhi[i]);
550  fOutput->Add(fMassHistK0st[i]);
551  fOutput->Add(fPtCandHist[i]);
552  if(fDoCutVarHistos){
553  fOutput->Add(fCosPHist[i]);
554  fOutput->Add(fCosPxyHist[i]);
555  fOutput->Add(fDLenHist[i]);
556  fOutput->Add(fDLenxyHist[i]);
557  fOutput->Add(fNDLenxyHist[i]);
558  fOutput->Add(fSumd02Hist[i]);
559  fOutput->Add(fSigVertHist[i]);
560  fOutput->Add(fPtMaxHist[i]);
561  fOutput->Add(fDCAHist[i]);
562  fOutput->Add(fNormIPHist[i]);
563  fOutput->Add(fCosPiDsHist[i]);
564  fOutput->Add(fCosPiKPhiHist[i]);
565  fOutput->Add(fPtProng0Hist[i]);
566  fOutput->Add(fPtProng1Hist[i]);
567  fOutput->Add(fPtProng2Hist[i]);
568  fOutput->Add(fDalitz[i]);
569  fOutput->Add(fDalitzPhi[i]);
570  fOutput->Add(fDalitzK0st[i]);
571  }
572  }
573 
574  fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",64,-0.5,63.5);
575  fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",64,-0.5,63.5);
576  fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",64,-0.5,63.5);
577  fChanHist[3] = new TH1F("hChanReflSig", "KKpi and piKK candidates",64,-0.5,63.5);
578  for(Int_t i=0;i<4;i++){
579  fChanHist[i]->SetMinimum(0);
580  fOutput->Add(fChanHist[i]);
581  }
582 
583  fCosPHist3D = new TH3F("fCosPHist3D","CosP vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.5,1.);
584  fCosPxyHist3D = new TH3F("fCosPxyHist3D","CosPxy vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.5,1.);
585  fDLenHist3D = new TH3F("fDLenHist3D","DLen vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,0.5);
586  fDLenxyHist3D = new TH3F("fDLenxyHist3D","DLenxy vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,0.5);
587  fNDLenxyHist3D = new TH3F("fNDLenxyHist3D","NDLenxy vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,11.);
588  fSigVertHist3D = new TH3F("fSigVertHist3D","SigVert vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,0.1);
589  fDCAHist3D = new TH3F("fDCAHist3D","DCA vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,0.1);
590  fNormIPHist3D = new TH3F("fNormIPHist3D","nIP vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,6.);
591  fCosPiDsHist3D = new TH3F("fCosPiDsHist3D","CosPiDs vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.5,1.);
592  fCosPiKPhiHist3D = new TH3F("fCosPiKPhiHist3D","CosPiKPhi vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.,0.5);
593  fPtProng0Hist3D = new TH3F("fPtProng0Hist3D","Pt prong0 vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.0,20.);
594  fPtProng1Hist3D = new TH3F("fPtProng1Hist3D","Pt prong1 vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.0,20.);
595  fPtProng2Hist3D = new TH3F("fPtProng2Hist3D","Pt prong2 vs Ds mass",nInvMassBins,minMass,maxMass,20,0.,20.,100,0.0,20.);
596 
597  if(!fReadMC && fDoCutVarHistos) {
598  fOutput->Add(fCosPHist3D);
599  fOutput->Add(fCosPxyHist3D);
600  fOutput->Add(fDLenHist3D);
601  fOutput->Add(fDLenxyHist3D);
602  fOutput->Add(fNDLenxyHist3D);
603  fOutput->Add(fSigVertHist3D);
604  fOutput->Add(fDCAHist3D);
605  fOutput->Add(fNormIPHist3D);
606  fOutput->Add(fCosPiDsHist3D);
608  fOutput->Add(fPtProng0Hist3D);
609  fOutput->Add(fPtProng1Hist3D);
610  fOutput->Add(fPtProng2Hist3D);
611  }
612 
613  fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nInvMassBins,minMass,maxMass,40,0.,20.);
614  fPtVsMassPhi=new TH2F("hPtVsMassPhi","PtVsMass (phi selection)",nInvMassBins,minMass,maxMass,200,0.,20.);
615  fPtVsMassK0st=new TH2F("hPtVsMassK0st","PtVsMass (K0* selection)",nInvMassBins,minMass,maxMass,200,0.,20.);
616  fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
617  fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
618 
619  for(Int_t i=0;i<fNPtBins;i++){
620  hisname.Form("hMassKKPt%d",i);
621  fMassHistKK[i]=new TH1F(hisname.Data(),hisname.Data(),200,0.95,1.35);
622  fMassHistKK[i]->Sumw2();
623  fOutput->Add(fMassHistKK[i]);
624  hisname.Form("hMassKpiPt%d",i);
625  fMassHistKpi[i]=new TH1F(hisname.Data(),hisname.Data(),200,0.7,1.1);
626  fMassHistKpi[i]->Sumw2();
627  fOutput->Add(fMassHistKpi[i]);
628  if(fDoRotBkg) {
629  hisname.Form("hMassAllPt%dphi_RotBkg",i);
630  fMassRotBkgHistPhi[i]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
631  fMassRotBkgHistPhi[i]->Sumw2();
632  fOutput->Add(fMassRotBkgHistPhi[i]);
633  }
634  if(fDoBkgPhiSB) {
635  hisname.Form("fMassLSBkgHistPhiPt%d",i);
636  fMassLSBkgHistPhi[i]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
637  fMassLSBkgHistPhi[i]->Sumw2();
638  fOutput->Add(fMassLSBkgHistPhi[i]);
639  hisname.Form("fMassRSBkgHistPhiPt%d",i);
640  fMassRSBkgHistPhi[i]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
641  fMassRSBkgHistPhi[i]->Sumw2();
642  fOutput->Add(fMassRSBkgHistPhi[i]);
643  }
644  }
645 
646  fOutput->Add(fPtVsMass);
647  fOutput->Add(fPtVsMassPhi);
648  fOutput->Add(fPtVsMassK0st);
649  fOutput->Add(fYVsPt);
650  fOutput->Add(fYVsPtSig);
651 
652  //Sparses for Cut variation and IP studies
653  if(fFillSparse) {
655  if(fReadMC) CreateIPSparse();
656  }
657 
658  //Sparses for Impact parameter fits
660 
661  //Counter for Normalization
662  fCounter = new AliNormalizationCounter("NormalizationCounter");
663  fCounter->Init();
664 
665  PostData(1,fOutput);
666  PostData(3,fCounter);
667 
668  if(fFillNtuple>0 && fFillNtuple<4){
669  OpenFile(4); // 4 is the slot number of the ntuple
670 
671  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");
672  }
673  else if(fFillNtuple==4) {
674  OpenFile(4); // 4 is the slot number of the ntuple
675 
676  fNtupleDs = new TNtuple("fNtupleDs","Ds","Pt:InvMass:d0:origin");
677  }
678 
679  return;
680 }
681 
682 //________________________________________________________________________
684 {
687 
688  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
689 
690  fHistNEvents->Fill(0); // all events
691  if(fAODProtection>=0){
692  // Protection against different number of events in the AOD and deltaAOD
693  // In case of discrepancy the event is rejected.
694  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
695  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
696  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
697  fHistNEvents->Fill(2);
698  return;
699  }
700  fHistNEvents->Fill(1);
701  }
702 
703  TClonesArray *array3Prong = 0;
704  if(!aod && AODEvent() && IsStandardAOD()) {
705  // In case there is an AOD handler writing a standard AOD, use the AOD
706  // event in memory rather than the input (ESD) event.
707  aod = dynamic_cast<AliAODEvent*> (AODEvent());
708  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
709  // have to taken from the AOD event hold by the AliAODExtension
710  AliAODHandler* aodHandler = (AliAODHandler*)
711  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
712  if(aodHandler->GetExtensions()) {
713  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
714  AliAODEvent *aodFromExt = ext->GetAOD();
715  array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
716  }
717  } else if(aod) {
718  array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
719  }
720 
721  if(!aod || !array3Prong) {
722  printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
723  return;
724  }
725 
726 
727  // fix for temporary bug in ESDfilter
728  // the AODs with null vertex pointer didn't pass the PhysSel
729  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
730 
731 
732 
733  fHistNEvents->Fill(3); // count event
734  // Post the data already here
735  PostData(1,fOutput);
736 
738 
739 
740  Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
741  Float_t ntracks=aod->GetNumberOfTracks();
742  Float_t evCentr=fAnalysisCuts->GetCentrality(aod);
743 
744  fHistCentrality[0]->Fill(evCentr);
745  fHistCentralityMult[0]->Fill(ntracks,evCentr);
752  fHistNEvents->Fill(10);
753  fHistCentrality[2]->Fill(evCentr);
754  fHistCentralityMult[2]->Fill(ntracks,evCentr);
755  }
756 
758  Int_t runNumber=aod->GetRunNumber();
759 
760 
761 
762  TClonesArray *arrayMC=0;
763  AliAODMCHeader *mcHeader=0;
764 
765  // AOD primary vertex
766  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
767  // vtx1->Print();
768 
769  // load MC particles
770  if(fReadMC){
771 
772  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
773  if(!arrayMC) {
774  printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
775  return;
776  }
777 
778  // load MC header
779  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
780  if(!mcHeader) {
781  printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
782  return;
783  }
784  }
785 
787 
788  if(fReadMC && fFillSparse){
789  if(aod->GetTriggerMask()==0 && (runNumber>=195344 && runNumber<=195677))
790  // protection for events with empty trigger mask in p-Pb
791  return;
793  // events not passing the centrality selection can be removed immediately.
794  return;
795  Double_t zMCVertex = mcHeader->GetVtxZ();
796  if (TMath::Abs(zMCVertex) > fAnalysisCuts->GetMaxVtxZ())
797  return;
798  FillMCGenAccHistos(arrayMC, mcHeader, nTracklets);
799  }
800 
801  if(!isEvSel)return;
802 
803  fHistNEvents->Fill(4);
804  fHistCentrality[1]->Fill(evCentr);
805  fHistCentralityMult[1]->Fill(ntracks,evCentr);
806 
807  if(fDoCutV0multTPCout) {
808  //cut on V0mult. vs #tracks kTPCout
809  Float_t V0mult=0.;
810  AliAODVZERO *aodVZERO = (AliAODVZERO*)aod->GetVZEROData();
811  if(aodVZERO) {
812  for(int ich=0; ich<64;ich++) V0mult += aodVZERO->GetMultiplicity(ich);
813  }
814  Int_t nTPCout=0;
815  for (Int_t i=0;i<aod->GetNumberOfTracks();++i) {
816  AliVTrack *track = aod->GetTrack(i);
817  if (!track) continue;
818  if((track->GetStatus() & AliVTrack::kTPCout)) nTPCout++;
819  }
820  fHistAllV0multNTPCout->Fill(V0mult,nTPCout);
821  if(nTPCout > (0.32*V0mult+750)) return;
822  else fHistSelV0multNTPCout->Fill(V0mult,nTPCout);
823  }
824 
825  Int_t n3Prong = array3Prong->GetEntriesFast();
826  if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
827 
828  Int_t pdgDstoKKpi[3]={321,321,211};
829  Int_t nSelected=0;
830  Int_t nFiltered=0;
831  Double_t massPhi=TDatabasePDG::Instance()->GetParticle(333)->Mass();
832 
833  // vHF object is needed to call the method that refills the missing info of the candidates
834  // if they have been deleted in dAOD reconstruction phase
835  // in order to reduce the size of the file
837 
838  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
839 
840  AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
841  fHistNEvents->Fill(11);
842 
844  continue;
845  }
846  nFiltered++;
847  fHistNEvents->Fill(12);
848 
849  if(!(vHF->FillRecoCand(aod,d))) {
850  fHistNEvents->Fill(14); //monitor how often this fails
851  continue;
852  }
853 
854  Bool_t unsetvtx=kFALSE;
855  if(!d->GetOwnPrimaryVtx()){
856  d->SetOwnPrimaryVtx(vtx1);
857  unsetvtx=kTRUE;
858  // NOTE: the ow primary vertex should be unset, otherwise there is a memory leak
859  // Pay attention if you use continue inside this loop!!!
860  }
861 
862  Bool_t recVtx=kFALSE;
863  AliAODVertex *origownvtx=0x0;
864 
865  Double_t ptCand = d->Pt();
866  Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
867  Double_t rapid=d->YDs();
868  fYVsPt->Fill(ptCand,rapid);
869  Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);
870 
871  if(isFidAcc){
872 
873  Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
874  Int_t retCodeNoRes=retCodeAnalysisCuts;
876  if(origRes){
878  retCodeNoRes=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
880  }
881  Double_t massKK = 0.;
882 
883  if(retCodeNoRes&1){ //KKpi
884  massKK=d->InvMass2Prongs(0,1,321,321);
885  Double_t massKp=d->InvMass2Prongs(1,2,321,211);
886  fMassHistKK[iPtBin]->Fill(massKK);
887  fMassHistKpi[iPtBin]->Fill(massKp);
888  }
889  if(retCodeNoRes&2){ //piKK
890  massKK=d->InvMass2Prongs(1,2,321,321);
891  Double_t massKp=d->InvMass2Prongs(0,1,211,321);
892  fMassHistKK[iPtBin]->Fill(massKK);
893  fMassHistKpi[iPtBin]->Fill(massKp);
894  }
895 
896  Int_t isKKpi=retCodeAnalysisCuts&1;
897  Int_t ispiKK=retCodeAnalysisCuts&2;
898  Int_t isPhiKKpi=retCodeAnalysisCuts&4;
899  Int_t isPhipiKK=retCodeAnalysisCuts&8;
900  Int_t isK0starKKpi=retCodeAnalysisCuts&16;
901  Int_t isK0starpiKK=retCodeAnalysisCuts&32;
902 
903  if(retCodeAnalysisCuts>0){
905  if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
906  if(fAnalysisCuts->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
907  else fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
908  }
909 
910 
911  fHistNEvents->Fill(13);
912  nSelected++;
913 
914  Int_t index=GetHistoIndex(iPtBin);
915  fPtCandHist[index]->Fill(ptCand);
916 
917 
918  Double_t weightKKpi=1.;
919  Double_t weightpiKK=1.;
921  weightKKpi=fAnalysisCuts->GetWeightForKKpi();
922  weightpiKK=fAnalysisCuts->GetWeightForpiKK();
923  if(weightKKpi>1. || weightKKpi<0.) weightKKpi=0.;
924  if(weightpiKK>1. || weightpiKK<0.) weightpiKK=0.;
925  }
926 
927  fChanHist[0]->Fill(retCodeAnalysisCuts);
928 
929 
930  Double_t invMass = 0.;
931  Int_t indexMCKKpi=-1;
932  Int_t indexMCpiKK=-1;
933  Int_t labDs=-1;
934  Int_t labDplus=-1;
935  Int_t pdgCode0=-999;
936  Int_t isMCSignal=-1;
937 
938 
939  if(fReadMC){
940 
941  labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
942  labDplus = d->MatchToMC(411,arrayMC,3,pdgDstoKKpi);
943  if(labDs>=0){
944  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
945  AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
946  pdgCode0=TMath::Abs(p->GetPdgCode());
947 
948  if(isKKpi){
949  if(pdgCode0==321) {
950  indexMCKKpi=GetSignalHistoIndex(iPtBin);
951  fYVsPtSig->Fill(ptCand,rapid);
952  fChanHist[1]->Fill(retCodeAnalysisCuts);
953  isMCSignal=1;
954  }else{
955  indexMCKKpi=GetReflSignalHistoIndex(iPtBin);
956  fChanHist[3]->Fill(retCodeAnalysisCuts);
957  isMCSignal=0;
958  }
959  }
960  if(ispiKK){
961  if(pdgCode0==211) {
962  indexMCpiKK=GetSignalHistoIndex(iPtBin);
963  fYVsPtSig->Fill(ptCand,rapid);
964  fChanHist[1]->Fill(retCodeAnalysisCuts);
965  isMCSignal=1;
966  }else{
967  indexMCpiKK=GetReflSignalHistoIndex(iPtBin);
968  fChanHist[3]->Fill(retCodeAnalysisCuts);
969  isMCSignal=0;
970  }
971  }
972  }else{
973  indexMCpiKK=GetBackgroundHistoIndex(iPtBin);
974  indexMCKKpi=GetBackgroundHistoIndex(iPtBin);
975  fChanHist[2]->Fill(retCodeAnalysisCuts);
976  if(labDplus>=0) {
977  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
978  AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
979  pdgCode0=TMath::Abs(p->GetPdgCode());
980  }
981  }
982  }
983 
984  Double_t candType = 0.5; //for bkg
985  Float_t trueImpParDsFromB = 99999.;
986  if(isKKpi){
987  if(fDoRotBkg && TMath::Abs(massKK-massPhi)<=fMaxDeltaPhiMass4Rot)GenerateRotBkg(d,1,iPtBin);
988 
989  invMass=d->InvMassDsKKpi();
990  fMassHist[index]->Fill(invMass,weightKKpi);
991  fPtVsMass->Fill(invMass,ptCand,weightKKpi);
992 
993  if(fDoBkgPhiSB && (0.010<TMath::Abs(massKK-massPhi)) && (TMath::Abs(massKK-massPhi)<0.030) ) {
994  if(massKK<massPhi)fMassLSBkgHistPhi[iPtBin]->Fill(invMass);
995  else fMassRSBkgHistPhi[iPtBin]->Fill(invMass);
996  }
997 
998  if(isPhiKKpi){
999  fMassHistPhi[index]->Fill(invMass,weightKKpi);
1000  fPtVsMassPhi->Fill(invMass,ptCand,weightKKpi);
1001  }
1002  if(isK0starKKpi){
1003  fMassHistK0st[index]->Fill(invMass,weightKKpi);
1004  fPtVsMassK0st->Fill(invMass,ptCand,weightKKpi);
1005  }
1006  if(fReadMC && indexMCKKpi!=-1){
1007  fMassHist[indexMCKKpi]->Fill(invMass,weightKKpi);
1008  if(isPhiKKpi) {
1009  fMassHistPhi[indexMCKKpi]->Fill(invMass,weightKKpi);
1010  if(fFillSparse) {
1011  if(indexMCKKpi==GetSignalHistoIndex(iPtBin) || labDplus >= 0) {
1012  AliAODMCParticle *partDs;
1013  if(indexMCKKpi==GetSignalHistoIndex(iPtBin)) {
1014  partDs = (AliAODMCParticle*)arrayMC->At(labDs);
1015  }
1016  if(labDplus >= 0) partDs = (AliAODMCParticle*)arrayMC->At(labDplus);
1017  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partDs,kTRUE);
1018  if(orig==4) {
1019  candType = 1.5;
1020  }
1021  if(orig==5) {
1022  candType = 2.5;
1023  if(isPhiKKpi && fFillImpParSparse) {trueImpParDsFromB = GetTrueImpactParameterDstoPhiPi(mcHeader,arrayMC,partDs)*10000;}
1024  }
1025  }
1026  }
1027  }
1028  if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass,weightKKpi);
1029  }
1030  if(isPhiKKpi && fFillImpParSparse) {
1031  Double_t impParxy = d->ImpParXY()*10000.;
1032  Double_t array4ImpPar[3] = {invMass,ptCand,impParxy};
1033  if(!fReadMC) fImpParSparse->Fill(array4ImpPar);
1034  else {
1035  if(candType == 1.5 && indexMCKKpi==GetSignalHistoIndex(iPtBin)) fImpParSparseMC[0]->Fill(array4ImpPar);
1036  else if(candType == 2.5 && indexMCKKpi==GetSignalHistoIndex(iPtBin)) {
1037  fImpParSparseMC[1]->Fill(array4ImpPar);
1038  Double_t array4ImpParTrueB[3] = {invMass,ptCand,trueImpParDsFromB};
1039  fImpParSparseMC[2]->Fill(array4ImpParTrueB);
1040  }
1041  else fImpParSparseMC[3]->Fill(array4ImpPar);
1042  }
1043  }
1044  if(isPhiKKpi && fFillNtuple==4) {
1045  Float_t impParxy = d->ImpParXY()*10000.;
1046  Float_t tmp[4] = {(Float_t)ptCand,(Float_t)invMass,impParxy,(Float_t)candType};
1047  fNtupleDs->Fill(tmp);
1048  PostData(4,fNtupleDs);
1049  }
1050  }
1051  if(ispiKK){
1052  if(fDoRotBkg && TMath::Abs(massKK-massPhi)<=fMaxDeltaPhiMass4Rot)GenerateRotBkg(d,2,iPtBin);
1053 
1054  invMass=d->InvMassDspiKK();
1055  fMassHist[index]->Fill(invMass,weightpiKK);
1056  fPtVsMass->Fill(invMass,ptCand,weightpiKK);
1057 
1058  if(fDoBkgPhiSB && (0.010<TMath::Abs(massKK-massPhi)) && (TMath::Abs(massKK-massPhi)<0.030) ) {
1059  if(massKK<massPhi)fMassLSBkgHistPhi[iPtBin]->Fill(invMass);
1060  else fMassRSBkgHistPhi[iPtBin]->Fill(invMass);
1061  }
1062 
1063  if(isPhipiKK){
1064  fMassHistPhi[index]->Fill(invMass,weightpiKK);
1065  fPtVsMassPhi->Fill(invMass,ptCand,weightpiKK);
1066  }
1067  if(isK0starpiKK){
1068  fMassHistK0st[index]->Fill(invMass,weightpiKK);
1069  fPtVsMassK0st->Fill(invMass,ptCand,weightpiKK);
1070  }
1071  if(fReadMC && indexMCpiKK!=-1){
1072  fMassHist[indexMCpiKK]->Fill(invMass,weightpiKK);
1073  if(isPhipiKK) {
1074  fMassHistPhi[indexMCpiKK]->Fill(invMass,weightpiKK);
1075  if(fFillSparse) {
1076  if(indexMCpiKK==GetSignalHistoIndex(iPtBin) || labDplus >= 0) {
1077  AliAODMCParticle *partDs;
1078  if(indexMCpiKK==GetSignalHistoIndex(iPtBin)) partDs = (AliAODMCParticle*)arrayMC->At(labDs);
1079  if(labDplus >= 0) partDs = (AliAODMCParticle*)arrayMC->At(labDplus);
1080  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partDs,kTRUE);
1081  if(orig==4) {
1082  candType = 1.5;
1083  }
1084  if(orig==5) {
1085  candType = 2.5;
1086  if(isPhipiKK && fFillImpParSparse) {trueImpParDsFromB = GetTrueImpactParameterDstoPhiPi(mcHeader,arrayMC,partDs)*10000;}
1087  }
1088  }
1089  }
1090  }
1091  if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass,weightpiKK);
1092  }
1093  if(isPhipiKK && fFillImpParSparse) {
1094  Double_t impParxy = d->ImpParXY()*10000.;
1095  Double_t array4ImpPar[3] = {invMass,ptCand,impParxy};
1096  if(!fReadMC) fImpParSparse->Fill(array4ImpPar);
1097  else {
1098  if(candType == 1.5 && indexMCpiKK==GetSignalHistoIndex(iPtBin)) fImpParSparseMC[0]->Fill(array4ImpPar);
1099  else if(candType == 2.5 && indexMCpiKK==GetSignalHistoIndex(iPtBin)) {
1100  fImpParSparseMC[1]->Fill(array4ImpPar);
1101  Double_t array4ImpParTrueB[3] = {invMass,ptCand,trueImpParDsFromB};
1102  fImpParSparseMC[2]->Fill(array4ImpParTrueB);
1103  }
1104  else fImpParSparseMC[3]->Fill(array4ImpPar);
1105  }
1106  }
1107  if(isPhipiKK && fFillNtuple==4) {
1108  Float_t impParxy = d->ImpParXY()*10000.;
1109  Float_t tmp[4] = {(Float_t)ptCand,(Float_t)invMass,impParxy,(Float_t)candType};
1110  fNtupleDs->Fill(tmp);
1111  PostData(4,fNtupleDs);
1112  }
1113  }
1114 
1116 
1117  const Int_t nProng = 3;
1118  Double_t deltaMassKK=999.;
1119  Double_t dlen=d->DecayLength();
1120  Double_t dlenxy=d->DecayLengthXY();
1121  // Double_t normdl=d->NormalizedDecayLength();
1122  Double_t normdlxy=d->NormalizedDecayLengthXY();
1123  Double_t cosp=d->CosPointingAngle();
1124  Double_t cospxy=d->CosPointingAngleXY();
1125  Double_t pt0=d->PtProng(0);
1126  Double_t pt1=d->PtProng(1);
1127  Double_t pt2=d->PtProng(2);
1128  Double_t sigvert=d->GetSigmaVert();
1129  Double_t cosPiDs=-99.;
1130  Double_t cosPiKPhi=-99.;
1131  Double_t normIP=-999.; //to store the maximum topomatic var. among the 3 prongs
1132  Double_t normIPprong[nProng]; //to store IP of k,k,pi
1133  Double_t absimpparxy=TMath::Abs(d->ImpParXY());
1134  for(Int_t ijp=0; ijp<nProng; ijp++) normIPprong[ijp]=-999.;
1135 
1136  Double_t ptWeight = 1.;
1137  if(fFillSparse) {
1138 
1139  AliMultSelection *multSelection = (AliMultSelection*)aod->FindListObject(fMultSelectionObjectName);
1140  Double_t cent = -999;
1141  if(multSelection && fUseCentrAxis) {
1142  cent = multSelection->GetMultiplicityPercentile(Form("%s",fCentEstName.Data()));
1143  Int_t qual = multSelection->GetEvSelCode();
1144  if(qual == 199 ) cent=-999;
1145  }
1146  else cent = 0.5;
1147 
1148  if (fUseWeight && fHistoPtWeight){
1149  AliDebug(2,"Using Histogram as Pt weight function");
1150  ptWeight = GetPtWeightFromHistogram(ptCand);
1151  }
1152  if(isPhiKKpi) {
1153  Double_t tmpNormIP[nProng];
1154 
1155  invMass = d->InvMassDsKKpi();
1156  massKK = d->InvMass2Prongs(0,1,321,321);
1157  deltaMassKK = TMath::Abs(massKK-massPhi);
1158  cosPiDs = d->CosPiDsLabFrameKKpi();
1159  cosPiKPhi = d->CosPiKPhiRFrameKKpi();
1160  cosPiKPhi = TMath::Abs(cosPiKPhi*cosPiKPhi*cosPiKPhi);
1161 
1162  for(Int_t ip=0; ip<nProng; ip++) {
1163  Double_t diffIP, errdiffIP;
1164  d->Getd0MeasMinusExpProng(ip,aod->GetMagneticField(),diffIP,errdiffIP);
1165  tmpNormIP[ip] = diffIP/errdiffIP;
1166  if(ip==0) normIP = tmpNormIP[ip];
1167  else if(TMath::Abs(tmpNormIP[ip])>TMath::Abs(normIP)) normIP = tmpNormIP[ip];
1168  }
1169  normIPprong[0] = tmpNormIP[0];
1170  normIPprong[1] = tmpNormIP[1];
1171  normIPprong[2] = tmpNormIP[2];
1172 
1173  Double_t var4nSparse[knVarForSparse] = {invMass,ptCand,deltaMassKK*1000,dlen*1000,dlenxy*1000,normdlxy,cosp*100,cospxy*100,
1174  sigvert*1000,cosPiDs*10,cosPiKPhi*10,TMath::Abs(normIP),nTracklets,cent,absimpparxy*10000};
1175 
1176  if(!fReadMC) {
1177  fnSparse->Fill(var4nSparse);
1178  }
1179  else {
1180  if(indexMCKKpi==GetSignalHistoIndex(iPtBin)) {
1181  if(candType==1.5) fnSparseMC[2]->Fill(var4nSparse,ptWeight);
1182  if(candType==2.5) fnSparseMC[3]->Fill(var4nSparse,ptWeight);
1183  }
1184  else if(fFillSparseDplus && labDplus>=0 && pdgCode0==321) {
1185  if(candType==1.5) fnSparseMCDplus[2]->Fill(var4nSparse,ptWeight);
1186  if(candType==2.5) fnSparseMCDplus[3]->Fill(var4nSparse,ptWeight);
1187  }
1188  }
1189  }
1190  if(isPhipiKK) {
1191  Double_t tmpNormIP[nProng];
1192 
1193  invMass = d->InvMassDspiKK();
1194  massKK = d->InvMass2Prongs(1,2,321,321);
1195  deltaMassKK = TMath::Abs(massKK-massPhi);
1196  cosPiDs = d->CosPiDsLabFramepiKK();
1197  cosPiKPhi = d->CosPiKPhiRFramepiKK();
1198  cosPiKPhi = TMath::Abs(cosPiKPhi*cosPiKPhi*cosPiKPhi);
1199 
1200  for(Int_t ip=0; ip<nProng; ip++) {
1201  Double_t diffIP, errdiffIP;
1202  d->Getd0MeasMinusExpProng(ip,aod->GetMagneticField(),diffIP,errdiffIP);
1203  tmpNormIP[ip] = diffIP/errdiffIP;
1204  if(ip==0) normIP = tmpNormIP[ip];
1205  else if(TMath::Abs(tmpNormIP[ip])>TMath::Abs(normIP)) normIP = tmpNormIP[ip];
1206  }
1207 
1208  normIPprong[0] = tmpNormIP[2];
1209  normIPprong[1] = tmpNormIP[1];
1210  normIPprong[2] = tmpNormIP[0];
1211 
1212  Double_t var4nSparse[knVarForSparse] = {invMass,ptCand,deltaMassKK*1000,dlen*1000,dlenxy*1000,normdlxy,cosp*100,cospxy*100,
1213  sigvert*1000,cosPiDs*10,cosPiKPhi*10,TMath::Abs(normIP),nTracklets,cent,absimpparxy*10000};
1214 
1215  if(!fReadMC) {
1216  fnSparse->Fill(var4nSparse);
1217  }
1218  else {
1219  if(indexMCpiKK==GetSignalHistoIndex(iPtBin)) {
1220  if(candType==1.5) fnSparseMC[2]->Fill(var4nSparse,ptWeight);
1221  if(candType==2.5) fnSparseMC[3]->Fill(var4nSparse,ptWeight);
1222  }
1223  else if(fFillSparseDplus && labDplus>=0 && pdgCode0==211) {
1224  if(candType==1.5) fnSparseMCDplus[2]->Fill(var4nSparse,ptWeight);
1225  if(candType==2.5) fnSparseMCDplus[3]->Fill(var4nSparse,ptWeight);
1226  }
1227  }
1228  }
1229 
1230  if(fReadMC && (isPhiKKpi || isPhiKKpi)) {
1231  Double_t var[6] = {ptCand,normIP,normIPprong[0],normIPprong[1],normIPprong[2],candType};
1232  fnSparseIP->Fill(var);
1233  }
1234  }
1236 
1237  if(fDoCutVarHistos){
1238  Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
1239  Double_t dca=d->GetDCA();
1240  Double_t ptmax=0;
1241 
1242  for(Int_t i=0;i<3;i++){
1243  if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
1244  }
1245  fCosPHist[index]->Fill(cosp);
1246  fCosPxyHist[index]->Fill(cospxy);
1247  fDLenHist[index]->Fill(dlen);
1248  fDLenxyHist[index]->Fill(dlenxy);
1249  fNDLenxyHist[index]->Fill(normdlxy);
1250  fSigVertHist[index]->Fill(sigvert);
1251  fSumd02Hist[index]->Fill(sumD02);
1252  fPtMaxHist[index]->Fill(ptmax);
1253  fDCAHist[index]->Fill(dca);
1254  fNormIPHist[index]->Fill(normIP);
1255  fCosPiDsHist[index]->Fill(cosPiDs);
1256  fCosPiKPhiHist[index]->Fill(cosPiKPhi);
1257  fPtProng0Hist[index]->Fill(pt0);
1258  fPtProng1Hist[index]->Fill(pt1);
1259  fPtProng2Hist[index]->Fill(pt2);
1260  if(!fReadMC) {
1261  fCosPHist3D->Fill(invMass,ptCand,cosp);
1262  fCosPxyHist3D->Fill(invMass,ptCand,cospxy);
1263  fDLenHist3D->Fill(invMass,ptCand,dlen);
1264  fDLenxyHist3D->Fill(invMass,ptCand,dlenxy);
1265  fNDLenxyHist3D->Fill(invMass,ptCand,normdlxy);
1266  fSigVertHist3D->Fill(invMass,ptCand,sigvert);
1267  fDCAHist3D->Fill(invMass,ptCand,dca);
1268  fNormIPHist3D->Fill(invMass,ptCand,normIP);
1269  fCosPiDsHist3D->Fill(invMass,ptCand,cosPiDs);
1270  fCosPiKPhiHist3D->Fill(invMass,ptCand,cosPiKPhi);
1271  fPtProng0Hist3D->Fill(invMass,ptCand,pt0);
1272  fPtProng1Hist3D->Fill(invMass,ptCand,pt1);
1273  fPtProng2Hist3D->Fill(invMass,ptCand,pt2);
1274  }
1275  if(isKKpi){
1276  Double_t massKK=d->InvMass2Prongs(0,1,321,321);
1277  Double_t massKp=d->InvMass2Prongs(1,2,321,211);
1278  fDalitz[index]->Fill(massKK,massKp);
1279  if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);
1280  if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);
1281  if(fReadMC && indexMCKKpi!=-1){
1282  fDalitz[indexMCKKpi]->Fill(massKK,massKp);
1283  if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);
1284  if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);
1285  fCosPHist[indexMCKKpi]->Fill(cosp);
1286  fCosPxyHist[indexMCKKpi]->Fill(cospxy);
1287  fDLenHist[indexMCKKpi]->Fill(dlen);
1288  fDLenxyHist[indexMCKKpi]->Fill(dlenxy);
1289  fNDLenxyHist[indexMCKKpi]->Fill(normdlxy);
1290  fSigVertHist[indexMCKKpi]->Fill(sigvert);
1291  fSumd02Hist[indexMCKKpi]->Fill(sumD02);
1292  fPtMaxHist[indexMCKKpi]->Fill(ptmax);
1293  fPtCandHist[indexMCKKpi]->Fill(ptCand);
1294  fDCAHist[indexMCKKpi]->Fill(dca);
1295  fNormIPHist[indexMCKKpi]->Fill(normIP);
1296  fCosPiDsHist[indexMCKKpi]->Fill(cosPiDs);
1297  fCosPiKPhiHist[indexMCKKpi]->Fill(cosPiKPhi);
1298  fPtProng0Hist[indexMCKKpi]->Fill(pt0);
1299  fPtProng1Hist[indexMCKKpi]->Fill(pt1);
1300  fPtProng2Hist[indexMCKKpi]->Fill(pt2);
1301  }
1302  }
1303  if(ispiKK){
1304  Double_t massKK=d->InvMass2Prongs(1,2,321,321);
1305  Double_t massKp=d->InvMass2Prongs(0,1,211,321);
1306  fDalitz[index]->Fill(massKK,massKp);
1307  if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);
1308  if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);
1309 
1310 
1311  if(fReadMC && indexMCpiKK!=-1){
1312  fDalitz[indexMCpiKK]->Fill(massKK,massKp);
1313  if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);
1314  if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);
1315  fCosPHist[indexMCpiKK]->Fill(cosp);
1316  fCosPxyHist[indexMCpiKK]->Fill(cospxy);
1317  fDLenHist[indexMCpiKK]->Fill(dlen);
1318  fDLenxyHist[indexMCpiKK]->Fill(dlenxy);
1319  fNDLenxyHist[indexMCpiKK]->Fill(normdlxy);
1320  fSigVertHist[indexMCpiKK]->Fill(sigvert);
1321  fSumd02Hist[indexMCpiKK]->Fill(sumD02);
1322  fPtMaxHist[indexMCpiKK]->Fill(ptmax);
1323  fPtCandHist[indexMCpiKK]->Fill(ptCand);
1324  fDCAHist[indexMCpiKK]->Fill(dca);
1325  fNormIPHist[indexMCpiKK]->Fill(normIP);
1326  fCosPiDsHist[indexMCpiKK]->Fill(cosPiDs);
1327  fCosPiKPhiHist[indexMCpiKK]->Fill(cosPiKPhi);
1328  fPtProng0Hist[indexMCpiKK]->Fill(pt0);
1329  fPtProng1Hist[indexMCpiKK]->Fill(pt1);
1330  fPtProng2Hist[indexMCpiKK]->Fill(pt2);
1331  }
1332  }
1333  }
1334 
1335  Float_t tmp[37];
1336  if(fFillNtuple>0 && fFillNtuple<4){
1337 
1338  if ((fFillNtuple==1 && (isPhiKKpi || isPhipiKK)) || (fFillNtuple==2 && (isK0starKKpi || isK0starpiKK)) || (fFillNtuple==3 && (isKKpi || ispiKK))){
1339 
1340  AliAODTrack *track0=(AliAODTrack*)d->GetDaughter(0);
1341  AliAODTrack *track1=(AliAODTrack*)d->GetDaughter(1);
1342  AliAODTrack *track2=(AliAODTrack*)d->GetDaughter(2);
1343 
1344  UInt_t bitMapPIDTrack0=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track0);
1345  UInt_t bitMapPIDTrack1=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track1);
1346  UInt_t bitMapPIDTrack2=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track2);
1347 
1348  tmp[0]=Float_t(labDs);
1349  if(fReadMC && fWriteOnlySignal) tmp[0]=Float_t(isMCSignal);
1350  tmp[1]=Float_t(retCodeAnalysisCuts);
1351  tmp[2]=Float_t(pdgCode0);
1352  tmp[3]=d->PtProng(0);
1353  tmp[4]=d->PtProng(1);
1354  tmp[5]=d->PtProng(2);
1355  tmp[6]=d->Pt();
1356  tmp[7]=d->PProng(0);
1357  tmp[8]=d->PProng(1);
1358  tmp[9]=d->PProng(2);
1359  tmp[10]=Int_t(bitMapPIDTrack0);
1360  tmp[11]=Int_t(bitMapPIDTrack1);
1361  tmp[12]=Int_t(bitMapPIDTrack2);
1362  tmp[13]=d->CosPointingAngle();
1363  tmp[14]=d->CosPointingAngleXY();
1364  tmp[15]=d->DecayLength();
1365  tmp[16]=d->DecayLengthXY();
1366  tmp[17]=d->NormalizedDecayLength();
1367  tmp[18]=d->NormalizedDecayLengthXY();
1368  tmp[19]=d->InvMassDsKKpi();
1369  tmp[20]=d->InvMassDspiKK();
1370  tmp[21]=d->GetSigmaVert();
1371  tmp[22]=d->Getd0Prong(0);
1372  tmp[23]=d->Getd0Prong(1);
1373  tmp[24]=d->Getd0Prong(2);
1374  tmp[25]=d->GetDCA();
1375  tmp[26]=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
1376  tmp[27]=d->InvMass2Prongs(0,1,321,321);
1377  tmp[28]=d->InvMass2Prongs(1,2,321,321);
1378  tmp[29]=d->InvMass2Prongs(1,2,321,211);
1379  tmp[30]=d->InvMass2Prongs(0,1,211,321);
1380  tmp[31]=d->CosPiDsLabFrameKKpi();
1381  tmp[32]=d->CosPiDsLabFramepiKK();
1382  tmp[33]=d->CosPiKPhiRFrameKKpi();
1383  tmp[34]=d->CosPiKPhiRFramepiKK();
1384  tmp[35]=(Float_t)(centrality);
1385  tmp[36]=(Float_t)(runNumber);
1386 
1387  if(fReadMC && fWriteOnlySignal){
1388  if(isMCSignal>=0) fNtupleDs->Fill(tmp);
1389  }else{
1390  fNtupleDs->Fill(tmp);
1391  }
1392  PostData(4,fNtupleDs);
1393  }
1394  }
1395  } //if(retCodeAnalysisCuts
1396  } // if(isFidAcc)
1397 
1398  if(unsetvtx) d->UnsetOwnPrimaryVtx();
1399  if(recVtx)fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
1400  }
1401 
1402  fCounter->StoreCandidates(aod,nFiltered,kTRUE);
1403  fCounter->StoreCandidates(aod,nSelected,kFALSE);
1404 
1405  delete vHF;
1406 
1407  PostData(1,fOutput);
1408  PostData(3,fCounter);
1409 
1410  return;
1411 }
1412 
1413 
1414 //_________________________________________________________________
1415 
1417 {
1419  //
1420  if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
1421  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1422  if (!fOutput) {
1423  printf("ERROR: fOutput not available\n");
1424  return;
1425  }
1426  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
1427  if(fHistNEvents){
1428  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1429  }else{
1430  printf("ERROR: fHistNEvents not available\n");
1431  return;
1432  }
1433  return;
1434 }
1435 
1436 //_________________________________________________________________
1437 void AliAnalysisTaskSEDs::FillMCGenAccHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader, Double_t nTracklets){
1441 
1442  Int_t nProng = 3;
1443  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
1444  if(TMath::Abs(zMCVertex) <= fAnalysisCuts->GetMaxVtxZ()) {
1445  for(Int_t iPart=0; iPart<arrayMC->GetEntriesFast(); iPart++){
1446 
1447  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
1448 
1449  if(TMath::Abs(mcPart->GetPdgCode()) == 431) {
1450  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,kTRUE);//Prompt = 4, FeedDown = 5
1451 
1452  Int_t deca = 0;
1453  Bool_t isGoodDecay = kFALSE;
1454  Int_t labDau[3] = {-1,-1,-1};
1455  Bool_t isFidAcc = kFALSE;
1456  Bool_t isDaugInAcc = kFALSE;
1457 
1458  deca = AliVertexingHFUtils::CheckDsDecay(arrayMC,mcPart,labDau);
1459  if(deca == 1) isGoodDecay=kTRUE; // == 1 -> Phi pi -> kkpi
1460 
1461  if(labDau[0]==-1) continue; //protection against unfilled array of labels
1462 
1463  if(isGoodDecay) {
1464  Double_t pt = mcPart->Pt();
1465  Double_t rapid = mcPart->Y();
1466  isFidAcc = fAnalysisCuts->IsInFiducialAcceptance(pt,rapid);
1467  isDaugInAcc = CheckDaugAcc(arrayMC,nProng,labDau);
1468 
1469  if(isDaugInAcc) {
1470  if (fFillAcceptanceLevel || (!fFillAcceptanceLevel && isFidAcc)) {
1471  Double_t var4nSparseAcc[knVarForSparseAcc] = {pt,rapid*10,nTracklets};
1472  Double_t ptWeight = 1.;
1473  if (fUseWeight && fHistoPtWeight){
1474  AliDebug(2,"Using Histogram as Pt weight function");
1475  ptWeight = GetPtWeightFromHistogram(pt);
1476  }
1477  if(orig==4) fnSparseMC[0]->Fill(var4nSparseAcc,ptWeight);
1478  if(orig==5) fnSparseMC[1]->Fill(var4nSparseAcc,ptWeight);
1479  }
1480  }
1481  }
1482  }
1483 
1484  if(fFillSparseDplus && TMath::Abs(mcPart->GetPdgCode()) == 411) {
1485  Int_t orig = AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,kTRUE);//Prompt = 4, FeedDown = 5
1486 
1487  Int_t deca = 0;
1488  Bool_t isGoodDecay = kFALSE;
1489  Int_t labDau[3] = {-1,-1,-1};
1490  Bool_t isFidAcc = kFALSE;
1491  Bool_t isDaugInAcc = kFALSE;
1492 
1493  deca = AliVertexingHFUtils::CheckDplusKKpiDecay(arrayMC,mcPart,labDau);
1494  if(deca == 1) isGoodDecay=kTRUE; // == 1 -> Phi pi -> kkpi
1495 
1496  if(labDau[0]==-1) continue; //protection against unfilled array of labels
1497 
1498  if(isGoodDecay) {
1499  Double_t pt = mcPart->Pt();
1500  Double_t rapid = mcPart->Y();
1501  isFidAcc = fAnalysisCuts->IsInFiducialAcceptance(pt,rapid);
1502  isDaugInAcc = CheckDaugAcc(arrayMC,nProng,labDau);
1503 
1504  if(isDaugInAcc) {
1505  if (fFillAcceptanceLevel || (!fFillAcceptanceLevel && isFidAcc)) {
1506  Double_t var4nSparseAcc[knVarForSparseAcc] = {pt,rapid*10,nTracklets};
1507  Double_t ptWeight = 1.;
1508  if (fUseWeight && fHistoPtWeight){
1509  AliDebug(2,"Using Histogram as Pt weight function");
1510  ptWeight = GetPtWeightFromHistogram(pt);
1511  }
1512  if(orig==4) fnSparseMCDplus[0]->Fill(var4nSparseAcc,ptWeight);
1513  if(orig==5) fnSparseMCDplus[1]->Fill(var4nSparseAcc,ptWeight);
1514  }
1515  }
1516  }
1517  }
1518  }
1519  }
1520 }
1521 
1522 //_________________________________________________________________
1523 Bool_t AliAnalysisTaskSEDs::CheckDaugAcc(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1525 
1526  for (Int_t iProng = 0; iProng<nProng; iProng++){
1527  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1528  if(!mcPartDaughter) {
1529  return kFALSE;
1530  }
1531  Double_t eta = mcPartDaughter->Eta();
1532  Double_t pt = mcPartDaughter->Pt();
1533  if (TMath::Abs(eta) > 0.9 || pt < 0.1) {
1534  return kFALSE;
1535  }
1536  }
1537  return kTRUE;
1538 }
1539 
1540 //_________________________________________________________________
1542 
1543  const Int_t nprongs = 3;
1544  Double_t PxProng[nprongs], PyProng[nprongs], PzProng[nprongs], P2Prong[nprongs], mProng[nprongs];
1545  Double_t Px, Py, Pz, P2;
1546  UInt_t pdg[3]={321,321,211};
1547  int idPion = 2;
1548  if(dec==2) {
1549  pdg[0]=211;
1550  pdg[2]=321;
1551  idPion = 0;
1552  }
1553 
1554  for (Int_t ip=0; ip<nprongs; ip++) {
1555  PxProng[ip] = d->PxProng(ip);
1556  PyProng[ip] = d->PxProng(ip);
1557  PzProng[ip] = d->PzProng(ip);
1558  P2Prong[ip] = d->P2Prong(ip);
1559  mProng[ip] = TDatabasePDG::Instance()->GetParticle(pdg[ip])->Mass();
1560  }
1561 
1562  for(Int_t i=0; i<9; i++) { //9 rotations implemented for the pion track around Pi
1563  Px = 0.;
1564  Py = 0.;
1565  Pz = 0.;
1566 
1567  Double_t phirot=TMath::Pi()*(5/6.+1/27.*i);
1568 
1569  PxProng[idPion] = PxProng[idPion]*TMath::Cos(phirot)-PyProng[idPion]*TMath::Sin(phirot);
1570  PyProng[idPion] = PxProng[idPion]*TMath::Sin(phirot)+PyProng[idPion]*TMath::Cos(phirot);
1571 
1572  for (Int_t j=0; j<nprongs; j++) {
1573  Px += PxProng[j];
1574  Py += PyProng[j];
1575  Pz += PzProng[j];
1576  }
1577  P2 = Px*Px + Py*Py + Pz*Pz;
1578 
1579  Double_t energysum = 0.;
1580  for(Int_t j=0; j<nprongs; j++) {
1581  energysum += TMath::Sqrt(mProng[j]*mProng[j]+P2Prong[j]);
1582  }
1583  Double_t mass = TMath::Sqrt(energysum*energysum-P2);
1584  if( (fminMass<=mass) && (mass<fmaxMass) ) fMassRotBkgHistPhi[iPtBin]->Fill(mass);
1585  }
1586 }
1587 
1588 //_________________________________________________________________________
1590 
1591  Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1592  Int_t nInvMassBins=(Int_t)(0.7/fMassBinSize+0.5);
1593  Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
1594  Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
1595 
1596  Int_t nTrklBins = 1;
1597  if(fUseTrkl) nTrklBins = 300;
1598  Int_t nCentrBins = 1;
1599  if(fUseCentrAxis) nCentrBins = 101;
1600 
1601  Int_t nBinsReco[knVarForSparse] = {nInvMassBins,24, 30, 14, 14, 20, 10, 10, 14, 6, 6, 12, nTrklBins, nCentrBins, 30};
1602  Double_t xminReco[knVarForSparse] = {minMass, 0., 0., 0., 0., 0., 90., 90., 0., 7., 0., 0., 1., 0., 0.};
1603  Double_t xmaxReco[knVarForSparse] = {maxMass, 24., 15., 70., 70., 10., 100., 100., 70., 10., 3., 6., 301., 101., 300.};
1604  TString axis[knVarForSparse] = {"invMassDsAllPhi","p_{T}","#Delta Mass(KK)","dlen","dlen_{xy}","normdl_{xy}","cosP","cosP_{xy}","sigVert","cosPiDs","|cosPiKPhi^{3}|","normIP","N tracklets",Form("Percentile (%s)",fCentEstName.Data()),"ImpPar_{xy}"};
1605 
1606  if(fSystem == 1) { //pPb,PbPb
1607  nInvMassBins=(Int_t)(0.45/fMassBinSize+0.5);
1608  minMass=massDs-0.5*nInvMassBins*fMassBinSize;
1609  maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
1610  nBinsReco[0] = nInvMassBins; //Ds mass
1611  xminReco[0] = minMass;
1612  xmaxReco[0] = maxMass;
1613 
1614  nBinsReco[1] = 16; //pt
1615  xminReco[1] = 0.;
1616  xmaxReco[1] = 16.;
1617 
1618  nBinsReco[2] = 12; //#Delta Mass(KK)
1619  xmaxReco[2] = 12.;
1620 
1621  nBinsReco[3] = 7; //dlen
1622  nBinsReco[4] = 7; //dlenxy
1623  nBinsReco[5] = 10; //ndlenxy
1624 
1625  nBinsReco[6] = 6; //cosP
1626  xminReco[6] = 97.;
1627  xmaxReco[6] = 100.;
1628 
1629  nBinsReco[7] = 6; //cosPxy
1630  xminReco[7] = 97.;
1631  xmaxReco[7] = 100.;
1632  }
1633 
1634  Int_t nBinsAcc[knVarForSparseAcc] = {20, 20, nTrklBins};
1635  Double_t xminAcc[knVarForSparseAcc] = {0., -10., 1.};
1636  Double_t xmaxAcc[knVarForSparseAcc] = {20, 10., 301.};
1637 
1638  if(fReadMC) {
1639  TString label[2] = {"fromC","fromB"};
1640  for (Int_t i=0; i<2; i++) {
1641  TString titleSparse = Form("MC nSparse (%s)- %s", fFillAcceptanceLevel ? "Acc.Step" : "Gen.Acc.Step", label[i].Data());
1642  fnSparseMC[i] = new THnSparseF(Form("fnSparseAcc_%s",label[i].Data()), titleSparse.Data(), knVarForSparseAcc, nBinsAcc, xminAcc, xmaxAcc);
1643  fnSparseMC[i]->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1644  fnSparseMC[i]->GetAxis(1)->SetTitle("y");
1645  fnSparseMC[i]->GetAxis(2)->SetTitle("N tracklets");
1646  fOutput->Add(fnSparseMC[i]);
1647 
1648  //Dplus
1649  if(fFillSparseDplus) {
1650  titleSparse = Form("MC nSparse D^{+} (%s)- %s", fFillAcceptanceLevel ? "Acc.Step" : "Gen.Acc.Step", label[i].Data());
1651  fnSparseMCDplus[i] = new THnSparseF(Form("fnSparseAccDplus_%s",label[i].Data()), titleSparse.Data(), knVarForSparseAcc, nBinsAcc, xminAcc, xmaxAcc);
1652  fnSparseMCDplus[i]->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1653  fnSparseMCDplus[i]->GetAxis(1)->SetTitle("y");
1654  fnSparseMCDplus[i]->GetAxis(2)->SetTitle("N tracklets");
1655  fOutput->Add(fnSparseMCDplus[i]);
1656  }
1657  }
1658  for (Int_t i=2; i<4; i++) {
1659  fnSparseMC[i] = new THnSparseF(Form("fnSparseReco_%s",label[i-2].Data()),Form("MC nSparse (Reco Step)- %s",label[i-2].Data()), knVarForSparse, nBinsReco, xminReco, xmaxReco);
1660  for (Int_t j=0; j<knVarForSparse; j++) {
1661  fnSparseMC[i]->GetAxis(j)->SetTitle(Form("%s",axis[j].Data()));
1662  }
1663  fOutput->Add(fnSparseMC[i]);
1664 
1665  //Dplus
1666  if(fFillSparseDplus) {
1667  fnSparseMCDplus[i] = new THnSparseF(Form("fnSparseRecoDplus_%s",label[i-2].Data()),Form("MC nSparse D^{+} (Reco Step)- %s",label[i-2].Data()), knVarForSparse, nBinsReco, xminReco, xmaxReco);
1668  for (Int_t j=0; j<knVarForSparse; j++) {
1669  fnSparseMCDplus[i]->GetAxis(j)->SetTitle(Form("%s",axis[j].Data()));
1670  }
1671  fOutput->Add(fnSparseMCDplus[i]);
1672  }
1673  }
1674  } //end MC
1675  else {
1676  fnSparse = new THnSparseF("fnSparse","nSparse", knVarForSparse, nBinsReco, xminReco, xmaxReco);
1677  for (Int_t j=0; j<knVarForSparse; j++) {
1678  fnSparse->GetAxis(j)->SetTitle(Form("%s",axis[j].Data()));
1679  }
1680  fOutput->Add(fnSparse);
1681  }
1682 }
1683 
1684 //_________________________________________________________________________
1686 
1687  Int_t nBinsIP[knVarForSparseIP] = { 20, 400, 400, 400, 400, 3};
1688  Double_t xminIP[knVarForSparseIP] = { 0., -10., -10., -10., -10., 0.};
1689  Double_t xmaxIP[knVarForSparseIP] = {20., 10., 10., 10., 10., 3.};
1690  TString axisIP[knVarForSparseIP] = {"motherPt","maxNormImp","IP0","IP1","IP2","candType"};
1691 
1692  fnSparseIP = new THnSparseF("fnSparseIP","nSparseIP", knVarForSparseIP, nBinsIP, xminIP, xmaxIP);
1693  for (Int_t j=0; j<knVarForSparseIP; j++) {
1694  fnSparseIP->GetAxis(j)->SetTitle(Form("%s",axisIP[j].Data()));
1695  }
1696  fnSparseIP->GetAxis(5)->SetTitle("candType (0.5=bkg; 1.5=prompt; 2.5=FD)");
1697  fOutput->Add(fnSparseIP);
1698 }
1699 
1700 //_________________________________________________________________________
1702 
1704  Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
1705  Int_t nInvMassBins=(Int_t)(0.7/fMassBinSize+0.5);
1706  Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
1707  Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
1708 
1709  Int_t nptbins=48;
1710  Double_t ptmin=0.;
1711  Double_t ptmax=24.;
1712 
1713  //dimensions for THnSparse
1714  TString axTit[kVarForImpPar]={"M_{K#pi#pi} (GeV/c^{2})","p_{T} (GeV/c)","Imp Par (#mum)"};
1715 
1716  Int_t nbins[kVarForImpPar]={nInvMassBins,nptbins,1000};
1717  Double_t xmin[kVarForImpPar]={minMass,ptmin,-1000};
1718  Double_t xmax[kVarForImpPar]={maxMass,ptmax,1000};
1719 
1720  //mass, pt, imppar
1721  fImpParSparse=new THnSparseF("hMassPtImpParAll","Mass vs. pt vs. imppar - All",kVarForImpPar,nbins,xmin,xmax);
1722  fImpParSparseMC[0]=new THnSparseF("hMassPtImpParPrompt","Mass vs. pt vs. imppar - promptD",kVarForImpPar,nbins,xmin,xmax);
1723  fImpParSparseMC[1]=new THnSparseF("hMassPtImpParBfeed","Mass vs. pt vs. imppar - DfromB",kVarForImpPar,nbins,xmin,xmax);
1724  fImpParSparseMC[2]=new THnSparseF("hMassPtImpParTrueBfeed","Mass vs. pt vs. true imppar -DfromB",kVarForImpPar,nbins,xmin,xmax);
1725  fImpParSparseMC[3]=new THnSparseF("hMassPtImpParBkg","Mass vs. pt vs. imppar - backgr.",kVarForImpPar,nbins,xmin,xmax);
1726 
1727  if(!fReadMC) fOutput->Add(fImpParSparse);
1728  else {
1729  for(Int_t iSparse=0; iSparse<4; iSparse++) {
1730  fOutput->Add(fImpParSparseMC[iSparse]);
1731  }
1732  }
1733 }
1734 
1735 //_________________________________________________________________________________________________
1736 Float_t AliAnalysisTaskSEDs::GetTrueImpactParameterDstoPhiPi(const AliAODMCHeader *mcHeader, TClonesArray* arrayMC, const AliAODMCParticle *partDs) const {
1738 
1739  Double_t vtxTrue[3];
1740  mcHeader->GetVertex(vtxTrue);
1741  Double_t origD[3];
1742  partDs->XvYvZv(origD);
1743  Short_t charge=partDs->Charge();
1744  Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1745  for(Int_t iDau=0; iDau<3; iDau++){
1746  pXdauTrue[iDau]=0.;
1747  pYdauTrue[iDau]=0.;
1748  pZdauTrue[iDau]=0.;
1749  }
1750 
1751  Int_t nDau=partDs->GetNDaughters();
1752  Int_t labelFirstDau = partDs->GetDaughter(0);
1753  if(nDau==2){
1754  Int_t theDau=0;
1755  for(Int_t iDau=0; iDau<2; iDau++){
1756  Int_t ind = labelFirstDau+iDau;
1757  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1758  if(!part){
1759  AliError("Daughter particle not found in MC array");
1760  return 99999.;
1761  }
1762  Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1763  if(pdgCode==211){
1764  pXdauTrue[theDau]=part->Px();
1765  pYdauTrue[theDau]=part->Py();
1766  pZdauTrue[theDau]=part->Pz();
1767  ++theDau;
1768  }else{
1769  Int_t nDauRes=part->GetNDaughters();
1770  if(nDauRes==2){
1771  Int_t labelFirstDauRes = part->GetDaughter(0);
1772  for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
1773  Int_t indDR = labelFirstDauRes+iDauRes;
1774  AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
1775  if(!partDR){
1776  AliError("Daughter particle not found in MC array");
1777  return 99999.;
1778  }
1779 
1780  Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
1781  if(pdgCodeDR==321){
1782  pXdauTrue[theDau]=partDR->Px();
1783  pYdauTrue[theDau]=partDR->Py();
1784  pZdauTrue[theDau]=partDR->Pz();
1785  ++theDau;
1786  }
1787  }
1788  }
1789  }
1790  }
1791  }
1792  else {
1793  AliError("Wrong number of decay prongs");
1794  return 99999.;
1795  }
1796 
1797  Double_t d0dummy[3]={0.,0.,0.};
1798  AliAODRecoDecayHF aodDsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1799  return aodDsMC.ImpParXY();
1800 }
1801 
1802 //_________________________________________________________________________
1804  // weight function from the ratio of the LHC16i2a MC
1805  // 1.5-14 GeV/c using data and 1-1.5, 14-50 GeV/c using FONLL calculations
1806 
1807  if(fHistoPtWeight) delete fHistoPtWeight;
1808  fHistoPtWeight = new TH1F("histoWeight","histoWeight",500,0.,50.);
1809  fHistoPtWeight->Sumw2();
1810  Float_t binc[500]={ 1.695705, 1.743693, 1.790289, 1.835410, 1.878981, 1.920938, 1.961223, 1.999787, 2.036589, 2.071597, 2.104784, 2.136132, 2.165629, 2.193270, 2.219057, 2.174545, 2.064698, 1.959489, 1.858770, 1.762396, 1.670224, 1.582115, 1.497931, 1.417541, 1.340814, 1.267622, 1.197842, 1.131352, 1.068033, 1.007770, 0.950450, 0.895963, 0.844202, 0.795062, 0.748441, 0.704241, 0.662363, 0.622715, 0.585204, 0.549742, 0.516242, 0.484620, 0.454795, 0.426686, 0.400217, 0.375314, 0.351903, 0.329915, 0.309281, 0.289936, 0.271816, 0.254860, 0.239007, 0.224201, 0.210386, 0.197508, 0.185516, 0.174360, 0.163992, 0.154366, 0.145438, 0.137166, 0.129508, 0.122426, 0.115882, 0.109840, 0.104266, 0.099128, 0.094395, 0.090036, 0.086023, 0.082331, 0.078933, 0.075805, 0.072925, 0.070271, 0.067823, 0.065562, 0.063471, 0.061532, 0.059730, 0.058051, 0.056481, 0.055007, 0.053619, 0.052306, 0.051059, 0.049867, 0.048725, 0.047624, 0.046558, 0.045522, 0.044511, 0.043521, 0.042548, 0.041590, 0.040643, 0.039706, 0.038778, 0.037857, 0.036944, 0.036039, 0.035141, 0.034251, 0.033370, 0.032500, 0.031641, 0.030796, 0.029966, 0.029153, 0.028359, 0.027587, 0.026837, 0.026113, 0.025416, 0.024748, 0.024111, 0.023507, 0.022937, 0.022402, 0.021904, 0.021443, 0.021020, 0.020634, 0.020286, 0.019974, 0.019698, 0.019455, 0.019244, 0.019062, 0.018905, 0.018770, 0.018652, 0.018545, 0.018444, 0.018342, 0.018231, 0.018102, 0.017947, 0.017755, 0.017536, 0.017327, 0.017120, 0.016915, 0.016713, 0.016514, 0.016317, 0.016122, 0.015929, 0.015739, 0.015551, 0.015366, 0.015182, 0.015001, 0.014822, 0.014645, 0.014470, 0.014297, 0.014127, 0.013958, 0.013791, 0.013627, 0.013464, 0.013303, 0.013145, 0.012988, 0.012833, 0.012679, 0.012528, 0.012378, 0.012231, 0.012085, 0.011940, 0.011798, 0.011657, 0.011518, 0.011380, 0.011244, 0.011110, 0.010978, 0.010846, 0.010717, 0.010589, 0.010463, 0.010338, 0.010214, 0.010092, 0.009972, 0.009853, 0.009735, 0.009619, 0.009504, 0.009391, 0.009279, 0.009168, 0.009058, 0.008950, 0.008843, 0.008738, 0.008633, 0.008530, 0.008429, 0.008328, 0.008229, 0.008130, 0.008033, 0.007937, 0.007843, 0.007749, 0.007656, 0.007565, 0.007475, 0.007385, 0.007297, 0.007210, 0.007124, 0.007039, 0.006955, 0.006872, 0.006790, 0.006709, 0.006629, 0.006550, 0.006471, 0.006394, 0.006318, 0.006242, 0.006168, 0.006094, 0.006022, 0.005950, 0.005879, 0.005808, 0.005739, 0.005671, 0.005603, 0.005536, 0.005470, 0.005405, 0.005340, 0.005276, 0.005213, 0.005151, 0.005090, 0.005029, 0.004969, 0.004909, 0.004851, 0.004793, 0.004736, 0.004679, 0.004623, 0.004568, 0.004514, 0.004460, 0.004406, 0.004354, 0.004302, 0.004251, 0.004200, 0.004150, 0.004100, 0.004051, 0.004003, 0.003955, 0.003908, 0.003861, 0.003815, 0.003770, 0.003725, 0.003680, 0.003636, 0.003593, 0.003550, 0.003507, 0.003466, 0.003424, 0.003383, 0.003343, 0.003303, 0.003264, 0.003225, 0.003186, 0.003148, 0.003110, 0.003073, 0.003037, 0.003000, 0.002965, 0.002929, 0.002894, 0.002860, 0.002826, 0.002792, 0.002758, 0.002726, 0.002693, 0.002661, 0.002629, 0.002598, 0.002567, 0.002536, 0.002506, 0.002476, 0.002446, 0.002417, 0.002388, 0.002360, 0.002332, 0.002304, 0.002276, 0.002249, 0.002222, 0.002196, 0.002169, 0.002144, 0.002118, 0.002093, 0.002068, 0.002043, 0.002019, 0.001995, 0.001971, 0.001947, 0.001924, 0.001901, 0.001878, 0.001856, 0.001834, 0.001812, 0.001790, 0.001769, 0.001748, 0.001727, 0.001706, 0.001686, 0.001666, 0.001646, 0.001626, 0.001607, 0.001588, 0.001569, 0.001550, 0.001531, 0.001513, 0.001495, 0.001477, 0.001460, 0.001442, 0.001425, 0.001408, 0.001391, 0.001374, 0.001358, 0.001342, 0.001326, 0.001310, 0.001294, 0.001279, 0.001264, 0.001249, 0.001234, 0.001219, 0.001204, 0.001190, 0.001176, 0.001162, 0.001148, 0.001134, 0.001121, 0.001107, 0.001094, 0.001081, 0.001068, 0.001055, 0.001043, 0.001030, 0.001018, 0.001006, 0.000994, 0.000982, 0.000970, 0.000959, 0.000947, 0.000936, 0.000925, 0.000914, 0.000903, 0.000892, 0.000881, 0.000871, 0.000860, 0.000850, 0.000840, 0.000830, 0.000820, 0.000810, 0.000801, 0.000791, 0.000782, 0.000772, 0.000763, 0.000754, 0.000745, 0.000736, 0.000727, 0.000719, 0.000710, 0.000702, 0.000693, 0.000685, 0.000677, 0.000669, 0.000661, 0.000653, 0.000645, 0.000637, 0.000630, 0.000622, 0.000615, 0.000607, 0.000600, 0.000593, 0.000586, 0.000579, 0.000572, 0.000565, 0.000558, 0.000552, 0.000545, 0.000539, 0.000532, 0.000526, 0.000520, 0.000513, 0.000507, 0.000501, 0.000495, 0.000489, 0.000483, 0.000478, 0.000472, 0.000466, 0.000461, 0.000455, 0.000450, 0.000444, 0.000439, 0.000434, 0.000429, 0.000424, 0.000419, 0.000414, 0.000409, 0.000404, 0.000399, 0.000394, 0.000389, 0.000385, 0.000380, 0.000376, 0.000371, 0.000367, 0.000362, 0.000358, 0.000354, 0.000350, 0.000345, 0.000341, 0.000337, 0.000333, 0.000329, 0.000325, 0.000321, 0.000318, 0.000314, 0.000310, 0.000306, 0.000303, 0.000299, 0.000295, 0.000292, 0.000288, 0.000285, 0.000282, 0.000278, 0.000275, 0.000272, 0.000268, 0.000265, 0.000262, 0.000259, 0.000256, 0.000253, 0.000250, 0.000247, 0.000244, 0.000241, 0.000238, 0.000235};
1811  for(Int_t i=0; i<500; i++){
1812  fHistoPtWeight->SetBinContent(i+1,binc[i]);
1813  }
1814  //SetWeightHistogram();
1815  fUseWeight=kTRUE;
1816 }
1817 
1818 //_________________________________________________________________________
1820  // weight function from the ratio of the LHC16i2a+b+c MC
1821  // and FONLL calculations for pp data
1822 
1823  if(fHistoPtWeight) delete fHistoPtWeight;
1824  fHistoPtWeight = new TH1F("histoWeight","histoWeight",400,0.,40.);
1825  fHistoPtWeight->Sumw2();
1826  Float_t binc[400]={1.118416, 1.003458, 0.935514, 0.907222, 0.904359, 0.913668, 0.933906, 0.963898, 0.996388, 1.031708, 1.066404, 1.099683, 1.125805, 1.145181, 1.165910, 1.181905, 1.193425, 1.203891, 1.204726, 1.209411, 1.209943, 1.204763, 1.205291, 1.198912, 1.197390, 1.182005, 1.184194, 1.175994, 1.167881, 1.158348, 1.147190, 1.139833, 1.126940, 1.123322, 1.108389, 1.102199, 1.089464, 1.075874, 1.061964, 1.051429, 1.038113, 1.026668, 1.011441, 0.998567, 0.987658, 0.972434, 0.950068, 0.940758, 0.916880, 0.911931, 0.894512, 0.878691, 0.860589, 0.848025, 0.830774, 0.819399, 0.801134, 0.775276, 0.766382, 0.750495, 0.736935, 0.717529, 0.702637, 0.689152, 0.671334, 0.652030, 0.635696, 0.621365, 0.608362, 0.599019, 0.576024, 0.562136, 0.550938, 0.533587, 0.516410, 0.509744, 0.501655, 0.487402, 0.476469, 0.463762, 0.445979, 0.438088, 0.422214, 0.417467, 0.404357, 0.391450, 0.379996, 0.371201, 0.361497, 0.352912, 0.343189, 0.329183, 0.327662, 0.310783, 0.304525, 0.301007, 0.293306, 0.278332, 0.274419, 0.267361, 0.261459, 0.255514, 0.249293, 0.241129, 0.237600, 0.231343, 0.221982, 0.216872, 0.211094, 0.206954, 0.202333, 0.196572, 0.193274, 0.188240, 0.181817, 0.178364, 0.173614, 0.167135, 0.166055, 0.163423, 0.156557, 0.155821, 0.151985, 0.144909, 0.145062, 0.139720, 0.138873, 0.131892, 0.129969, 0.126509, 0.126978, 0.120451, 0.117661, 0.116300, 0.115604, 0.112215, 0.109237, 0.107720, 0.106419, 0.102050, 0.102777, 0.097406, 0.098447, 0.095964, 0.093868, 0.092430, 0.089329, 0.088249, 0.085881, 0.084417, 0.085498, 0.082444, 0.079151, 0.079565, 0.077811, 0.077293, 0.075218, 0.072445, 0.073054, 0.071545, 0.070279, 0.068046, 0.067854, 0.068092, 0.065378, 0.064405, 0.062060, 0.063391, 0.061718, 0.059616, 0.058913, 0.058895, 0.058311, 0.056320, 0.056527, 0.055349, 0.053701, 0.054735, 0.052264, 0.051277, 0.051554, 0.050545, 0.048995, 0.049507, 0.048466, 0.048156, 0.046809, 0.047600, 0.046078, 0.044801, 0.044113, 0.043700, 0.043530, 0.043396, 0.042556, 0.041048, 0.041657, 0.040394, 0.041314, 0.040720, 0.039656, 0.038478, 0.039276, 0.038777, 0.037730, 0.036918, 0.036466, 0.035827, 0.035285, 0.035963, 0.034371, 0.034757, 0.033205, 0.033666, 0.033266, 0.032583, 0.033570, 0.032102, 0.032107, 0.031464, 0.032160, 0.030091, 0.030564, 0.029464, 0.029613, 0.029626, 0.029512, 0.029324, 0.028607, 0.027628, 0.027251, 0.027072, 0.027077, 0.026724, 0.026961, 0.026303, 0.026237, 0.025454, 0.025133, 0.025365, 0.026014, 0.024807, 0.023901, 0.023459, 0.023405, 0.023654, 0.023981, 0.023675, 0.022493, 0.022781, 0.021801, 0.021704, 0.022372, 0.021189, 0.020681, 0.020779, 0.021324, 0.020558, 0.020901, 0.020586, 0.020808, 0.019276, 0.019516, 0.019706, 0.018935, 0.018632, 0.018516, 0.019187, 0.018916, 0.018039, 0.018208, 0.018045, 0.017628, 0.017916, 0.017711, 0.017838, 0.017222, 0.016565, 0.015733, 0.016264, 0.015826, 0.016090, 0.016622, 0.015802, 0.016621, 0.015441, 0.015309, 0.014860, 0.014935, 0.014968, 0.014443, 0.014485, 0.015136, 0.014078, 0.014414, 0.013908, 0.014071, 0.014078, 0.013766, 0.013436, 0.013507, 0.013480, 0.013224, 0.013041, 0.013935, 0.012885, 0.012453, 0.012528, 0.012492, 0.012225, 0.012542, 0.012706, 0.012136, 0.011902, 0.011560, 0.011448, 0.011861, 0.011271, 0.011831, 0.011159, 0.011171, 0.010966, 0.011311, 0.011002, 0.011130, 0.010995, 0.010450, 0.010663, 0.010678, 0.010492, 0.009861, 0.010507, 0.009916, 0.010121, 0.010029, 0.010046, 0.009370, 0.009647, 0.010104, 0.009282, 0.009830, 0.009403, 0.009148, 0.009172, 0.008893, 0.009158, 0.009019, 0.008780, 0.008579, 0.009063, 0.008634, 0.008988, 0.008265, 0.008581, 0.008575, 0.008690, 0.008181, 0.008352, 0.008150, 0.008430, 0.008256, 0.008119, 0.008453, 0.008447, 0.008021, 0.007938, 0.008025, 0.007718, 0.008127, 0.007651, 0.007590, 0.007316, 0.007839, 0.007504, 0.007341, 0.007527, 0.007263, 0.007668, 0.007306, 0.007271, 0.006910, 0.007257, 0.007260, 0.006810, 0.006967, 0.006887, 0.006867, 0.007202, 0.006829, 0.006370, 0.006710, 0.006417, 0.006361, 0.006800, 0.006410, 0.006323, 0.006790, 0.006322, 0.006673, 0.006547};
1827  for(Int_t i=0; i<400; i++){
1828  fHistoPtWeight->SetBinContent(i+1,binc[i]);
1829  }
1830  //SetWeightHistogram();
1831  fUseWeight=kTRUE;
1832 }
1833 
1834 //_________________________________________________________________________
1836  // weight function from the ratio of the LHC16i2a+b+c MC
1837  // and FONLL calculations for pp data
1838  // corrected by the BAMPS Raa calculation for 30-50% CC
1839  if(fHistoPtWeight) delete fHistoPtWeight;
1840  fHistoPtWeight = new TH1F("histoWeight","histoWeight",400,0.,40.);
1841  fHistoPtWeight->Sumw2();
1842  Float_t binc[400]={2.166180, 1.866117, 1.667595, 1.547176, 1.486661, 1.457891, 1.426949, 1.399055, 1.383278, 1.349383, 1.317009, 1.282321, 1.234257, 1.181136, 1.136655, 1.087523, 1.037912, 0.993256, 0.944746, 0.900948, 0.865869, 0.827193, 0.794424, 0.757723, 0.733020, 0.700164, 0.682189, 0.659872, 0.637918, 0.615749, 0.593020, 0.574402, 0.556158, 0.542663, 0.525494, 0.516038, 0.503629, 0.490980, 0.479143, 0.469005, 0.457749, 0.447668, 0.436803, 0.427073, 0.418282, 0.407867, 0.395093, 0.387861, 0.374742, 0.369462, 0.360146, 0.351991, 0.342990, 0.336259, 0.327730, 0.322382, 0.314602, 0.303874, 0.299820, 0.293049, 0.287539, 0.280329, 0.274866, 0.269939, 0.263299, 0.256057, 0.249215, 0.242170, 0.235704, 0.230709, 0.220529, 0.213921, 0.208394, 0.202424, 0.196700, 0.194943, 0.192620, 0.187894, 0.184411, 0.180204, 0.172915, 0.169077, 0.162201, 0.159636, 0.153904, 0.148296, 0.143282, 0.139306, 0.135561, 0.132342, 0.128696, 0.123444, 0.122873, 0.116544, 0.114197, 0.112878, 0.110018, 0.104547, 0.103222, 0.100707, 0.098622, 0.096513, 0.094295, 0.091334, 0.090122, 0.087870, 0.084894, 0.083729, 0.082265, 0.081404, 0.080323, 0.078750, 0.078132, 0.076781, 0.074823, 0.074050, 0.072614, 0.070093, 0.069828, 0.068907, 0.066189, 0.066054, 0.064600, 0.061757, 0.061986, 0.059862, 0.059656, 0.056807, 0.055956, 0.054386, 0.054507, 0.051629, 0.050358, 0.049702, 0.049331, 0.047814, 0.046476, 0.045762, 0.045142, 0.043224, 0.043484, 0.041282, 0.041794, 0.040809, 0.039985, 0.039439, 0.038181, 0.037782, 0.036831, 0.036264, 0.036790, 0.035535, 0.034173, 0.034409, 0.033659, 0.033308, 0.032290, 0.030981, 0.031121, 0.030361, 0.029708, 0.028653, 0.028461, 0.028449, 0.027208, 0.026697, 0.025623, 0.026069, 0.025279, 0.024332, 0.024341, 0.024629, 0.024677, 0.024117, 0.024490, 0.024257, 0.023804, 0.024537, 0.023692, 0.023502, 0.023888, 0.023673, 0.023193, 0.023684, 0.023429, 0.023521, 0.023014, 0.023346, 0.022544, 0.021866, 0.021477, 0.021224, 0.021089, 0.020972, 0.020515, 0.019739, 0.019982, 0.019328, 0.019719, 0.019387, 0.018833, 0.018227, 0.018558, 0.018276, 0.017738, 0.017460, 0.017365, 0.017178, 0.017033, 0.017478, 0.016817, 0.017119, 0.016463, 0.016802, 0.016711, 0.016475, 0.017083, 0.016441, 0.016548, 0.016320, 0.016786, 0.015804, 0.016153, 0.015668, 0.015843, 0.015810, 0.015651, 0.015454, 0.014981, 0.014376, 0.014089, 0.013906, 0.013818, 0.013549, 0.013580, 0.013160, 0.013040, 0.012566, 0.012324, 0.012353, 0.012582, 0.011915, 0.011401, 0.011112, 0.011008, 0.011046, 0.011119, 0.010954, 0.010439, 0.010604, 0.010179, 0.010163, 0.010507, 0.009981, 0.009771, 0.009846, 0.010134, 0.009798, 0.009991, 0.009869, 0.010005, 0.009295, 0.009438, 0.009557, 0.009210, 0.009088, 0.009057, 0.009412, 0.009306, 0.008899, 0.009009, 0.008952, 0.008764, 0.008926, 0.008842, 0.008924, 0.008634, 0.008322, 0.007920, 0.008205, 0.008000, 0.008151, 0.008438, 0.008037, 0.008472, 0.007886, 0.007835, 0.007621, 0.007675, 0.007707, 0.007452, 0.007489, 0.007841, 0.007308, 0.007497, 0.007248, 0.007348, 0.007367, 0.007227, 0.007097, 0.007179, 0.007209, 0.007115, 0.007059, 0.007588, 0.007058, 0.006862, 0.006945, 0.006965, 0.006856, 0.007075, 0.007209, 0.006925, 0.006830, 0.006672, 0.006645, 0.006923, 0.006615, 0.006982, 0.006622, 0.006666, 0.006579, 0.006823, 0.006673, 0.006786, 0.006740, 0.006440, 0.006606, 0.006650, 0.006568, 0.006206, 0.006646, 0.006305, 0.006468, 0.006442, 0.006486, 0.006080, 0.006291, 0.006622, 0.006113, 0.006506, 0.006254, 0.006114, 0.006161, 0.006002, 0.006211, 0.006146, 0.006012, 0.005902, 0.006264, 0.005996, 0.006271, 0.005793, 0.006043, 0.006067, 0.006177, 0.005842, 0.005991, 0.005872, 0.006102, 0.006003, 0.005930, 0.006201, 0.006224, 0.005937, 0.005901, 0.005992, 0.005788, 0.006121, 0.005787, 0.005766, 0.005582, 0.006006, 0.005774, 0.005672, 0.005841, 0.005660, 0.006000, 0.005741, 0.005737, 0.005475, 0.005773, 0.005799, 0.005462, 0.005610, 0.005569, 0.005574, 0.005871, 0.005589, 0.005234, 0.005535, 0.005314, 0.005288, 0.005676, 0.005371, 0.005319, 0.005734, 0.005360, 0.005679, 0.005593};
1843  for(Int_t i=0; i<400; i++){
1844  fHistoPtWeight->SetBinContent(i+1,binc[i]);
1845  }
1846  fUseWeight=kTRUE;
1847 }
1848 
1849 //_________________________________________________________________________
1851  // weight function from the ratio of the LHC16i2a+b+c MC
1852  // and FONLL calculations for pp data
1853  // corrected by the TAMU Raa calculation for 0-10% CC (not available in 30-50% CC)
1854  if(fHistoPtWeight) delete fHistoPtWeight;
1855  fHistoPtWeight = new TH1F("histoWeight","histoWeight",400,0.,40.);
1856  fHistoPtWeight->Sumw2();
1857  Float_t binc[400]={1.179906, 1.091249, 1.047774, 1.045579, 1.071679, 1.112413, 1.167414, 1.236240, 1.310301, 1.390289, 1.471711, 1.553389, 1.626886, 1.692115, 1.760647, 1.813658, 1.850817, 1.886699, 1.907671, 1.934832, 1.955433, 1.966727, 1.987262, 1.996316, 2.013326, 1.973926, 1.931144, 1.871654, 1.812942, 1.752718, 1.690846, 1.635303, 1.572611, 1.523510, 1.459790, 1.402510, 1.331908, 1.261575, 1.192241, 1.127915, 1.061798, 0.998830, 0.933514, 0.871774, 0.812936, 0.762844, 0.719340, 0.686587, 0.644108, 0.615714, 0.579512, 0.545254, 0.510508, 0.479884, 0.447423, 0.426154, 0.408934, 0.388264, 0.376424, 0.361389, 0.347757, 0.331685, 0.318029, 0.305285, 0.290922, 0.278523, 0.269807, 0.262025, 0.254878, 0.249325, 0.238179, 0.230899, 0.224792, 0.216253, 0.207879, 0.204465, 0.201153, 0.195373, 0.190926, 0.185773, 0.178589, 0.175371, 0.168959, 0.167004, 0.161705, 0.156809, 0.152788, 0.149806, 0.146429, 0.143478, 0.140037, 0.134813, 0.134679, 0.128205, 0.126078, 0.125038, 0.122214, 0.116329, 0.115044, 0.112427, 0.110279, 0.108098, 0.105784, 0.102628, 0.101429, 0.099101, 0.095464, 0.093631, 0.091491, 0.090045, 0.088374, 0.086188, 0.085067, 0.083168, 0.080636, 0.079414, 0.077610, 0.075013, 0.074825, 0.073932, 0.071106, 0.071050, 0.069574, 0.066593, 0.066924, 0.064876, 0.065064, 0.062345, 0.061980, 0.060859, 0.061616, 0.058952, 0.058079, 0.057894, 0.058031, 0.056604, 0.055180, 0.054490, 0.053909, 0.051768, 0.052210, 0.049552, 0.050152, 0.048955, 0.047953, 0.047224, 0.045588, 0.044985, 0.043728, 0.042934, 0.043434, 0.041834, 0.040118, 0.040281, 0.039348, 0.038987, 0.037793, 0.036258, 0.036420, 0.035528, 0.034761, 0.033524, 0.033296, 0.033280, 0.031825, 0.031351, 0.030329, 0.031103, 0.030401, 0.029481, 0.029247, 0.029352, 0.029174, 0.028286, 0.028500, 0.028017, 0.027293, 0.027932, 0.026779, 0.026379, 0.026628, 0.026211, 0.025508, 0.025877, 0.025433, 0.025328, 0.024636, 0.025069, 0.024282, 0.023625, 0.023278, 0.023074, 0.023000, 0.022943, 0.022514, 0.021767, 0.022180, 0.021594, 0.022175, 0.021944, 0.021456, 0.020901, 0.021419, 0.021230, 0.020738, 0.020322, 0.020055, 0.019686, 0.019371, 0.019725, 0.018835, 0.019029, 0.018163, 0.018398, 0.018163, 0.017719, 0.018126, 0.017208, 0.017086, 0.016622, 0.016865, 0.015663, 0.015791, 0.015108, 0.015069, 0.015033, 0.015006, 0.014940, 0.014604, 0.014133, 0.013968, 0.013904, 0.013934, 0.013780, 0.013930, 0.013727, 0.013940, 0.013763, 0.013826, 0.014192, 0.014801, 0.014347, 0.014048, 0.014009, 0.014197, 0.014571, 0.014999, 0.015030, 0.014491, 0.014891, 0.014456, 0.014596, 0.015256, 0.014648, 0.014492, 0.014756, 0.015344, 0.014986, 0.015433, 0.015394, 0.015756, 0.014778, 0.015145, 0.015478, 0.015051, 0.014986, 0.015067, 0.015793, 0.015748, 0.015188, 0.015502, 0.015533, 0.015340, 0.015759, 0.015745, 0.016026, 0.015635, 0.015194, 0.014579, 0.015225, 0.014963, 0.015365, 0.016030, 0.015387, 0.016341, 0.015327, 0.015340, 0.015030, 0.015246, 0.015420, 0.015015, 0.015195, 0.016021, 0.015034, 0.015528, 0.015114, 0.015423, 0.015564, 0.015348, 0.015107, 0.015314, 0.015411, 0.015243, 0.015154, 0.016324, 0.015215, 0.014823, 0.015030, 0.015104, 0.014896, 0.015400, 0.015721, 0.015131, 0.014951, 0.014630, 0.014597, 0.015235, 0.014583, 0.015418, 0.014648, 0.014769, 0.014601, 0.015167, 0.014857, 0.015134, 0.015053, 0.014405, 0.014800, 0.014921, 0.014760, 0.013966, 0.014979, 0.014230, 0.014620, 0.014581, 0.014701, 0.013799, 0.014299, 0.015071, 0.013931, 0.014846, 0.014290, 0.013988, 0.014113, 0.013767, 0.014263, 0.014131, 0.013840, 0.013604, 0.014456, 0.013853, 0.014505, 0.013416, 0.014010, 0.014081, 0.014352, 0.013589, 0.013952, 0.013690, 0.014241, 0.014024, 0.013868, 0.014517, 0.014587, 0.013927, 0.013857, 0.014084, 0.013619, 0.014417, 0.013644, 0.013607, 0.013185, 0.014200, 0.013665, 0.013437, 0.013849, 0.013431, 0.014252, 0.013648, 0.013652, 0.013039, 0.013761, 0.013836, 0.013043, 0.013408, 0.013319, 0.013344, 0.014065, 0.013400, 0.012560, 0.013294, 0.012773, 0.012721, 0.013663, 0.012939, 0.012823, 0.013835, 0.012942, 0.013723, 0.013525};
1858  for(Int_t i=0; i<400; i++){
1859  fHistoPtWeight->SetBinContent(i+1,binc[i]);
1860  }
1861  fUseWeight=kTRUE;
1862 }
1863 
1864 //_________________________________________________________________________
1866  // Weight function from the ratio of the LHC17c3a1+2 MC
1867  // and FONLL calculations for pp data at 13 TeV
1868  // (From D0, Susanna Costanza)
1869  if(fHistoPtWeight) delete fHistoPtWeight;
1870  fHistoPtWeight = new TH1F("histoWeight","histoWeight",400,0.,40.);
1871  fHistoPtWeight->Sumw2();
1872  Float_t binc[400]={1.489198, 1.386131, 1.328213, 1.309866, 1.324383, 1.364766, 1.424282, 1.496803, 1.576996, 1.660407, 1.743464, 1.823431, 1.898334, 1.966862, 2.028262, 2.082235, 2.128839, 2.168395, 2.201406, 2.228492, 2.250331, 2.267616, 2.281022, 2.291181, 2.298668, 2.303992, 2.307592, 2.309839, 2.311043, 2.311454, 2.311272, 2.310653, 2.309719, 2.308559, 2.307242, 2.305816, 2.304318, 2.302771, 2.301193, 2.299596, 2.297986, 2.296369, 2.294749, 2.293127, 2.291504, 2.289881, 2.288259, 2.286638, 2.285018, 2.283398, 2.281780, 2.280163, 2.278547, 2.276933, 2.275319, 2.273707, 2.272095, 2.270485, 2.268876, 2.267268, 2.265661, 2.264056, 2.262451, 2.260848, 2.259246, 2.257645, 2.256045, 2.254446, 2.252848, 2.251252, 2.249656, 2.248062, 2.246469, 2.244877, 2.243286, 2.241696, 2.240108, 2.238520, 2.236934, 2.235348, 2.233764, 2.232181, 2.230599, 2.229019, 2.227439, 2.225860, 2.224283, 2.222707, 2.221132, 2.219557, 2.217985, 2.216413, 2.214842, 2.213272, 2.211704, 2.210136, 2.208570, 2.207005, 2.205441, 2.203878, 2.202316, 2.200755, 2.199196, 2.197637, 2.196080, 2.194524, 2.192968, 2.192011, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000, 2.192000};
1873  for(Int_t i=0; i<400; i++){
1874  fHistoPtWeight->SetBinContent(i+1,binc[i]);
1875  }
1876  fUseWeight=kTRUE;
1877 }
1878 
1879 //_________________________________________________________________________
1881 {
1882  //
1883  // Using an histogram as weight function
1884  // weight = 0 in the range outside the function
1885  //
1886  Double_t weight = 0.0;
1887  Int_t histoNbins = fHistoPtWeight->GetNbinsX();
1888  Int_t bin2 = fHistoPtWeight->FindBin(pt);
1889  if( (bin2>0) && (bin2<=histoNbins) ) {
1890  Int_t bin1=bin2-1;
1891  Int_t bin3=bin2+1;
1892  if(bin2==1) bin1=bin2+2;
1893  if(bin2==histoNbins) bin3=bin2-2;
1894  Float_t x_1=fHistoPtWeight->GetXaxis()->GetBinCenter(bin1);
1895  Float_t x_2=fHistoPtWeight->GetXaxis()->GetBinCenter(bin2);
1896  Float_t x_3=fHistoPtWeight->GetXaxis()->GetBinCenter(bin3);
1897  Float_t y_1=fHistoPtWeight->GetBinContent(bin1);
1898  Float_t y_2=fHistoPtWeight->GetBinContent(bin2);
1899  Float_t y_3=fHistoPtWeight->GetBinContent(bin3);
1900  Double_t a=( (y_3-y_2)*(x_1-x_2) - (y_1-y_2)*(x_3-x_2) )/( (x_3*x_3-x_2*x_2)*(x_1-x_2) - (x_1*x_1-x_2*x_2)*(x_3-x_2) );
1901  Double_t b=((y_1-y_2)-a*(x_1*x_1-x_2*x_2))/(x_1-x_2);
1902  Double_t c=y_3-a*(x_3*x_3)-b*x_3;
1903  weight = a*pt*pt+b*pt+c;
1904  }
1905  return weight;
1906 }
1907 
1908 
1909 
1910 
1911 
1912 
1913 
Int_t charge
Double_t NormalizedDecayLengthXY() const
Bool_t IsEventRejectedDueToCentrality() const
Definition: AliRDHFCuts.h:330
Int_t pdg
THnSparseF * fnSparseMCDplus[4]
TH1F * fPtMaxHist[4 *kMaxPtBins]
! hist. for Pt Max (Prod Cuts)
Bool_t IsEventRejectedDueToZVertexOutsideFiducialRegion() const
Definition: AliRDHFCuts.h:324
Double_t NormalizedDecayLength() const
THnSparseF * fnSparseMC[4]
!<!THnSparse for topomatic variable
Bool_t IsEventRejectedDueToNotRecoVertex() const
Definition: AliRDHFCuts.h:318
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:321
TH2F * fPtVsMass
! hist. of pt vs. mass (prod. cuts)
void SetFillCentralityAxis(Int_t flag=0)
TCanvas * c
Definition: TestFitELoss.C:172
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
THnSparseF * fImpParSparseMC[4]
!<!THnSparse for imp. par. on data
Float_t GetTrueImpactParameterDstoPhiPi(const AliAODMCHeader *mcHeader, TClonesArray *arrayMC, const AliAODMCParticle *partDs) 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)
Functions to check the decay tree.
Bool_t fFillSparse
flag for usage of HasSelectionBit
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:256
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)
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
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
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 filling true reconstructed Ds at acceptance level (see FillMCGenAccHistos) ...
TH3F * fPtProng0Hist3D
! Pt prong0 vs Ds mass vs pt
TH1F * fHistCentrality[3]
!hist. for cent distr (all,sel ev, )
void SetPtWeightsFromFONLL5anddataoverLHC16i2a()
static Int_t GetNumberOfTrackletsInEtaRange(AliAODEvent *ev, Double_t mineta, Double_t maxeta)
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)
Double_t GetPtWeightFromHistogram(Double_t pt)
AliRDHFCutsDstoKKpi * fAnalysisCuts
TList * fOutput
! list send on output slot 0
TH2F * fHistAllV0multNTPCout
! histo for V0mult vs #tracks TPCout (all)
Double_t CosPiKPhiRFrameKKpi() const
Bool_t fUseTrkl
flag to decide whether to use pt-weights != 1 when filling the container or not
TH1F * fCosPxyHist[4 *kMaxPtBins]
! hist. of cosXY pointing angle (sig,bkg,tot)
void SetPtWeightsFromFONLL5andTAMUoverLHC16i2abc()
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)
see enum
TH2F * fDalitzPhi[4 *kMaxPtBins]
! dalitz plot via phi (sig,bkg,tot)
const Double_t ptmin
TH2F * fHistSelV0multNTPCout
! histo for V0mult vs #tracks TPCout (sel)
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:257
Bool_t IsEventRejectedDueToPileup() const
Definition: AliRDHFCuts.h:327
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)
short Short_t
Definition: External.C:23
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 fill sparse with Ntracklets
Bool_t GetIsPrimaryWithoutDaughters() const
Definition: AliRDHFCuts.h:277
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)
TString fCentEstName
name of the AliMultSelection object to be considered
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:247
Bool_t fUseWeight
flag to activate cut on V0mult vs #tracks TPCout
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:248
TH1F * fDCAHist[4 *kMaxPtBins]
! hist. for DCA (Prod Cuts)
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:315
Bool_t fDoCutVarHistos
flag to control ntuple writing in MC
const Int_t nbins
Double_t maxMass
TH1F * fPtProng2Hist[4 *kMaxPtBins]
! hist. for DCA (Prod Cuts)
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
void SetPtWeightsFromFONLL5andBAMPSoverLHC16i2abc()
void FillMCGenAccHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader, Double_t nTracklets)
TH2F * fDalitz[4 *kMaxPtBins]
! dalitz plot (sig,bkg,tot)
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:279
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)
Bool_t fFillAcceptanceLevel
flag for usage of sparse for imp. parameter
TString fMultSelectionObjectName
!<!THnSparse for imp. par. on MC
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]