AliPhysics  master (3d17d9d)
AliAnalysisTaskSECharmHadronvnTMVA.cxx
Go to the documentation of this file.
1 /* Copyright(c) 1998-2019, ALICE Experiment at CERN, All rights reserved. *
2  * See cxx source for full Copyright notice */
3 
4 /* $Id$ */
5 
6 //*************************************************************************************
7 // \class AliAnalysisTaskSECharmHadronvnTMVA
8 // \brief task for the analysis of D-meson vn
9 // \authors:
10 // F. Grosa, fabrizio.grosa@cern.ch
11 // F. Catalano, fabio.catalano@cern.ch
12 // A. Festanti, andrea.festanti@cern.ch
13 // G. Luparello, grazia.luparello@cern.ch
14 // F. Prino, prino@to.infn.it
15 // A. Rossi, andrea.rossi@cern.ch
16 // S. Trogolo, stefano.trogolo@cern.ch
18 
19 #include <Riostream.h>
20 #include <TClonesArray.h>
21 #include <TList.h>
22 #include <TFile.h>
23 #include <TDatabasePDG.h>
24 #include <TAxis.h>
25 #include <TSpline.h>
26 #include <TGrid.h>
27 
28 #include <AliLog.h>
29 #include "AliAnalysisManager.h"
30 #include "AliAODHandler.h"
31 #include "AliAODVertex.h"
32 #include "AliAODTrack.h"
33 #include "AliAODRecoDecayHF.h"
35 #include "AliAODRecoCascadeHF.h"
36 #include "AliESDtrack.h"
37 #include "AliVertexerTracks.h"
38 #include "AliAnalysisTaskSE.h"
40 #include "AliRDHFCutsDstoKKpi.h"
42 #include "AliRDHFCutsD0toKpi.h"
43 #include "AliAODEvent.h"
44 #include "AliAODVertex.h"
45 #include "AliAODTrack.h"
46 #include "AliAnalysisVertexingHF.h"
47 
48 #include "AliMultSelection.h"
49 #include "AliVertexingHFUtils.h"
50 
51 #include "AliAnalysisTaskSEHFTenderQnVectors.h"
53 
54 #include "AliAODPidHF.h"
55 
59 //________________________________________________________________________
62  fAOD(nullptr),
63  fOutput(nullptr),
64  fHistNEvents(nullptr),
65  fHistCandVsCent(nullptr),
66  fHistPercqnVsqnVsCentr(nullptr),
67  fHistNtrklVsqnVsCentr(nullptr),
68  fHistMassPtPhiqnCentr(nullptr),
69  fHistMassPtPhiqnCentr_1(nullptr),
70  fHistMassPtPhiqnCentr_2(nullptr),
71  fHistMassPtPhiqnCentr_3(nullptr),
72  fHistMassPtPhiqnCentr_4(nullptr),
73  fHistMassPtPhiqnCentr_5(nullptr),
74  fRDCuts(nullptr),
75  fTenderTaskName("HFTenderQnVectors"),
76  fMinCentr(0.),
77  fMaxCentr(100.),
78  fAODProtection(1),
79  fDaughterTracks(),
80  fHarmonic(2),
81  fCalibType(AliHFQnVectorHandler::kQnFrameworkCalib),
82  fNormMethod(AliHFQnVectorHandler::kQoverM),
83  fOADBFileName(""),
84  fFlowMethod(kEvShapeEP),
85  fEvPlaneDet(kFullV0),
86  fSubEvDetA(kPosTPC),
87  fSubEvDetB(kNegTPC),
88  fqnMeth(kq2TPC),
89  fScalProdLimit(0.4),
90  fPercentileqn(false),
91  fqnSplineFileName(""),
92  fLoadedSplines(false),
93  fDecChannel(kD0toKpi),
94  fLowmasslimit(1.669),
95  fUpmasslimit(2.069),
96  fNMassBins(200),
97  fEtaGapInTPCHalves(0),
98  fRemoveDauFromqn(0),
99  fListRDHFBDT(0),
100  fListBDTNtuple(0),
101  fRemoveSoftPion(false),
102  fEnableDownsamplqn(false),
103  fFracToKeepDownSamplqn(1.1)
104 {
105  // Default constructor
106  for (Int_t ih=0; ih<13; ih++) fBDT1Cut[ih] = -2;
107  for (Int_t in=0; in<13; in++) fBDT2Cut[in] = -2;
108  for(int iHisto=0; iHisto<3; iHisto++) {
109  fHistCentrality[iHisto] = nullptr;
110  fHistEPResolVsCentrVsqn[iHisto] = nullptr;
111  }
112  for(int iDet=0; iDet<6; iDet++) {
113  fHistEvPlaneQncorr[iDet] = nullptr;
114  fHistqnVsCentrPercCalib[iDet] = nullptr;
115  fqnSplinesList[iDet] = nullptr;
116  }
117 }
118 
119 //________________________________________________________________________
121  AliAnalysisTaskSE(name),
122  fAOD(nullptr),
123  fOutput(nullptr),
124  fHistNEvents(nullptr),
125  fHistCandVsCent(nullptr),
126  fHistPercqnVsqnVsCentr(nullptr),
127  fHistNtrklVsqnVsCentr(nullptr),
128  fHistMassPtPhiqnCentr(nullptr),
129  fHistMassPtPhiqnCentr_1(nullptr),
130  fHistMassPtPhiqnCentr_2(nullptr),
131  fHistMassPtPhiqnCentr_3(nullptr),
132  fHistMassPtPhiqnCentr_4(nullptr),
133  fHistMassPtPhiqnCentr_5(nullptr),
134  fRDCuts(rdCuts),
135  fTenderTaskName("HFTenderQnVectors"),
136  fMinCentr(0.),
137  fMaxCentr(100.),
138  fAODProtection(1),
139  fDaughterTracks(),
140  fHarmonic(2),
141  fCalibType(AliHFQnVectorHandler::kQnFrameworkCalib),
142  fNormMethod(AliHFQnVectorHandler::kQoverM),
143  fOADBFileName(""),
144  fFlowMethod(kEvShapeEP),
145  fEvPlaneDet(kFullV0),
146  fSubEvDetA(kPosTPC),
147  fSubEvDetB(kNegTPC),
148  fqnMeth(kq2TPC),
149  fScalProdLimit(0.4),
150  fPercentileqn(false),
151  fqnSplineFileName(""),
152  fLoadedSplines(false),
153  fDecChannel(decaychannel),
154  fLowmasslimit(1.669),
155  fUpmasslimit(2.069),
156  fNMassBins(200),
157  fEtaGapInTPCHalves(0),
158  fRemoveDauFromqn(0),
159  fRemoveSoftPion(false),
160  fEnableDownsamplqn(false),
161  fListRDHFBDT(0),
162  fListBDTNtuple(0),
163  fFracToKeepDownSamplqn(1.1)
164 {
165  // standard constructor
166  for (Int_t ih=0; ih<13; ih++) fBDT1Cut[ih] = -2;
167  for (Int_t in=0; in<13; in++) fBDT2Cut[in] = -2;
168  for(int iHisto=0; iHisto<3; iHisto++) {
169  fHistCentrality[iHisto] = nullptr;
170  fHistEPResolVsCentrVsqn[iHisto] = nullptr;
171  }
172  for(int iDet=0; iDet<6; iDet++) {
173  fHistEvPlaneQncorr[iDet] = nullptr;
174  fHistqnVsCentrPercCalib[iDet] = nullptr;
175  fqnSplinesList[iDet] = nullptr;
176  }
177 
178  int pdg=421;
179  switch(fDecChannel){
180  case kDplustoKpipi:
181  pdg=411;
182  break;
183  case kD0toKpi:
184  pdg=421;
185  break;
186  case kDstartoKpipi:
187  pdg=413;
188  break;
189  case kDstoKKpi:
190  pdg=431;
191  break;
192  }
193  if(pdg==413) SetMassLimits((float)0.1,(float)0.2);
194  else SetMassLimits((float)0.2,pdg); //check range
195 
196  // Output slot #1 writes into a TList container
197  DefineOutput(1,TList::Class()); //Main output
198  // Output slot #3 writes into a AliRDHFCuts container (cuts)
199  switch(fDecChannel){
200  case kDplustoKpipi:
201  DefineOutput(2,AliRDHFCutsDplustoKpipi::Class()); //Cut object for Dplus
202  break;
203  case kD0toKpi:
204  DefineOutput(2,AliRDHFCutsD0toKpi::Class()); //Cut object for D0
205  break;
206  case kDstartoKpipi:
207  DefineOutput(2,AliRDHFCutsDStartoKpipi::Class()); //Cut object for D*
208  break;
209  case kDstoKKpi:
210  DefineOutput(2,AliRDHFCutsDstoKKpi::Class()); //Cut object for Ds
211  break;
212  }
213 
214 }
215 
216 //________________________________________________________________________
218 {
219  // Destructor
220  if(fOutput && !fOutput->IsOwner()){
221  delete fHistNEvents;
222  delete fHistCandVsCent;
223  delete fHistPercqnVsqnVsCentr;
224  delete fHistNtrklVsqnVsCentr;
225  delete fHistMassPtPhiqnCentr;
231 
232  for(int iHisto=0; iHisto<3; iHisto++){
233  delete fHistCentrality[iHisto];
234  delete fHistEPResolVsCentrVsqn[iHisto];
235  }
236  for(int iDet=0; iDet<6; iDet++) {
237  delete fHistEvPlaneQncorr[iDet];
238  delete fHistqnVsCentrPercCalib[iDet];
239  }
240  }
241  delete fOutput;
242  delete fRDCuts;
243 
244  for(int iDet=0; iDet<6; iDet++) {
245  if(fqnSplinesList[iDet] && fLoadedSplines) delete fqnSplinesList[iDet];
246  }
247  if (fListRDHFBDT) {
248  delete fListRDHFBDT;
249  fListRDHFBDT = 0;
250  }
251  if (fListBDTNtuple) {
252  delete fListBDTNtuple;
253  fListBDTNtuple = 0;
254  }
255 }
256 //_________________________________________________________________
258  // Set limits for mass spectra plots
259  float mass=0;
260  int abspdg=TMath::Abs(pdg);
261  mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
262  fUpmasslimit = mass+range;
263  fLowmasslimit = mass-range;
264 }
265 
266 //_________________________________________________________________
267 void AliAnalysisTaskSECharmHadronvnTMVA::SetMassLimits(float lowlimit, float uplimit){
268  // Set limits for mass spectra plots
269  if(uplimit>lowlimit) {
270  fUpmasslimit = uplimit;
271  fLowmasslimit = lowlimit;
272  }
273 }
274 
275 //________________________________________________________________________
277 {
278  // Initialization
281 
282  switch(fDecChannel) {
283  case kDplustoKpipi:
284  {
285  AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCuts)));
286  PostData(2,copycut);
287  }
288  break;
289  case kD0toKpi:
290  {
291  AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCuts)));
292  PostData(2,copycut);
293  }
294  break;
295  case kDstartoKpipi:
296  {
297  AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCuts)));
298  PostData(2,copycut);
299  }
300  break;
301  case kDstoKKpi:
302  {
303  AliRDHFCutsDstoKKpi* copycut=new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fRDCuts)));
304  PostData(2,copycut);
305  }
306  break;
307  }
308 
309  return;
310 }
311 
312 //________________________________________________________________________
314 {
315  // Create the output container
316  fOutput = new TList();
317  fOutput->SetOwner();
318  fOutput->SetName("MainOutput");
319  // fListRDHFBDT->SetOwner();
320  fHistNEvents = new TH1F("fHistNEvents", "Number of AODs scanned",16,-0.5,15.5);
321  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsRead");
322  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents Matched dAOD");
323  fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents Mismatched dAOD");
324  fHistNEvents->GetXaxis()->SetBinLabel(4,"nEventsAnal");
325  fHistNEvents->GetXaxis()->SetBinLabel(5,"n. passing IsEvSelected");
326  fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected due to trigger");
327  fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected due to not reco vertex");
328  fHistNEvents->GetXaxis()->SetBinLabel(8,"n. rejected for contr vertex");
329  fHistNEvents->GetXaxis()->SetBinLabel(9,"n. rejected for vertex out of accept");
330  fHistNEvents->GetXaxis()->SetBinLabel(10,"n. rejected for pileup events");
331  fHistNEvents->GetXaxis()->SetBinLabel(11,Form("no. of out %.0f-%.0f%% centrality events",fMinCentr,fMaxCentr));
332  fHistNEvents->GetXaxis()->SetBinLabel(12,"n. rejected for bad cent corr");
333  fHistNEvents->GetXaxis()->SetBinLabel(13,"non valid TPC EP");
334  fHistNEvents->GetXaxis()->SetBinLabel(14,"bad event plane");
335  fHistNEvents->GetXaxis()->SetBinLabel(15,"no. of sel. candidates");
336  fHistNEvents->GetXaxis()->SetBinLabel(16,"no. cand. out of pt bounds");
337  fHistNEvents->GetXaxis()->SetNdivisions(1,false);
338  fOutput->Add(fHistNEvents);
339 
340  fHistCentrality[0]=new TH1F("hCentr","centrality",10000,0.,100.);
341  fHistCentrality[1]=new TH1F("hCentr(selectedCent)","centrality(selectedCent)",10000,0.,100.);
342  fHistCentrality[2]=new TH1F("hCentr(OutofCent)","centrality(OutofCent)",10000,0.,100.);
343  for(int iHisto=0; iHisto<3; iHisto++){
344  fOutput->Add(fHistCentrality[iHisto]);
345  }
346 
347  const int ncentbins = static_cast<int>(fMaxCentr-fMinCentr);
348 
349  fHistCandVsCent=new TH2F("hCandVsCent","number of selected candidates vs. centrality;centrality(%);number of candidates",ncentbins,fMinCentr,fMaxCentr,101,-0.5,100.5);
350  fOutput->Add(fHistCandVsCent);
351 
352  TString detConfName[6] = {"TPC","TPCPosEta","TPCNegEta","V0","V0A","V0C"};
353 
354  TString qnaxisname=Form("#it{q}_{%d}",fHarmonic);
355  switch(fqnMeth) {
356  case kq2TPC:
357  qnaxisname=Form("#it{q}_{%d}^{TPC}",fHarmonic);
358  break;
359  case kq2PosTPC:
360  qnaxisname=Form("#it{q}_{%d}^{TPCPosEta}",fHarmonic);
361  break;
362  case kq2NegTPC:
363  qnaxisname=Form("#it{q}_{%d}^{TPCNegEta}",fHarmonic);
364  break;
365  case kq2VZERO:
366  qnaxisname=Form("#it{q}_{%d}^{V0}",fHarmonic);
367  break;
368  case kq2VZEROA:
369  qnaxisname=Form("#it{q}_{%d}^{V0A}",fHarmonic);
370  break;
371  case kq2VZEROC:
372  qnaxisname=Form("#it{q}_{%d}^{V0C}",fHarmonic);
373  break;
374  }
375 
376  TString qnpercaxisname = qnaxisname + " (%)";
377  TString qnaxisnamefill = qnaxisname;
378  if(fPercentileqn)
379  qnaxisnamefill = qnpercaxisname;
380 
381  int nqnbins=1; //single bin if unbiased analysis
382  double qnmin = 0.;
383  double qnmax = 15.;
385  if(!fPercentileqn) {
386  nqnbins=300;
387  }
388  else {
389  nqnbins=100;
390  qnmin = 0.;
391  qnmax = 100.;
392  }
393  }
394 
395  for(int iDet = 0; iDet < 6; iDet++) {
396  fHistEvPlaneQncorr[iDet] = new TH3F(Form("fHistEvPlaneQncorr%sVsqnVsCent",detConfName[iDet].Data()),Form("hEvPlaneQncorr%s;centrality(%%);%s;#psi_{%d}",detConfName[iDet].Data(),qnaxisnamefill.Data(),fHarmonic),ncentbins,fMinCentr,fMaxCentr,nqnbins,qnmin,qnmax,200,0.,TMath::Pi());
397  fOutput->Add(fHistEvPlaneQncorr[iDet]);
398 
399  // histos for qn vs. centrality with fine binning (for qn percentiles calibration)
401  fHistqnVsCentrPercCalib[iDet] = new TH2F(Form("fHistqnVsCentr%s",detConfName[iDet].Data()),Form("#it{q}_{%d}^{%s} vs. centrality;centrality(%%);#it{q}_{%d}^{%s}",fHarmonic,detConfName[iDet].Data(),fHarmonic,detConfName[iDet].Data()),ncentbins,fMinCentr,fMaxCentr,15000,0.,15.);
402  fOutput->Add(fHistqnVsCentrPercCalib[iDet]);
403  }
404  }
405 
406  //qn percentile vs. qn vs. centrality
408  fHistNtrklVsqnVsCentr = new TH3F("fHistNtrklVsqnVsCentr",Form("#it{N}_{tracklets} vs. %s vs. centrality;centrality (%%);%s;#it{N}_{tracklets}",qnpercaxisname.Data(),qnpercaxisname.Data()),ncentbins,fMinCentr,fMaxCentr,nqnbins,qnmin,qnmax,500,-0.5,4999.5);
410  if(fPercentileqn) {
411  fHistPercqnVsqnVsCentr = new TH3F("fHistPercqnVsqnVsCentr",Form("%s vs. %s vs. centrality;centrality (%%);%s;%s",qnpercaxisname.Data(),qnaxisname.Data(),qnaxisname.Data(),qnpercaxisname.Data()),ncentbins,fMinCentr,fMaxCentr,300,0,15,nqnbins,qnmin,qnmax);
413  }
414  }
415 
416  //EP / Qn resolutions
417  TString detLabels[3][2] = {{"A","B"},{"A","C"},{"B","C"}};
418  for(int iResoHisto=0; iResoHisto<3; iResoHisto++) {
420  fHistEPResolVsCentrVsqn[iResoHisto] = new TH3F(Form("fHistEvPlaneReso%d",iResoHisto+1),Form("Event plane angle Resolution;centrality (%%);%s;cos2(#psi_{%s}-#psi_{%s})",qnaxisnamefill.Data(),detLabels[iResoHisto][0].Data(),detLabels[iResoHisto][1].Data()),ncentbins,fMinCentr,fMaxCentr,nqnbins,qnmin,qnmax,220,-1.1,1.1);
421  }
422  else if(fFlowMethod==kSP || fFlowMethod==kEvShapeSP) {
423  fHistEPResolVsCentrVsqn[iResoHisto] = new TH3F(Form("hScalProdQnVectors%d",iResoHisto+1),Form("Scalar product between Q-vectors;centrality (%%);%s;(Q{%s}Q{%s})",qnaxisnamefill.Data(),detLabels[iResoHisto][0].Data(),detLabels[iResoHisto][1].Data()),ncentbins,fMinCentr,fMaxCentr,nqnbins,qnmin,qnmax,200,-fScalProdLimit*fScalProdLimit,fScalProdLimit*fScalProdLimit);
424  }
425  fOutput->Add(fHistEPResolVsCentrVsqn[iResoHisto]);
426  }
427 
428  //Sparse with D-meson candidates and useful azimuthal info
429  int nptbins = 200;
430  double ptmin = 0.;
431  double ptmax = 50.;
432 
433  int ndeltaphibins = 1;
434  double mindeltaphi = -1.;
435  double maxdeltaphi = -1.;
436 
437  int nfphibins = 100;
438  double fphimin = -1.;
439  double fphimax = 1.;
440 
441  TString deltaphiname = "";
442  TString nfphiname1 = Form("Cos(%d#varphi_{D})",fHarmonic);
443  TString nfphiname2 = Form("Sin(%d#varphi_{D})",fHarmonic);
445  ndeltaphibins = 96;
446  mindeltaphi = 0.;
447  maxdeltaphi = TMath::Pi();
448  deltaphiname = "#Delta#varphi";
449  }
450  else if(fFlowMethod==kSP || fFlowMethod==kEvShapeSP) {
451  ndeltaphibins = 100;
452  mindeltaphi = -fScalProdLimit*fScalProdLimit;
453  maxdeltaphi = fScalProdLimit*fScalProdLimit;
454  deltaphiname = Form("u_{%d,D}Q_{%d,A}",fHarmonic,fHarmonic);
455  fphimin = -fScalProdLimit*fScalProdLimit;
456  fphimax = fScalProdLimit*fScalProdLimit;
457  nfphiname1 = Form("Cos(%d#varphi_{D})Q_{%d,y,A}",fHarmonic,fHarmonic);
458  nfphiname2 = Form("Sin(%d#varphi_{D})Q_{%d,y,A}",fHarmonic,fHarmonic);
459  }
461  ndeltaphibins = 100;
462  mindeltaphi = -1.;
463  maxdeltaphi = 1.;
464  deltaphiname = Form("cos(%d#Delta#varphi)",fHarmonic);
465  }
466 
467  int nphibins = 18;
468  double phimin = 0.;
469  double phimax = 2*TMath::Pi();
470 
471  int nNtrkBins = 100;
472  double Ntrkmin = 0.;
473  double Ntrkmax = 5000.;
474 
475  TString massaxisname = "";
476  if(fDecChannel==0) massaxisname = "#it{M}(K#pi#pi) (GeV/#it{c}^{2})";
477  else if(fDecChannel==1) massaxisname = "#it{M}(K#pi) (GeV/#it{c}^{2})";
478  else if(fDecChannel==2) massaxisname = "#it{M}(K#pi#pi)-#it{M}(K#pi) (GeV/#it{c}^{2})";
479  else if(fDecChannel==3) massaxisname = "#it{M}(KK#pi) (GeV/#it{c}^{2})";
480 
481  int naxes=kVarForSparse;
482 
483  int nbins[kVarForSparse] = {fNMassBins, nptbins, ndeltaphibins, nfphibins, nfphibins, nphibins, ncentbins, nNtrkBins, nqnbins};
484  double xmin[kVarForSparse] = {fLowmasslimit, ptmin, mindeltaphi, fphimin, fphimin, phimin, fMinCentr, Ntrkmin, qnmin};
485  double xmax[kVarForSparse] = {fUpmasslimit, ptmax, maxdeltaphi, fphimax, fphimax, phimax, fMaxCentr, Ntrkmax, qnmax};
486  TString axTit[kVarForSparse] = {massaxisname, "#it{p}_{T} (GeV/#it{c})", deltaphiname, nfphiname1, nfphiname2, "#varphi_{D}", "Centrality (%)", "#it{N}_{tracklets}", qnaxisnamefill};
487 
488  fHistMassPtPhiqnCentr = new THnSparseF("fHistMassPtPhiqnCentr",Form("InvMass vs. #it{p}_{T} vs. %s vs. centr vs. #it{q}_{%d} ",deltaphiname.Data(),fHarmonic),naxes,nbins,xmin,xmax);
489  fHistMassPtPhiqnCentr_1 = new THnSparseF("fHistMassPtPhiqnCentr_1",Form("InvMass vs. #it{p}_{T} vs. %s vs. centr vs. #it{q}_{%d} ",deltaphiname.Data(),fHarmonic),naxes,nbins,xmin,xmax);
490  fHistMassPtPhiqnCentr_2 = new THnSparseF("fHistMassPtPhiqnCentr_2",Form("InvMass vs. #it{p}_{T} vs. %s vs. centr vs. #it{q}_{%d} ",deltaphiname.Data(),fHarmonic),naxes,nbins,xmin,xmax);
491  fHistMassPtPhiqnCentr_3 = new THnSparseF("fHistMassPtPhiqnCentr_3",Form("InvMass vs. #it{p}_{T} vs. %s vs. centr vs. #it{q}_{%d} ",deltaphiname.Data(),fHarmonic),naxes,nbins,xmin,xmax);
492  fHistMassPtPhiqnCentr_4 = new THnSparseF("fHistMassPtPhiqnCentr_4",Form("InvMass vs. #it{p}_{T} vs. %s vs. centr vs. #it{q}_{%d} ",deltaphiname.Data(),fHarmonic),naxes,nbins,xmin,xmax);
493  fHistMassPtPhiqnCentr_5 = new THnSparseF("fHistMassPtPhiqnCentr_5",Form("InvMass vs. #it{p}_{T} vs. %s vs. centr vs. #it{q}_{%d} ",deltaphiname.Data(),fHarmonic),naxes,nbins,xmin,xmax);
494 
495  for(int iAx=0; iAx<naxes; iAx++){
496  fHistMassPtPhiqnCentr->GetAxis(iAx)->SetTitle(axTit[iAx].Data());
497  fHistMassPtPhiqnCentr_1->GetAxis(iAx)->SetTitle(axTit[iAx].Data());
498  fHistMassPtPhiqnCentr_2->GetAxis(iAx)->SetTitle(axTit[iAx].Data());
499  fHistMassPtPhiqnCentr_3->GetAxis(iAx)->SetTitle(axTit[iAx].Data());
500  fHistMassPtPhiqnCentr_4->GetAxis(iAx)->SetTitle(axTit[iAx].Data());
501  fHistMassPtPhiqnCentr_5->GetAxis(iAx)->SetTitle(axTit[iAx].Data());
502  }
503 
510 
511  // fListBDTNtuple = new TList();fListBDTNtuple->SetOwner(); fListBDTNtuple->SetName("NtupleList");
512  // fListRDHFBDT->SetOwner(); fListBDTNtuple->SetName("BDTList");
513  //ML model
514 
515  PostData(1,fOutput);
516  // PostData(3,fListBDTNtuple);
517  return;
518 }
519 
520 //________________________________________________________________________
522 {
523 
524  // Execute analysis for current event:
525  // heavy flavor candidates association to MC truth
526  fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
527  fHistNEvents->Fill(0);
528 
529  // Protection against the mismatch of candidate TRefs:
530  // Check if AOD and corresponding deltaAOD files contain the same number of events.
531  // In case of discrepancy the event is rejected.
532  if(fAODProtection>=0) {
533  // Protection against different number of events in the AOD and deltaAOD
534  // In case of discrepancy the event is rejected.
535  int matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
536  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
537  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
538  fHistNEvents->Fill(2);
539  return;
540  }
541  fHistNEvents->Fill(1);
542  }
543 
544  // Post the data already here
545  PostData(1,fOutput);
546 
547  TClonesArray *arrayProng = nullptr, *arrayD0toKpi = nullptr;
548  int absPdgMom=0;
549  int nDau=0;
550  if(!fAOD && AODEvent() && IsStandardAOD()) {
551  // In case there is an AOD handler writing a standard AOD, use the AOD
552  // event in memory rather than the input (ESD) event.
553  fAOD = dynamic_cast<AliAODEvent*> (AODEvent());
554  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
555  // have to taken from the AOD event hold by the AliAODExtension
556  AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
557  if(aodHandler->GetExtensions()) {
558 
559  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
560  AliAODEvent *aodFromExt = ext->GetAOD();
561 
562  switch(fDecChannel) {
563  case kDplustoKpipi:
564  absPdgMom=411;
565  nDau=3;
566  arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
567  break;
568  case kD0toKpi:
569  absPdgMom=421;
570  nDau=2;
571  arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
572  break;
573  case kDstartoKpipi:
574  absPdgMom=413;
575  nDau=3;
576  arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
577  arrayD0toKpi=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
578  break;
579  case kDstoKKpi:
580  absPdgMom=431;
581  nDau=3;
582  arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
583  break;
584  }
585  }
586  }
587  else if(fAOD){
588  switch(fDecChannel) {
589  case kDplustoKpipi:
590  absPdgMom=411;
591  nDau = 3;
592  arrayProng=(TClonesArray*)fAOD->GetList()->FindObject("Charm3Prong");
593  break;
594  case kD0toKpi:
595  absPdgMom=421;
596  nDau = 2;
597  arrayProng=(TClonesArray*)fAOD->GetList()->FindObject("D0toKpi");
598  break;
599  case kDstartoKpipi:
600  absPdgMom=413;
601  nDau = 3;
602  arrayProng=(TClonesArray*)fAOD->GetList()->FindObject("Dstar");
603  arrayD0toKpi=(TClonesArray*)fAOD->GetList()->FindObject("D0toKpi");
604  break;
605  case kDstoKKpi:
606  absPdgMom=431;
607  nDau = 3;
608  arrayProng=(TClonesArray*)fAOD->GetList()->FindObject("Charm3Prong");
609  break;
610  }
611  }
612 
613  if(!fAOD || !arrayProng || (!arrayD0toKpi && fDecChannel==kDstartoKpipi)) {
614  AliError("AliAnalysisTaskSECharmHadronvnTMVA::UserExec:Branch not found!\n");
615  return;
616  }
617 
618  // fix for temporary bug in ESDfilter
619  // the AODs with null vertex pointer didn't pass the PhysSel
620  if(!fAOD->GetPrimaryVertex() || TMath::Abs(fAOD->GetMagneticField())<0.001) return;
621 
622  fHistNEvents->Fill(3);
623 
624  double evCentr = fRDCuts->GetCentrality(fAOD);
625  bool isEvSel = fRDCuts->IsEventSelected(fAOD);
626  fHistCentrality[0]->Fill(evCentr);
627 
628  if(!isEvSel) {
635  fHistNEvents->Fill(10);
636  fHistCentrality[2]->Fill(evCentr);
637  }
638  return;
639  }
640 
641  fHistNEvents->Fill(4);
642  fHistCentrality[1]->Fill(evCentr);
643 
645 
646  //Get Qn-vectors from tender task
647  AliHFQnVectorHandler *HFQnVectorHandler = nullptr;
648  bool isHandlerFound = false;
649  AliAnalysisTaskSEHFTenderQnVectors *HFQnVectorTask = dynamic_cast<AliAnalysisTaskSEHFTenderQnVectors*>(AliAnalysisManager::GetAnalysisManager()->GetTask(fTenderTaskName.Data()));
650  if(HFQnVectorTask) {
651  HFQnVectorHandler = HFQnVectorTask->GetQnVectorHandler();
652 
653  if(fPercentileqn && !fqnSplinesList[0]) {
654  for(int iDet=0; iDet<6; iDet++) {
655  fqnSplinesList[iDet] = dynamic_cast<TList*>(HFQnVectorTask->GetSplineForqnPercentileList(iDet));
656  }
657  }
658  }
659 
660  if(HFQnVectorHandler) {
661  isHandlerFound = true;
662  if(HFQnVectorHandler->GetHarmonic() != fHarmonic) {
663  AliWarning("Harmonic of task and Qn-vector handler not consistent!");
664  return;
665  }
666  if(HFQnVectorHandler->GetCalibrationType() != fCalibType) {
667  AliWarning("Calibration strategy of task and Qn-vector handler not consistent!");
668  return;
669  }
670  if(HFQnVectorHandler->GetNormalisationMethod() != fNormMethod) {
671  AliWarning("Normalisation method of task and Qn-vector handler not consistent!");
672  return;
673  }
674  if(fCalibType==AliHFQnVectorHandler::kQnCalib && HFQnVectorHandler->GetCalibrationsOADBFileName() != fOADBFileName) {
675  AliWarning("OADB file name for calibrations of task and Qn-vector handler not consistent!");
676  return;
677  }
678  }
679  else { //create a new handler if not found in tender task
680  AliWarning("Qn-vector tender task not found! Create a new one");
681  HFQnVectorHandler = new AliHFQnVectorHandler(fCalibType,fNormMethod,fHarmonic,fOADBFileName);
682  HFQnVectorHandler->SetAODEvent(fAOD);
683  HFQnVectorHandler->ComputeCalibratedQnVectorTPC();
684  HFQnVectorHandler->ComputeCalibratedQnVectorV0();
685  }
686 
687  if(fPercentileqn && !fqnSplinesList[0]) { //probably qn splines not found in tender task, let's load them here
689  }
690 
691  double QnFullTPC[2], QnPosTPC[2], QnNegTPC[2];
692  double QnFullV0[2], QnV0A[2], QnV0C[2];
693  double MultQnFullTPC = -1., MultQnPosTPC = -1., MultQnNegTPC = -1.;
694  double MultQnFullV0 = -1., MultQnV0A = -1., MultQnV0C = -1.;
695  double PsinFullTPC = -1., PsinPosTPC = -1., PsinNegTPC = -1.;
696  double PsinFullV0 = -1., PsinV0A = -1., PsinV0C = -1.;
697 
698  //get the unnormalised Qn-vectors --> normalisation can be done in the task
699  HFQnVectorHandler->GetUnNormQnVecTPC(QnFullTPC,QnPosTPC,QnNegTPC);
700  HFQnVectorHandler->GetUnNormQnVecV0(QnFullV0,QnV0A,QnV0C);
701 
702  HFQnVectorHandler->GetMultQnVecTPC(MultQnFullTPC,MultQnPosTPC,MultQnNegTPC);
703  HFQnVectorHandler->GetMultQnVecV0(MultQnFullV0,MultQnV0A,MultQnV0C);
704 
705  HFQnVectorHandler->GetEventPlaneAngleTPC(PsinFullTPC,PsinPosTPC,PsinNegTPC);
706  HFQnVectorHandler->GetEventPlaneAngleV0(PsinFullV0,PsinV0A,PsinV0C);
707 
708  double mainQn[2], SubAQn[2], SubBQn[2], SubCQn[2];
709  double mainPsin = -1., SubAPsin = -1., SubBPsin = -1., SubCPsin = -1.;
710  double mainMultQn = -1., SubAMultQn = -1., SubBMultQn = -1., SubCMultQn = -1.;
711  GetMainQnVectorInfo(mainPsin, mainMultQn, mainQn, SubAPsin, SubAMultQn, SubAQn, SubBPsin, SubBMultQn, SubBQn, HFQnVectorHandler);
712  SubCPsin = mainPsin;
713  SubCMultQn = mainMultQn;
714  SubCQn[0] = mainQn[0];
715  SubCQn[1] = mainQn[1];
716 
717  int nsubevents = 3;
719  nsubevents = 2;
720 
721  double qnFullTPC = -1., qnPosTPC = -1., qnNegTPC = -1.;
722  double qnFullV0 = -1., qnV0A = -1., qnV0C = -1.;
723  double mainqn = -1., mainpercqn = -1.;
724  TSpline3* qnspline = nullptr;
725 
726  //define vars for candidate loop already here, in case daughter tracks removed from qn
728  AliAODRecoDecayHF *d = nullptr;
729  AliAODRecoDecayHF2Prong *dD0 = nullptr;
730 
731  int nCand = arrayProng->GetEntriesFast();
732  bool alreadyLooped = false;
733  vector<int> isSelByPrevLoop;
734  vector<int> isBDTByPrevLoop;
735  vector<AliAODTrack*> trackstoremove;
736 
738  HFQnVectorHandler->GetqnTPC(qnFullTPC,qnPosTPC,qnNegTPC);
739  HFQnVectorHandler->GetqnV0(qnFullV0,qnV0A,qnV0C);
741 
742  //random downsampling of tracks for qn
744  for(int iTrack=0; iTrack<fAOD->GetNumberOfTracks(); iTrack++) {
745  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTrack));
746  double pseudoRand = track->Pt()*1000.-(long)(track->Pt()*1000);
747  if(pseudoRand>fFracToKeepDownSamplqn)
748  trackstoremove.push_back(track);
749  }
750  }
751 
752  //remove all daughter tracks of the analysed event, if enabled
753  if(fRemoveDauFromqn==2) {
754  alreadyLooped = true;
755  for (int iCand = 0; iCand < nCand; iCand++) { //loop already here on candidates
756  d = dynamic_cast<AliAODRecoDecayHF*>(arrayProng->UncheckedAt(iCand));
757  AliAODRecoDecayHF2Prong *dd0 = (AliAODRecoDecayHF2Prong*)arrayProng->UncheckedAt(iCand);
759  if(d->GetIsFilled()<1)
760  dD0 = dynamic_cast<AliAODRecoDecayHF2Prong*>(arrayD0toKpi->At(d->GetProngID(1)));
761  else
762  dD0 = (dynamic_cast<AliAODRecoCascadeHF*>(d))->Get2Prong();
763  }
764  double modelPred[2] = {-1., -1.};
765  int isSel = IsCandidateSelected(d, nDau, absPdgMom, vHF, dD0);
766  int isBDT = ProcessBDT(fAOD,dd0,isSel,vHF);
767  isSelByPrevLoop.push_back(isSel);
768  isBDTByPrevLoop.push_back(isBDT);
769  if(!isSel) continue;
770  if(!isBDT) continue;
771 
772  GetDaughterTracksToRemove(d,nDau,trackstoremove);
773  }
774  }
775 
776  double QnFullTPCDownSampl[2], QnPosTPCDownSampl[2], QnNegTPCDownSampl[2];
777  double MultQnFullTPCDownSampl = -1., MultQnPosTPCDownSampl = -1., MultQnNegTPCDownSampl = -1.;
778  HFQnVectorHandler->RemoveTracksFromQnTPC(trackstoremove, QnFullTPCDownSampl, QnPosTPCDownSampl, QnNegTPCDownSampl, MultQnFullTPCDownSampl, MultQnPosTPCDownSampl, MultQnNegTPCDownSampl, true);
779  qnFullTPC = TMath::Sqrt((QnFullTPCDownSampl[0]*QnFullTPCDownSampl[0]+QnFullTPCDownSampl[1]*QnFullTPCDownSampl[1])/MultQnFullTPCDownSampl);
780  qnPosTPC = TMath::Sqrt((QnPosTPCDownSampl[0]*QnPosTPCDownSampl[0]+QnPosTPCDownSampl[1]*QnPosTPCDownSampl[1])/MultQnPosTPCDownSampl);
781  qnNegTPC = TMath::Sqrt((QnNegTPCDownSampl[0]*QnNegTPCDownSampl[0]+QnNegTPCDownSampl[1]*QnNegTPCDownSampl[1])/MultQnNegTPCDownSampl);
782  }
783 
784  TAxis* centraxis = fHistqnVsCentrPercCalib[0]->GetXaxis();
785  int centrbin = centraxis->FindBin(evCentr);
786  double centrbinmin = centraxis->GetBinLowEdge(centrbin);
787  double centrbinmax = centrbinmin+centraxis->GetBinWidth(centrbin);
788 
789  switch(fqnMeth) {
790  case kq2TPC:
791  mainqn = qnFullTPC;
792  if(fPercentileqn) qnspline = static_cast<TSpline3*>(fqnSplinesList[0]->FindObject(Form("sq2Int_centr_%0.f_%0.f",centrbinmin,centrbinmax)));
793  break;
794  case kq2PosTPC:
795  mainqn = qnPosTPC;
796  if(fPercentileqn) qnspline = static_cast<TSpline3*>(fqnSplinesList[1]->FindObject(Form("sq2Int_centr_%0.f_%0.f",centrbinmin,centrbinmax)));
797  break;
798  case kq2NegTPC:
799  mainqn = qnNegTPC;
800  if(fPercentileqn) qnspline = static_cast<TSpline3*>(fqnSplinesList[2]->FindObject(Form("sq2Int_centr_%0.f_%0.f",centrbinmin,centrbinmax)));
801  break;
802  case kq2VZERO:
803  mainqn = qnFullV0;
804  if(fPercentileqn) qnspline = static_cast<TSpline3*>(fqnSplinesList[3]->FindObject(Form("sq2Int_centr_%0.f_%0.f",centrbinmin,centrbinmax)));
805  break;
806  case kq2VZEROA:
807  mainqn = qnV0A;
808  if(fPercentileqn) qnspline = static_cast<TSpline3*>(fqnSplinesList[4]->FindObject(Form("sq2Int_centr_%0.f_%0.f",centrbinmin,centrbinmax)));
809  break;
810  case kq2VZEROC:
811  mainqn = qnV0C;
812  if(fPercentileqn) qnspline = static_cast<TSpline3*>(fqnSplinesList[5]->FindObject(Form("sq2Int_centr_%0.f_%0.f",centrbinmin,centrbinmax)));
813  break;
814  }
815  if(fPercentileqn) {
816  if(qnspline)
817  mainpercqn = qnspline->Eval(mainqn);
818  else
819  AliWarning("Centrality binning and centrality intervals of qn splines do not match!");
820  }
821  else {
822  mainpercqn = mainqn;
823  }
824  }
825  //Fill event-based histos
826  //EP and q2 histos
827  fHistEvPlaneQncorr[0]->Fill(evCentr,mainpercqn,PsinFullTPC);
828  fHistEvPlaneQncorr[1]->Fill(evCentr,mainpercqn,PsinPosTPC);
829  fHistEvPlaneQncorr[2]->Fill(evCentr,mainpercqn,PsinNegTPC);
830  fHistEvPlaneQncorr[3]->Fill(evCentr,mainpercqn,PsinFullV0);
831  fHistEvPlaneQncorr[4]->Fill(evCentr,mainpercqn,PsinV0A);
832  fHistEvPlaneQncorr[5]->Fill(evCentr,mainpercqn,PsinV0C);
833 
835  fHistqnVsCentrPercCalib[0]->Fill(evCentr,qnFullTPC);
836  fHistqnVsCentrPercCalib[1]->Fill(evCentr,qnPosTPC);
837  fHistqnVsCentrPercCalib[2]->Fill(evCentr,qnNegTPC);
838  fHistqnVsCentrPercCalib[3]->Fill(evCentr,qnFullV0);
839  fHistqnVsCentrPercCalib[4]->Fill(evCentr,qnV0A);
840  fHistqnVsCentrPercCalib[5]->Fill(evCentr,qnV0C);
841 
842  fHistNtrklVsqnVsCentr->Fill(evCentr,mainpercqn,tracklets);
843 
844  if(fPercentileqn) {
845  fHistPercqnVsqnVsCentr->Fill(evCentr,mainqn,mainpercqn);
846  }
847  }
848 
849  //EP / Qn resolution histos
851  fHistEPResolVsCentrVsqn[0]->Fill(evCentr,mainpercqn,TMath::Cos(fHarmonic*GetDeltaPsiSubInRange(SubAPsin,SubBPsin)));
852  if(nsubevents==3) {
853  fHistEPResolVsCentrVsqn[1]->Fill(evCentr,mainpercqn,TMath::Cos(fHarmonic*GetDeltaPsiSubInRange(SubAPsin,SubCPsin)));
854  fHistEPResolVsCentrVsqn[2]->Fill(evCentr,mainpercqn,TMath::Cos(fHarmonic*GetDeltaPsiSubInRange(SubBPsin,SubCPsin)));
855  }
856  }
857  else if(fFlowMethod==kSP || fFlowMethod==kEvShapeSP) {
858  fHistEPResolVsCentrVsqn[0]->Fill(evCentr,mainpercqn,(SubAQn[0]*SubBQn[0]+SubAQn[1]*SubBQn[1])/(SubAMultQn * SubBMultQn));
859  if(nsubevents==3) {
860  fHistEPResolVsCentrVsqn[1]->Fill(evCentr,mainpercqn,(SubAQn[0]*SubCQn[0]+SubAQn[1]*SubCQn[1])/(SubAMultQn * SubCMultQn));
861  fHistEPResolVsCentrVsqn[2]->Fill(evCentr,mainpercqn,(SubBQn[0]*SubCQn[0]+SubBQn[1]*SubCQn[1])/(SubCMultQn * SubBMultQn));
862  }
863  }
864  //Loop on D candidates
865  int nSelCand = 0;
866  for (int iCand = 0; iCand < nCand; iCand++) {
867  d = dynamic_cast<AliAODRecoDecayHF*>(arrayProng->UncheckedAt(iCand));
868  AliAODRecoDecayHF2Prong *dd0 = (AliAODRecoDecayHF2Prong*)arrayProng->UncheckedAt(iCand);
870  if(d->GetIsFilled()<1)
871  dD0 = dynamic_cast<AliAODRecoDecayHF2Prong*>(arrayD0toKpi->At(d->GetProngID(1)));
872  else
873  dD0 = (dynamic_cast<AliAODRecoCascadeHF*>(d))->Get2Prong();
874  }
875 
876  int isSelected = 0;
877  int isBDT = 0;
878  double modelPred[2] = {-1., -1.};
879  if(alreadyLooped) {//check already here to avoid application of cuts twice
880  isSelected = isSelByPrevLoop[iCand];
881  isBDT = isBDTByPrevLoop[iCand];
882  }
883  else {
884  isSelected = IsCandidateSelected(d, nDau, absPdgMom, vHF, dD0);
885  isBDT = ProcessBDT(fAOD,dd0,isSelected,vHF);
886  // printf("isBDT = %d\n",isBDT);
887  // printf("isSelected = %d\n",isSelected);
888 
889  }
890  if(!isSelected) continue;
891  if(!isBDT) continue;
892 
893  nSelCand++;
894  fHistNEvents->Fill(14); // candidate selected (including ML selection if used)
895 
896  float* invMass = nullptr;
897  int nmasses = 0;
898  CalculateInvMasses(d,invMass,nmasses);
899 
900  double ptD = d->Pt();
901  double phiD = d->Phi();
902 
903  double deltaphi = GetPhiInRange(phiD-mainPsin);
904  double scalprod = (TMath::Cos(fHarmonic*phiD)*mainQn[0]+TMath::Sin(fHarmonic*phiD)*mainQn[1]) / mainMultQn;
905  double cosndeltaphi = TMath::Cos(fHarmonic*deltaphi);
906  double vnfunc = -999.;
907  double phifunc1 = -999.;
908  double phifunc2 = -999.;
910  vnfunc = deltaphi;
911  phifunc1 = TMath::Cos(fHarmonic*phiD);
912  phifunc2 = TMath::Sin(fHarmonic*phiD);
913  }
914  else if(fFlowMethod==kSP || fFlowMethod==kEvShapeSP) {
915  vnfunc = scalprod;
916  phifunc1 = TMath::Cos(fHarmonic*phiD)*mainQn[1] / mainMultQn;
917  phifunc2 = TMath::Sin(fHarmonic*phiD)*mainQn[0] / mainMultQn;
918  }
920  vnfunc = cosndeltaphi;
921  phifunc1 = TMath::Cos(fHarmonic*phiD);
922  phifunc2 = TMath::Sin(fHarmonic*phiD);
923  }
924 
925  //remove daughter tracks from qn (on top of random downsampling, if enabled)
926  double candQnFullTPC[2], candQnPosTPC[2], candQnNegTPC[2];
927  double candMultQnFullTPC = -1., candMultQnPosTPC = -1., candMultQnNegTPC = -1.;
928  double candqn = -1., candpercqn = -1.;
930 
931  vector<AliAODTrack*> trackstoremoveDau(trackstoremove);
932  GetDaughterTracksToRemove(d,nDau,trackstoremoveDau);
933  HFQnVectorHandler->RemoveTracksFromQnTPC(trackstoremoveDau, candQnFullTPC, candQnPosTPC, candQnNegTPC, candMultQnFullTPC, candMultQnPosTPC, candMultQnNegTPC, true);
934  switch(fqnMeth) {
935  case kq2TPC:
936  candqn = TMath::Sqrt((candQnFullTPC[0]*candQnFullTPC[0]+candQnFullTPC[1]*candQnFullTPC[1])/candMultQnFullTPC);
937  break;
938  case kq2PosTPC:
939  candqn = TMath::Sqrt((candQnPosTPC[0]*candQnPosTPC[0]+candQnPosTPC[1]*candQnPosTPC[1])/candMultQnPosTPC);
940  break;
941  case kq2NegTPC:
942  candqn = TMath::Sqrt((candQnNegTPC[0]*candQnNegTPC[0]+candQnNegTPC[1]*candQnNegTPC[1])/candMultQnNegTPC);
943  break;
944  }
945  trackstoremoveDau.clear();
946  if(fPercentileqn) {
947  if(qnspline)
948  candpercqn = qnspline->Eval(candqn);
949  else
950  AliWarning("Centrality binning and centrality intervals of qn splines do not match!");
951  }
952  else {
953  candpercqn = candqn;
954  }
955  }
956  else {
957  candpercqn = mainpercqn;
958  }
959 
960  switch(fDecChannel) {
961  case kDplustoKpipi:
962  {
963  double sparsearray[9] = {invMass[0],ptD,vnfunc,phifunc1,phifunc2,phiD,evCentr,static_cast<double>(tracklets),candpercqn};
964  fHistMassPtPhiqnCentr->Fill(sparsearray);
965  break;
966  }
967  case kD0toKpi:
968  {
969  if(isSelected==1 || isSelected==3) {
970  double sparsearray[9] = {invMass[0],ptD,vnfunc,phifunc1,phifunc2,phiD,evCentr,static_cast<double>(tracklets),candpercqn};
971  // cout<< isBDT;
972  if (((isBDT/1000000)%10)==1)fHistMassPtPhiqnCentr->Fill(sparsearray);
973  if (((isBDT/100000)%10)==1)fHistMassPtPhiqnCentr_1->Fill(sparsearray);
974  if (((isBDT/10000)%10)==1)fHistMassPtPhiqnCentr_2->Fill(sparsearray);
975  if (((isBDT/1000)%10)==1)fHistMassPtPhiqnCentr_3->Fill(sparsearray);
976  if (((isBDT/100)%10)==1)fHistMassPtPhiqnCentr_4->Fill(sparsearray);
977  if (((isBDT/10)%10)==1)fHistMassPtPhiqnCentr_5->Fill(sparsearray);
978 
979  }
980  if(isSelected==2 || isSelected==3) {
981  // cout<< isBDT;
982  double sparsearray[9] = {invMass[1],ptD,vnfunc,phifunc1,phifunc2,phiD,evCentr,static_cast<double>(tracklets),candpercqn};
983  if (((isBDT/1000000)%10)==1)fHistMassPtPhiqnCentr->Fill(sparsearray);
984  if (((isBDT/100000)%10)==1)fHistMassPtPhiqnCentr_1->Fill(sparsearray);
985  if (((isBDT/10000)%10)==1)fHistMassPtPhiqnCentr_2->Fill(sparsearray);
986  if (((isBDT/1000)%10)==1)fHistMassPtPhiqnCentr_3->Fill(sparsearray);
987  if (((isBDT/100)%10)==1)fHistMassPtPhiqnCentr_4->Fill(sparsearray);
988  if (((isBDT/10)%10)==1)fHistMassPtPhiqnCentr_5->Fill(sparsearray);
989  }
990  break;
991  }
992  case kDstartoKpipi:
993  {
994  double sparsearray[9] = {invMass[0],ptD,vnfunc,phifunc1,phifunc2,phiD,evCentr,static_cast<double>(tracklets),candpercqn};
995  fHistMassPtPhiqnCentr->Fill(sparsearray);
996  break;
997  }
998  case kDstoKKpi:
999  {
1000  if(isSelected&4) {
1001  double sparsearray[9] = {invMass[0],ptD,vnfunc,phifunc1,phifunc2,phiD,evCentr,static_cast<double>(tracklets),candpercqn};
1002  fHistMassPtPhiqnCentr->Fill(sparsearray);
1003  }
1004  if(isSelected&8) {
1005  double sparsearray[9] = {invMass[1],ptD,vnfunc,phifunc1,phifunc2,phiD,evCentr,static_cast<double>(tracklets),candpercqn};
1006  fHistMassPtPhiqnCentr->Fill(sparsearray);
1007  }
1008  break;
1009  }
1010  }
1011  }
1012 
1013  fHistCandVsCent->Fill(evCentr,nSelCand);
1014 
1015  trackstoremove.clear();
1016  isSelByPrevLoop.clear();
1017 
1018  fDaughterTracks.Clear();
1019 
1020  delete vHF;
1021  vHF = nullptr;
1022  if(!isHandlerFound) { // if not found in the tender task, allocated memory with new --> to be deleted
1023  delete HFQnVectorHandler;
1024  HFQnVectorHandler = nullptr;
1025  }
1026 
1027  PostData(1,fOutput);
1028  // PostData(3,fListBDTNtuple);
1029  return;
1030 }
1031 
1032 //________________________________________________________________________
1034 {
1035  // Interface method to set some specific cases
1036 
1037  switch (detconf) {
1038  case kTPC:
1042  break;
1043  case kTPCVZERO:
1045  fSubEvDetA=kV0A;
1046  fSubEvDetB=kV0C;
1047  break;
1048  case kVZERO:
1052  break;
1053  case kVZEROA:
1054  fEvPlaneDet=kV0A;
1055  fSubEvDetA=kV0C;
1057  break;
1058  case kVZEROC:
1059  fEvPlaneDet=kV0C;
1060  fSubEvDetA=kV0A;
1062  break;
1063  case kPosTPCVZERO:
1065  fSubEvDetA=kV0A;
1066  fSubEvDetB=kV0C;
1067  break;
1068  case kNegTPCVZERO:
1070  fSubEvDetA=kV0A;
1071  fSubEvDetB=kV0C;
1072  break;
1073  default:
1074  AliWarning("Case not defined -> use full V0");
1078  break;
1079  }
1080 }
1081 
1082 //________________________________________________________________________
1084 {
1085  // Sets the phi angle in the range [0,2*pi/harmonic]
1086 
1087  double result = phi;
1088  while(result < 0) {
1089  result = result + 2. * TMath::Pi() / fHarmonic;
1090  }
1091  while(result > 2.*TMath::Pi() / fHarmonic){
1092  result = result - 2. * TMath::Pi() / fHarmonic;
1093  }
1094  return result;
1095 }
1096 
1097 //________________________________________________________________________
1099 {
1100  // difference of subevents reaction plane angle cannot be bigger than pi / n
1101 
1102  double delta = psi1 - psi2;
1103  if(TMath::Abs(delta) > TMath::Pi() / fHarmonic) {
1104  if(delta>0.) delta -= 2.*TMath::Pi() / fHarmonic;
1105  else delta += 2.*TMath::Pi() / fHarmonic;
1106  }
1107 
1108  return delta;
1109 }
1110 
1111 //________________________________________________________________________
1113 {
1114  //Calculates all the possible invariant masses for each candidate
1115  switch(fDecChannel) {
1116  case kDplustoKpipi:
1117  {
1118  nmasses=1;
1119  masses=new float[nmasses];
1120  unsigned int pdgdaughters[3] = {211,321,211};
1121  masses[0]=d->InvMass(3,pdgdaughters);
1122  }
1123  break;
1124  case kD0toKpi:
1125  {
1126  nmasses=2;
1127  masses=new float[nmasses];
1128  unsigned int pdgdaughtersD0[2]={211,321};//pi,K
1129  masses[0]=d->InvMass(2,pdgdaughtersD0); //D0
1130  unsigned int pdgdaughtersD0bar[2]={321,211};//K,pi
1131  masses[1]=d->InvMass(2,pdgdaughtersD0bar); //D0bar
1132  }
1133  break;
1134  case kDstartoKpipi:
1135  {
1136  nmasses=1;
1137  masses=new float[nmasses];
1138  masses[0]=((AliAODRecoCascadeHF*)d)->DeltaInvMass();
1139  }
1140  break;
1141  case kDstoKKpi:
1142  {
1143  nmasses=2;
1144  masses=new float[nmasses];
1145  unsigned int pdgdaughtersKKpi[3] = {321,321,211};
1146  unsigned int pdgdaughterspiKK[3] = {211,321,321};
1147  masses[0]=d->InvMass(3,pdgdaughtersKKpi);
1148  masses[1]=d->InvMass(3,pdgdaughterspiKK);
1149  }
1150  break;
1151  }
1152 }
1153 
1154 //________________________________________________________________________
1155 void AliAnalysisTaskSECharmHadronvnTMVA::GetMainQnVectorInfo(double &mainPsin, double &mainMultQn, double mainQn[2], double &SubAPsin, double &SubAMultQn, double SubAQn[2], double &SubBPsin, double &SubBMultQn, double SubBQn[2], AliHFQnVectorHandler* HFQnVectorHandler)
1156 {
1157  double QnFullTPC[2], QnPosTPC[2], QnNegTPC[2];
1158  double QnFullV0[2], QnV0A[2], QnV0C[2];
1159  double MultQnFullTPC = -1., MultQnPosTPC = -1., MultQnNegTPC = -1.;
1160  double MultQnFullV0 = -1., MultQnV0A = -1., MultQnV0C = -1.;
1161  double PsinFullTPC = -1., PsinPosTPC = -1., PsinNegTPC = -1.;
1162  double PsinFullV0 = -1., PsinV0A = -1., PsinV0C = -1.;
1163 
1164  //get the unnormalised Qn-vectors --> normalisation can be done in the task
1165  HFQnVectorHandler->GetUnNormQnVecTPC(QnFullTPC,QnPosTPC,QnNegTPC);
1166  HFQnVectorHandler->GetUnNormQnVecV0(QnFullV0,QnV0A,QnV0C);
1167 
1168  HFQnVectorHandler->GetMultQnVecTPC(MultQnFullTPC,MultQnPosTPC,MultQnNegTPC);
1169  HFQnVectorHandler->GetMultQnVecV0(MultQnFullV0,MultQnV0A,MultQnV0C);
1170 
1171  HFQnVectorHandler->GetEventPlaneAngleTPC(PsinFullTPC,PsinPosTPC,PsinNegTPC);
1172  HFQnVectorHandler->GetEventPlaneAngleV0(PsinFullV0,PsinV0A,PsinV0C);
1173 
1174  if(fEtaGapInTPCHalves>0.) {
1175  vector<AliAODTrack*> trackstoremove;
1176  for(int iTrack=0; iTrack<fAOD->GetNumberOfTracks(); iTrack++) {
1177  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTrack));
1178  if(TMath::Abs(track->Eta())<fEtaGapInTPCHalves/2)
1179  trackstoremove.push_back(track);
1180  }
1181  HFQnVectorHandler->RemoveTracksFromQnTPC(trackstoremove, QnFullTPC, QnPosTPC, QnNegTPC, MultQnFullTPC, MultQnPosTPC, MultQnNegTPC, true);
1182  PsinFullTPC = (TMath::Pi()+TMath::ATan2(-QnFullTPC[1],-QnFullTPC[0]))/2;
1183  PsinPosTPC = (TMath::Pi()+TMath::ATan2(-QnPosTPC[1],-QnPosTPC[0]))/2;
1184  PsinNegTPC = (TMath::Pi()+TMath::ATan2(-QnNegTPC[1],-QnNegTPC[0]))/2;
1185  }
1186 
1187  switch(fEvPlaneDet) {
1188  case kFullTPC:
1189  mainPsin = PsinFullTPC;
1190  mainMultQn = MultQnFullTPC;
1191  mainQn[0] = QnFullTPC[0];
1192  mainQn[1] = QnFullTPC[1];
1193  break;
1194  case kPosTPC:
1195  mainPsin = PsinPosTPC;
1196  mainMultQn = MultQnPosTPC;
1197  mainQn[0] = QnPosTPC[0];
1198  mainQn[1] = QnPosTPC[1];
1199  break;
1200  case kNegTPC:
1201  mainPsin = PsinNegTPC;
1202  mainMultQn = MultQnNegTPC;
1203  mainQn[0] = QnNegTPC[0];
1204  mainQn[1] = QnNegTPC[1];
1205  break;
1206  case kFullV0:
1207  mainPsin = PsinFullV0;
1208  mainMultQn = MultQnFullV0;
1209  mainQn[0] = QnFullV0[0];
1210  mainQn[1] = QnFullV0[1];
1211  break;
1212  case kV0A:
1213  mainPsin = PsinV0A;
1214  mainMultQn = MultQnV0A;
1215  mainQn[0] = QnV0A[0];
1216  mainQn[1] = QnV0A[1];
1217  break;
1218  case kV0C:
1219  mainPsin = PsinV0C;
1220  mainMultQn = MultQnV0C;
1221  mainQn[0] = QnV0C[0];
1222  mainQn[1] = QnV0C[1];
1223  break;
1224  }
1225 
1226  switch(fSubEvDetA) {
1227  case kFullTPC:
1228  SubAPsin = PsinFullTPC;
1229  SubAMultQn = MultQnFullTPC;
1230  SubAQn[0] = QnFullTPC[0];
1231  SubAQn[1] = QnFullTPC[1];
1232  break;
1233  case kPosTPC:
1234  SubAPsin = PsinPosTPC;
1235  SubAMultQn = MultQnPosTPC;
1236  SubAQn[0] = QnPosTPC[0];
1237  SubAQn[1] = QnPosTPC[1];
1238  break;
1239  case kNegTPC:
1240  SubAPsin = PsinNegTPC;
1241  SubAMultQn = MultQnNegTPC;
1242  SubAQn[0] = QnNegTPC[0];
1243  SubAQn[1] = QnNegTPC[1];
1244  break;
1245  case kFullV0:
1246  SubAPsin = PsinFullV0;
1247  SubAMultQn = MultQnFullV0;
1248  SubAQn[0] = QnFullV0[0];
1249  SubAQn[1] = QnFullV0[1];
1250  break;
1251  case kV0A:
1252  SubAPsin = PsinV0A;
1253  SubAMultQn = MultQnV0A;
1254  SubAQn[0] = QnV0A[0];
1255  SubAQn[1] = QnV0A[1];
1256  break;
1257  case kV0C:
1258  SubAPsin = PsinV0C;
1259  SubAMultQn = MultQnV0C;
1260  SubAQn[0] = QnV0C[0];
1261  SubAQn[1] = QnV0C[1];
1262  break;
1263  }
1264 
1265  switch(fSubEvDetB) {
1266  case kFullTPC:
1267  SubBPsin = PsinFullTPC;
1268  SubBMultQn = MultQnFullTPC;
1269  SubBQn[0] = QnFullTPC[0];
1270  SubBQn[1] = QnFullTPC[1];
1271  break;
1272  case kPosTPC:
1273  SubBPsin = PsinPosTPC;
1274  SubBMultQn = MultQnPosTPC;
1275  SubBQn[0] = QnPosTPC[0];
1276  SubBQn[1] = QnPosTPC[1];
1277  break;
1278  case kNegTPC:
1279  SubBPsin = PsinNegTPC;
1280  SubBMultQn = MultQnNegTPC;
1281  SubBQn[0] = QnNegTPC[0];
1282  SubBQn[1] = QnNegTPC[1];
1283  break;
1284  case kFullV0:
1285  SubBPsin = PsinFullV0;
1286  SubBMultQn = MultQnFullV0;
1287  SubBQn[0] = QnFullV0[0];
1288  SubBQn[1] = QnFullV0[1];
1289  break;
1290  case kV0A:
1291  SubBPsin = PsinV0A;
1292  SubBMultQn = MultQnV0A;
1293  SubBQn[0] = QnV0A[0];
1294  SubBQn[1] = QnV0A[1];
1295  break;
1296  case kV0C:
1297  SubBPsin = PsinV0C;
1298  SubBMultQn = MultQnV0C;
1299  SubBQn[0] = QnV0C[0];
1300  SubBQn[1] = QnV0C[1];
1301  break;
1302  }
1303 }
1304 
1305 //________________________________________________________________________
1306 void AliAnalysisTaskSECharmHadronvnTMVA::GetDaughterTracksToRemove(AliAODRecoDecayHF* d, int nDau, vector<AliAODTrack*> &trackstoremove)
1307 {
1308  if(fDecChannel!=kDstartoKpipi) {
1309  for(int iDau=0; iDau<nDau; iDau++) {
1310  AliAODTrack* dau = dynamic_cast<AliAODTrack*>(d->GetDaughter(iDau));
1311  trackstoremove.push_back(dau);
1312  }
1313  }
1314  else {
1315  AliAODRecoDecayHF2Prong *dauD0 = dynamic_cast<AliAODRecoCascadeHF*>(d)->Get2Prong();
1316  trackstoremove.push_back(dynamic_cast<AliAODTrack*>(dauD0->GetDaughter(0)));
1317  trackstoremove.push_back(dynamic_cast<AliAODTrack*>(dauD0->GetDaughter(1)));
1318  if(fRemoveSoftPion) trackstoremove.push_back(dynamic_cast<AliAODRecoCascadeHF*>(d)->GetBachelor());
1319  }
1320 }
1321 
1322 //________________________________________________________________________
1324  if(!d || !vHF || (fDecChannel==kDstartoKpipi && !dD0)) return false;
1325  //Preselection to speed up task
1326  TObjArray arrDauTracks(nDau);
1327  AliAODTrack *track = nullptr;
1328  if(fDecChannel!=kDstartoKpipi) {
1329  for(int iDau=0; iDau<nDau; iDau++){
1330  AliAODTrack *track = vHF->GetProng(fAOD,d,iDau);
1331  arrDauTracks.AddAt(track,iDau);
1332  }
1333  }
1334  else {
1335  for(int iDau=0; iDau<nDau; iDau++){
1336  if(iDau == 0)
1337  track=vHF->GetProng(fAOD,d,iDau); //soft pion
1338  else
1339  track=vHF->GetProng(fAOD,dD0,iDau-1); //D0 daughters
1340  arrDauTracks.AddAt(track,iDau);
1341  }
1342  }
1343  if(!fRDCuts->PreSelect(arrDauTracks)){
1344  return 0;
1345  }
1346 
1347  bool isSelBit=true;
1348  switch(fDecChannel) {
1349  case kDplustoKpipi:
1351  if(!isSelBit || !vHF->FillRecoCand(fAOD,(AliAODRecoDecayHF3Prong*)d)) return 0;
1352  break;
1353  case kD0toKpi:
1355  if(!isSelBit || !vHF->FillRecoCand(fAOD,(AliAODRecoDecayHF2Prong*)d)) return 0;
1356  break;
1357  case kDstartoKpipi:
1358  if(!vHF->FillRecoCasc(fAOD,((AliAODRecoCascadeHF*)d),true)) return 0;
1359  break;
1360  case kDstoKKpi:
1361  isSelBit = d->HasSelectionBit(AliRDHFCuts::kDsCuts);
1362  if(!isSelBit || !vHF->FillRecoCand(fAOD,(AliAODRecoDecayHF3Prong*)d)) return 0;
1363  break;
1364  }
1365 
1366  double ptD = d->Pt();
1367  double yD = d->Y(absPdgMom);
1368  int ptbin=fRDCuts->PtBin(ptD);
1369  if(ptbin<0) {
1370  fHistNEvents->Fill(15);
1371  return 0;
1372  }
1373  bool isFidAcc = fRDCuts->IsInFiducialAcceptance(ptD,yD);
1374  if(!isFidAcc) return 0;
1375 
1376  int isSelected = fRDCuts->IsSelected(d,AliRDHFCuts::kAll,fAOD);
1377 
1378  //ML application
1379 
1380  if(!isSelected) return 0;
1381  if(fDecChannel==kDstoKKpi)
1382  if(!(isSelected&4) && !(isSelected&8)) return 0;
1383 
1384  return isSelected;
1385 }
1386 
1387 //________________________________________________________________________
1389 {
1390  // load splines from file
1391 
1392  TString listname[6] = {"SplineListq2TPC", "SplineListq2TPCPosEta", "SplineListq2TPCNegEta", "SplineListq2V0", "SplineListq2V0A", "SplineListq2V0C"};
1393 
1394  if (!gGrid) {
1395  TGrid::Connect("alien://");
1396  }
1397  TFile* splinesfile = TFile::Open(fqnSplineFileName.Data());
1398  if(!splinesfile) {
1399  AliFatal("File with splines for qn percentiles not found!");
1400  return false;
1401  }
1402 
1403  for(int iDet=0; iDet<6; iDet++) {
1404  fqnSplinesList[iDet] = (TList*)splinesfile->Get(listname[iDet].Data());
1405  if(!fqnSplinesList[iDet]) {
1406  AliFatal("TList with splines for qn percentiles not found in the spline file!");
1407  return false;
1408  }
1409  fqnSplinesList[iDet]->SetOwner();
1410  }
1411  splinesfile->Close();
1412 
1413  return true;
1414 }
1415 
1416 //________________________________________________________________________
1418  if(!vHF||!dD0||!fAOD) return 0;
1419  if(!isSelected) return 0;
1420  AliAODTrack *prong2 = vHF->GetProng(fAOD,dD0,0);
1421  AliAODTrack *prong3 = vHF->GetProng(fAOD,dD0,1);
1422  Double_t normIP[2];
1423  Double_t d0Prong[2];
1424  Double_t ptProng[2];
1425 
1426  Double_t invmassD0 = dD0->InvMassD0(); Double_t invmassD0bar = dD0->InvMassD0bar();
1427  Double_t ptB;
1428 
1429  Double_t cosPointingAngle;
1430  Double_t cosThetaStarD0 = 99;
1431  Double_t cosThetaStarD0bar = 99;
1432  Double_t diffIP[2], errdiffIP[2];
1433  dD0->Getd0MeasMinusExpProng(0,fAOD->GetMagneticField(),diffIP[0],errdiffIP[0]);
1434  dD0->Getd0MeasMinusExpProng(1,fAOD->GetMagneticField(),diffIP[1],errdiffIP[1]);
1435  if(!diffIP[0]||!errdiffIP[0])return 0;
1436  normIP[0]=diffIP[0]/errdiffIP[0];
1437  normIP[1]=diffIP[1]/errdiffIP[1];
1438 
1439  AliAODVertex *secvtx = dD0->GetSecondaryVtx();
1440  AliAODVertex *primvtx = dD0->GetPrimaryVtx();
1441  Double_t err2decaylength = secvtx->Error2DistanceXYToVertex(primvtx);
1442  Double_t lxy = dD0->AliAODRecoDecay::DecayLengthXY(primvtx);
1443  Bool_t isusepid = fRDCuts->GetIsUsePID();
1444  //peng
1445  d0Prong[0]=dD0->Getd0Prong(0); d0Prong[1]=dD0->Getd0Prong(1); //d0*d0 d0Prong[1]*d0Prong[0]
1446  cosPointingAngle = dD0->CosPointingAngle();
1447  ptProng[0]=dD0->Pt2Prong(0); ptProng[1]=dD0->Pt2Prong(1);
1448  cosThetaStarD0 = dD0->CosThetaStarD0();
1449  cosThetaStarD0bar = dD0->CosThetaStarD0bar();
1450  // if (part->Pt() > 10) return;
1451  // DCA part->GetDCA();
1452  Float_t tmp[22];
1453  tmp[8]= -99; // Invariant Mass
1454  tmp[18] = -99; // ptB (if accessible)
1455  tmp[19] = 0; // PDGCode
1456  tmp[20] = -99; // Rapidity YD0
1457  tmp[21] = -99; // Azimuthal Phi
1458 
1459  tmp[0] = dD0->Pt(); // ptD
1460  tmp[1] = normIP[0]; // Normalized d0N-<d0N> (topo1)
1461  tmp[2] = normIP[1]; // Normalized d0P-<d0P> (topo2)
1462  tmp[3] = lxy; // Decay length
1463  tmp[4] = lxy/TMath::Sqrt(err2decaylength); // Normalized decay length
1464  tmp[9] = d0Prong[0]*d0Prong[1]; // d0N*d0P
1465  tmp[10] = cosPointingAngle; // CosThetaPointing
1466  tmp[11] = dD0->GetDCA(); // DCAtracks
1467  tmp[12] = prong2->Pt(); // pt1
1468  tmp[13] = prong3->Pt(); // pt2
1469  tmp[14] = dD0->CosPointingAngleXY(); // CosThetaPointingXY
1470  tmp[15] = d0Prong[0]; // d01
1471  tmp[16] = d0Prong[1]; // d02
1472  tmp[20] = dD0->YD0();
1473  tmp[21] = dD0->Phi();
1474 
1475  // PID and Cuts
1476  if(isusepid)fRDCuts->SetUsePID(kFALSE);// if PID on, switch it off
1477  Int_t isCuts=fRDCuts->IsSelected(dD0,AliRDHFCuts::kAll,fAOD);
1478  if(isusepid)fRDCuts->SetUsePID(kTRUE);//if PID was on, switch it on
1479  Int_t isPid= fRDCuts->IsSelectedPID(dD0);
1480  tmp[5] = (Double_t)isCuts;
1481  tmp[6] = (Double_t)isPid;
1482 
1483  //~ std::vector<TString> vnameString; vnameString.resize(23);
1485  //~ vnameString.push_back("ptD"); vnameString.push_back("topo1"); vnameString.push_back("topo2"); vnameString.push_back("lxy"); vnameString.push_back("nlxy");
1486  //~ vnameString.push_back("iscut"); vnameString.push_back("ispid"); vnameString.push_back("type"); vnameString.push_back("mass"); vnameString.push_back("d0d0");
1487  //~ vnameString.push_back("cosp"); vnameString.push_back("dca"); vnameString.push_back("ptk"); vnameString.push_back("ptpi"); vnameString.push_back("cospxy");
1488  //~ vnameString.push_back("d0k"); vnameString.push_back("d0pi"); vnameString.push_back("cosstar"); vnameString.push_back("ptB"); vnameString.push_back("pdgcode");
1489  //~ vnameString.push_back("YD0"); vnameString.push_back("phi");
1490  // float BDTmatrix[13][2]={{0.11,0.13},/* 1<pt<2*/
1491  // {0.08,0.13},/* 2<pt<3*/
1492  // {0.11,0.15},/* 3<pt<4 */
1493  // {0.10,0.16},/* 4<pt<5 */
1494  // {0.07,0.16},/* 5<pt<6 */
1495  // {0.08,0.12},/* 6<pt<7 */
1496  // {0.08,0.08},/* 7<pt<8 */
1497  // {0.05,0.08},/* 8<pt<10 */
1498  // {0.05,0.06},/* 10<pt<12 */
1499  // {0.04,0.09},/* 12<pt<16 */
1500  // {0.03,0.05},/* 16<pt<24 */
1501  // {0.01,0.04},/* 24<pt<36 */
1502  // {0.01,0.04}};/* 36<pt<50 */
1503  Double_t BDTmatrix_1[13][2]={{-0.01,0.095},/* 1<pt<2*/
1504  {-0.01,0.16},/* 2<pt<3*/
1505  {-0.01,0.155},/* 3<pt<4 */
1506  {-0.01,0.165},/* 4<pt<5 */
1507  {-0.01,0.145},/* 5<pt<6 */
1508  {-0.01,0.14},/* 6<pt<7 */
1509  {-0.01,0.105},/* 7<pt<8 */
1510  {-0.01,0.11},/* 8<pt<10 */
1511  {0,0.065},/* 10<pt<12 */
1512  {-0.01,0.08},/* 12<pt<16 */
1513  {-0.01,0.095},/* 16<pt<24 */
1514  {-0.01,0.1},/* 24<pt<36 */
1515  {-0.01,0.1}};/* 36<pt<50 */
1516 
1517  Double_t BDTmatrix_2[13][2]={{0.01,0.095},/* 1<pt<2*/
1518  {0.01,0.16},/* 2<pt<3*/
1519  {0.01,0.16},/* 3<pt<4 */
1520  {0.01,0.16},/* 4<pt<5 */
1521  {0.01,0.15},/* 5<pt<6 */
1522  {0.01,0.14},/* 6<pt<7 */
1523  {0.01,0.11},/* 7<pt<8 */
1524  {0.01,0.115},/* 8<pt<10 */
1525  {0.01,0.065},/* 10<pt<12 */
1526  {0,0.08},/* 12<pt<16 */
1527  {0,0.095},/* 16<pt<24 */
1528  {0,0.1},/* 24<pt<36 */
1529  {0,0.1}};/* 36<pt<50 */
1530  Double_t BDTmatrix_3[13][2]={{0.03,0.135},/* 1<pt<2*/
1531  {0.03,0.16},/* 2<pt<3*/
1532  {0.02,0.16},/* 3<pt<4 */
1533  {0.02,0.145},/* 4<pt<5 */
1534  {0.02,0.145},/* 5<pt<6 */
1535  {0.03,0.12},/* 6<pt<7 */
1536  {0.02,0.095},/* 7<pt<8 */
1537  {0.02,0.09},/* 8<pt<10 */
1538  {0.02,0.05},/* 10<pt<12 */
1539  {0.015,0.095},/* 12<pt<16 */
1540  {0.015,0.09},/* 16<pt<24 */
1541  {0.005,0.08},/* 24<pt<36 */
1542  {0.005,0.08}};/* 36<pt<50 */
1543  Double_t BDTmatrix_4[13][2]={{0.05,0.1},/* 1<pt<2*/
1544  {0.05,0.16},/* 2<pt<3*/
1545  {0.04,0.16},/* 3<pt<4 */
1546  {0.04,0.145},/* 4<pt<5 */
1547  {0.04,0.145},/* 5<pt<6 */
1548  {0.05,0.115},/* 6<pt<7 */
1549  {0.04,0.1},/* 7<pt<8 */
1550  {0.04,0.09},/* 8<pt<10 */
1551  {0.04,0.065},/* 10<pt<12 */
1552  {0.025,0.095},/* 12<pt<16 */
1553  {0.025,0.09},/* 16<pt<24 */
1554  {0.015,0.1},/* 24<pt<36 */
1555  {0.015,0.1}};/* 36<pt<50 */
1556  Double_t BDTmatrix_5[13][2]={{0.07,0.14},/* 1<pt<2*/
1557  {0.07,0.135},/* 2<pt<3*/
1558  {0.06,0.16},/* 3<pt<4 */
1559  {0.06,0.155},/* 4<pt<5 */
1560  {0.05,0.16},/* 5<pt<6 */
1561  {0.06,0.12},/* 6<pt<7 */
1562  {0.05,0.11},/* 7<pt<8 */
1563  {0.05,0.08},/* 8<pt<10 */
1564  {0.055,0.055},/* 10<pt<12 */
1565  {0.035,0.095},/* 12<pt<16 */
1566  {0.03,0.095},/* 16<pt<24 */
1567  {0.02,0.1},/* 24<pt<36 */
1568  {0.02,0.1}};/* 36<pt<50 */
1569  Double_t BDTmatrix_6[13][2]={{0.085,0.135},/* 1<pt<2*/
1570  {0.08,0.135},/* 2<pt<3*/
1571  {0.08,0.16},/* 3<pt<4 */
1572  {0.08,0.155},/* 4<pt<5 */
1573  {0.07,0.15},/* 5<pt<6 */
1574  {0.07,0.12},/* 6<pt<7 */
1575  {0.07,0.13},/* 7<pt<8 */
1576  {0.07,0.075},/* 8<pt<10 */
1577  {0.065,0.055},/* 10<pt<12 */
1578  {0.045,0.095},/* 12<pt<16 */
1579  {0.035,0.1},/* 16<pt<24 */
1580  {0.03,0.1},/* 24<pt<36 */
1581  {0.03,0.1}};/* 36<pt<50 */
1582 
1583 
1584 
1585 
1586  // if(!isSelected) return 0;
1587  // printf("isSelected = %d\n",isSelected);
1588  int isBDTSelected = 0;
1589  std::vector<Double_t> BDTClsVar;// BDT cls input
1590  BDTClsVar.resize(10);
1591  // Data fill this
1592  Int_t thisptbin = fRDCuts->PtBin(tmp[0]);
1593  if(thisptbin<0) return 0;
1594  Float_t *ptbin = fRDCuts->GetPtBinLimits();
1595  TString ptstring = Form("_%.0f_%.0f",ptbin[thisptbin],ptbin[thisptbin+1]);
1596  // printf("ptstring = _%.0f_%.0f\n",ptbin[thisptbin],ptbin[thisptbin+1]);
1597  Double_t bdt1resp1 = -3; Double_t bdt2resp1_0 = -3; Double_t bdt2resp1_1 = -3; Double_t bdt2resp1_2 = -3;
1598  Double_t bdt1resp = -3; Double_t bdt2resp_0 = -3; Double_t bdt2resp_1 = -3; Double_t bdt2resp_2 = -3;
1599  if(isSelected==1 || isSelected==3){
1600  tmp[7] = 1; tmp[8] = invmassD0; tmp[17] = cosThetaStarD0;
1601  if(tmp[8]>2.12||tmp[8]<1.65) return 0;
1602 
1603  // Link variables to be used as classifier
1604  BDTClsVar[0] = tmp[1]; BDTClsVar[1] = tmp[2]; BDTClsVar[2] = tmp[4]; BDTClsVar[3] = tmp[9]; BDTClsVar[4] = tmp[10];
1605  BDTClsVar[5] = tmp[11]; BDTClsVar[6] = tmp[14]; BDTClsVar[7] = tmp[15]; BDTClsVar[8] = tmp[16]; BDTClsVar[9] = tmp[17];
1606 
1607 
1608  AliRDHFBDT *thisbdt1 = (AliRDHFBDT*)fListRDHFBDT->FindObject(Form("_BDT1%s",ptstring.Data()));
1609  AliRDHFBDT *thisbdt2ll = (AliRDHFBDT*)fListRDHFBDT->FindObject(Form("_BDT2%s_0",ptstring.Data()));
1610  AliRDHFBDT *thisbdt2mm = (AliRDHFBDT*)fListRDHFBDT->FindObject(Form("_BDT2%s_1",ptstring.Data()));
1611  AliRDHFBDT *thisbdt2hh = (AliRDHFBDT*)fListRDHFBDT->FindObject(Form("_BDT2%s_2",ptstring.Data()));
1612 
1613  // Float_t bdt1resp1; Float_t bdt2resp1;
1614  bdt1resp1 = thisbdt1->GetResponse(BDTClsVar);
1615  bdt2resp1_0 = thisbdt2ll->GetResponse(BDTClsVar);
1616  bdt2resp1_1 = thisbdt2mm->GetResponse(BDTClsVar);
1617  bdt2resp1_2 = thisbdt2hh->GetResponse(BDTClsVar);
1618  // cout<< bdt1resp1;
1619  // cout<< bdt2resp1;
1620  // printf("bdt1resp1 = %lf\n",thisbdt1);
1621  // printf("bdt1resp2 = %lf\n",thisbdt2ll);
1622  //~ bdt1resp1=0; bdt2resp1=0; bdt2resp[1]=0; bdt2resp[2]=0; bdt2resp[3]=0; bdt2resp[4]=0; bdt2resp[5]=0;
1623  // BDT Responses ready
1624 
1625 
1626  }
1627  if (isSelected>1){
1628  tmp[7] = 2; tmp[8] = invmassD0bar; tmp[17] = cosThetaStarD0bar;
1629 
1630  if(tmp[8]>2.12||tmp[8]<1.65) return 0;
1631  // Link variables to be used as classifier
1632  BDTClsVar[0] = tmp[1]; BDTClsVar[1] = tmp[2]; BDTClsVar[2] = tmp[4]; BDTClsVar[3] = tmp[9]; BDTClsVar[4] = tmp[10];
1633  BDTClsVar[5] = tmp[11]; BDTClsVar[6] = tmp[14]; BDTClsVar[7] = tmp[15]; BDTClsVar[8] = tmp[16]; BDTClsVar[9] = tmp[17];
1634 
1635  // Data application
1636  AliRDHFBDT *thisbdt1 = (AliRDHFBDT*)fListRDHFBDT->FindObject(Form("_BDT1%s",ptstring.Data()));
1637  AliRDHFBDT *thisbdt2ll = (AliRDHFBDT*)fListRDHFBDT->FindObject(Form("_BDT2%s_0",ptstring.Data()));
1638  AliRDHFBDT *thisbdt2mm = (AliRDHFBDT*)fListRDHFBDT->FindObject(Form("_BDT2%s_1",ptstring.Data()));
1639  AliRDHFBDT *thisbdt2hh = (AliRDHFBDT*)fListRDHFBDT->FindObject(Form("_BDT2%s_2",ptstring.Data()));
1640 
1641  // thisbdt2ll->Print();
1642  // Float_t bdt1resp; Float_t bdt2resp;
1643  bdt1resp = thisbdt1->GetResponse(BDTClsVar);
1644  bdt2resp_0 = thisbdt2ll->GetResponse(BDTClsVar);
1645  bdt2resp_1 = thisbdt2mm->GetResponse(BDTClsVar);
1646  bdt2resp_2 = thisbdt2hh->GetResponse(BDTClsVar);
1647  // cout<< bdt1resp;
1648  // cout<< bdt2resp;
1649  // printf("bdt1resp1 = %lf\n",thisbdt1);
1650  // printf("bdt1resp2 = %lf\n",thisbdt2ll);
1651  //~ bdt1resp=0; bdt2resp[0]=0; bdt2resp[1]=0; bdt2resp[2]=0; bdt2resp[3]=0; bdt2resp[4]=0; bdt2resp[5]=0;
1652  // BDT Responses ready
1653 
1654 
1655  }
1656  // printf("thisptbin = %d\n",thisptbin);
1657  // printf("BDTmatrix1 = %f\n",fBDT1Cut[thisptbin]);
1658  // printf("BDTmatrix2 = %f\n",fBDT2Cut[thisptbin]);
1659  int rep1 = 0; int rep2 = 0; int rep3 = 0; int rep4 = 0; int rep5 = 0; int rep6 = 0;
1660  if((bdt1resp1>BDTmatrix_1[thisptbin][0] && bdt2resp1_0>BDTmatrix_1[thisptbin][1])||(bdt1resp>BDTmatrix_1[thisptbin][0]&& bdt2resp_0>BDTmatrix_1[thisptbin][1])) rep1 = 1;
1661  if((bdt1resp1>BDTmatrix_2[thisptbin][0] && bdt2resp1_0>BDTmatrix_2[thisptbin][1])||(bdt1resp>BDTmatrix_2[thisptbin][0]&& bdt2resp_0>BDTmatrix_2[thisptbin][1])) rep2 = 1;
1662  if((bdt1resp1>BDTmatrix_3[thisptbin][0] && bdt2resp1_1>BDTmatrix_3[thisptbin][1])||(bdt1resp>BDTmatrix_3[thisptbin][0]&& bdt2resp_1>BDTmatrix_3[thisptbin][1])) rep3 = 1;
1663  if((bdt1resp1>BDTmatrix_4[thisptbin][0] && bdt2resp1_1>BDTmatrix_4[thisptbin][1])||(bdt1resp>BDTmatrix_4[thisptbin][0]&& bdt2resp_1>BDTmatrix_4[thisptbin][1])) rep4 = 1;
1664  if((bdt1resp1>BDTmatrix_5[thisptbin][0] && bdt2resp1_2>BDTmatrix_5[thisptbin][1])||(bdt1resp>BDTmatrix_5[thisptbin][0]&& bdt2resp_2>BDTmatrix_5[thisptbin][1])) rep5 = 1;
1665  if((bdt1resp1>BDTmatrix_6[thisptbin][0] && bdt2resp1_2>BDTmatrix_6[thisptbin][1])||(bdt1resp>BDTmatrix_6[thisptbin][0]&& bdt2resp_2>BDTmatrix_6[thisptbin][1])) rep6 = 1;
1666  isBDTSelected = 10*rep6+100*rep5+1000*rep4+10000*rep3+100000*rep2+1000000*rep1;
1667  return isBDTSelected;
1668 
1669 
1670 }
Bool_t IsEventRejectedDueToCentrality() const
Definition: AliRDHFCuts.h:361
double fEtaGapInTPCHalves
number of bins in the mass histograms
Int_t pdg
AliAODTrack * GetProng(AliVEvent *event, AliAODRecoDecayHF *rd, Int_t iprong)
Bool_t IsEventRejectedDueToZVertexOutsideFiducialRegion() const
Definition: AliRDHFCuts.h:355
THnSparseF * fHistMassPtPhiqnCentr_5
! THnSparse for the analysis of vn
TH3F * fHistEPResolVsCentrVsqn[3]
! histos needed to compute EP resolution vs. centrality vs. qn
Bool_t IsEventRejectedDueToNotRecoVertex() const
Definition: AliRDHFCuts.h:346
Int_t PtBin(Float_t pt) const
Definition: AliRDHFCuts.h:330
virtual Int_t PreSelect(TObjArray)
Definition: AliRDHFCuts.h:314
double Double_t
Definition: External.C:58
TH3F * fHistEvPlaneQncorr[6]
! histograms for EP angle vs. centrality vs. qn
Definition: External.C:260
Definition: External.C:236
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
Bool_t HasSelectionBit(Int_t i) const
TList * fqnSplinesList[6]
flag to know if splines loaded with task or read from tender task
Double_t mass
TH1F * fHistCentrality[3]
! hist. for cent distr (all,sel ev,out of cent)
TH1F * fHistNEvents
! histogram send on output slot 1
static Int_t CheckMatchingAODdeltaAODevents()
Bool_t IsEventRejectedDueToVertexContributors() const
Definition: AliRDHFCuts.h:349
TObjArray fDaughterTracks
lists of splines used to compute the qn percentile
Double_t CosPointingAngleXY() const
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
THnSparseF * fHistMassPtPhiqnCentr_4
! THnSparse for the analysis of vn
void CalculateInvMasses(AliAODRecoDecayHF *d, float *&masses, int &nmasses)
int fFlowMethod
AODB file name for calibrations (if Qn-framework not used)
Bool_t FillRecoCasc(AliVEvent *event, AliAODRecoCascadeHF *rc, Bool_t isDStar, Bool_t recoSecVtx=kFALSE)
virtual Int_t IsSelectedPID(AliAODRecoDecayHF *)
Definition: AliRDHFCuts.h:321
int ProcessBDT(AliAODEvent *fAOD, AliAODRecoDecayHF2Prong *dD0, int isSelected, AliAnalysisVertexingHF *vHF)
Class for cuts on AOD reconstructed D+->Kpipi.
static Int_t GetNumberOfTrackletsInEtaRange(AliAODEvent *ev, Double_t mineta, Double_t maxeta)
TH2F * fHistqnVsCentrPercCalib[6]
! histograms for qn percentile calibration with qn fine binning
int Int_t
Definition: External.C:63
void SetMinCentrality(Float_t minCentrality=0.)
Definition: AliRDHFCuts.h:54
Int_t GetIsFilled() const
Double_t CosThetaStarD0bar() const
angle of K
float Float_t
Definition: External.C:68
const Double_t ptmax
THnSparseF * fHistMassPtPhiqnCentr
! THnSparse for the analysis of vn
bool fPercentileqn
max value for the scalar product histograms
TList * fListRDHFBDT
fraction of tracks to keep in qn with random downsampling
bool fLoadedSplines
qn spline file name for qn percentile calibrations
const Double_t ptmin
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:275
UShort_t GetProngID(Int_t ip) const
Bool_t IsEventRejectedDueToPileup() const
Definition: AliRDHFCuts.h:358
bool fEnableDownsamplqn
flag to enable removal of soft pion too (only D*)
TH3F * fHistNtrklVsqnVsCentr
! histo of Ntracklets vs. qn vs. centrality
TH2F * fHistCandVsCent
! hist. for number of selected candidates vs. centrality
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
void SetMaxCentrality(Float_t maxCentrality=100.)
Definition: AliRDHFCuts.h:55
THnSparseF * fHistMassPtPhiqnCentr_2
! THnSparse for the analysis of vn
TString fqnSplineFileName
flag to replace qn with its percentile in the histograms
int IsCandidateSelected(AliAODRecoDecayHF *&d, int nDau, int absPdgMom, AliAnalysisVertexingHF *vHF, AliAODRecoDecayHF2Prong *dD0)
THnSparseF * fHistMassPtPhiqnCentr_3
! THnSparse for the analysis of vn
Bool_t IsEventSelected(AliVEvent *event)
int fRemoveDauFromqn
eta gap between TPC subevents
void SetUsePID(Bool_t flag=kTRUE)
Definition: AliRDHFCuts.h:221
Float_t * GetPtBinLimits() const
Definition: AliRDHFCuts.h:265
TString fTenderTaskName
cut values (saved in slot 3)
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:312
bool fRemoveSoftPion
flag to enable removal of D-meson daughter tracks from qn (1->remove single cand, 2->remove all cand ...
Bool_t GetIsUsePID() const
Definition: AliRDHFCuts.h:287
const char Option_t
Definition: External.C:48
AliAODVertex * GetPrimaryVtx() const
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:343
const Int_t nbins
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
THnSparseF * fHistMassPtPhiqnCentr_1
! THnSparse for the analysis of vn
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:338
int fNMassBins
upper inv mass limit for histos
float fUpmasslimit
lower inv mass limit for histos
TH3F * fHistPercqnVsqnVsCentr
! histo of qn percentile vs. qn vs. centrality
TString fOADBFileName
Qn-vector normalisation method.
Int_t nptbins
void GetMainQnVectorInfo(double &mainPsin, double &mainMultQn, double mainQn[2], double &SubAPsin, double &SubAMultQn, double SubAQn[2], double &SubBPsin, double &SubBMultQn, double SubBQn[2], AliHFQnVectorHandler *HFQnVectorHandler)
void GetDaughterTracksToRemove(AliAODRecoDecayHF *d, int nDau, vector< AliAODTrack * > &trackstoremove)
const Double_t phimin
float fMinCentr
name of tender task needed to get the calibrated Qn vectors
double fFracToKeepDownSamplqn
flag to enable random downsampling for qn