AliPhysics  608b256 (608b256)
AliRDHFCutsBPlustoD0Pi.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2010, 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: AliRDHFCutsBPlustoD0Pi.cxx $ */
17 
19 //
20 // Class for cuts on AOD reconstructed BPlus->D0Pi->KPiPi
21 //
22 //
23 // Author Lennart van Doremalen
24 // Utrecht University - l.v.r.vandoremalen@uu.nl
25 //
26 // Several AliPhysics classes have been used as a basis for this code
27 //
28 //
30 
31 #include <TDatabasePDG.h>
32 #include <Riostream.h>
34 #include "AliRDHFCutsBPlustoD0Pi.h"
35 #include "AliAODTrack.h"
36 #include "AliESDtrack.h"
37 #include "AliAODPid.h"
38 #include "AliTPCPIDResponse.h"
39 #include "AliAODVertex.h"
40 #include "AliESDVertex.h"
41 
42 using std::cout;
43 using std::endl;
44 
46 ClassImp(AliRDHFCutsBPlustoD0Pi);
48 
49 
50 
51 //--------------------------------------------------------------------------
53  AliRDHFCuts(name),
54  fTrackCutsSoftPi(0),
55  fMaxPtPid(9999.),
56  fTPCflag(999.),
57  fCircRadius(0.),
58  fGetCutInfo(0),
59 
60  fIsCutUsed(0x0),
61 
62  fnVarsD0forD0ptbin(0),
63  fnPtBinsD0forD0ptbin(1),
64  fGlobalIndexD0forD0ptbin(0),
65  fCutsRDD0forD0ptbin(0x0),
66  fnPtBinLimitsD0forD0ptbin(1),
67  fPtBinLimitsD0forD0ptbin(0x0),
68  fIsUpperCutD0forD0ptbin(0x0),
69  fIsCutUsedD0forD0ptbin(0x0),
70  fVarNamesD0forD0ptbin(0x0),
71 
72  fMinITSNclsD0FirstDaughter(0),
73  fMinTPCNclsD0FirstDaughter(0),
74  fUseITSRefitD0FirstDaughter(0),
75  fUseTPCRefitD0FirstDaughter(0),
76  fUseFilterBitD0FirstDaughter(0),
77  fFilterBitD0FirstDaughter(0),
78  fMinPtD0FirstDaughter(0),
79  fMaxAbsEtaD0FirstDaughter(999.),
80  fHardSelectionArrayITSD0FirstDaughter(),
81  fSoftSelectionArrayITSD0FirstDaughter(),
82  fNSoftITSCutD0FirstDaughter(0),
83 
84  fMinITSNclsD0SecondDaughter(0),
85  fMinTPCNclsD0SecondDaughter(0),
86  fUseITSRefitD0SecondDaughter(0),
87  fUseTPCRefitD0SecondDaughter(0),
88  fUseFilterBitD0SecondDaughter(0),
89  fFilterBitD0SecondDaughter(0),
90  fMinPtD0SecondDaughter(0),
91  fMaxAbsEtaD0SecondDaughter(999.),
92  fHardSelectionArrayITSD0SecondDaughter(),
93  fSoftSelectionArrayITSD0SecondDaughter(),
94  fNSoftITSCutD0SecondDaughter(0),
95 
96  fMinITSNclsBPlusPion(0),
97  fMinTPCNclsBPlusPion(0),
98  fUseITSRefitBPlusPion(0),
99  fUseTPCRefitBPlusPion(0),
100  fUseFilterBitBPlusPion(0),
101  fFilterBitBPlusPion(0),
102  fMinPtBPlusPion(0),
103  fMaxAbsEtaBPlusPion(999.),
104  fHardSelectionArrayITSBPlusPion(),
105  fSoftSelectionArrayITSBPlusPion(),
106  fNSoftITSCutBPlusPion(0),
107 
108  fMind0D0FirstDaughter(0),
109  fMind0D0SecondDaughter(0),
110  fMind0BPlusPion(0.),
111  fFiducialYCut(0.8),
112 
113  fnVariablesForCutOptimization(0),
114  fnCutsForOptimization(0),
115  fGlobalIndexCutOptimization(0),
116  fCutsRDForCutOptimization(0x0),
117  fIsUpperCutForCutOptimization(0x0),
118  fCutIndexForCutOptimization(0x0),
119  fSigmaForCutOptimization(0x0)
120 {
121  //
122  // Default Constructor
123  //
124 
125  // Main cut setup as function of BPlus pt bins
126  const Int_t nvars = 68;
127  SetNVars(nvars);
128 
129  TString varNames[nvars];
130  Int_t iterator = 0;
131 
132  //D0 cut variables
133  varNames[iterator++] = /*-00-*/ "inv. mass width[GeV]";
134  varNames[iterator++] = /*-01-*/ "delta mass width [GeV]"; //not used for D0
135  varNames[iterator++] = /*-02-*/ "pointing angle [Cos(theta)]";
136  varNames[iterator++] = /*-03-*/ "dca [cm]";
137  varNames[iterator++] = /*-04-*/ "Pt D0 [GeV/c]";
138  varNames[iterator++] = /*-05-*/ "Pt first daughter [GeV/c]";
139  varNames[iterator++] = /*-06-*/ "Pt second daughter [GeV/c]";
140  varNames[iterator++] = /*-07-*/ "d0 D0 [cm]";
141  varNames[iterator++] = /*-08-*/ "d0 first daughter [cm]";
142  varNames[iterator++] = /*-09-*/ "d0 second daughter [cm]";
143  varNames[iterator++] = /*-10-*/ "d0d0 [cm^2]";
144  varNames[iterator++] = /*-11-*/ "d0d0 XY [cm^2]";
145 
146  varNames[iterator++] = /*-12-*/ "angle between both daughters";
147  varNames[iterator++] = /*-13-*/ "angle mother with first daughter";
148  varNames[iterator++] = /*-14-*/ "angle mother with second daughter";
149  varNames[iterator++] = /*-15-*/ "cosThetaStar";
150  varNames[iterator++] = /*-16-*/ "vertexDistance";
151  varNames[iterator++] = /*-17-*/ "pseudoProperDecayTime";
152  varNames[iterator++] = /*-18-*/ "DecayTime";
153  varNames[iterator++] = /*-19-*/ "normalizedDecayTime";
154  varNames[iterator++] = /*-20-*/ "normDecayLength";
155  varNames[iterator++] = /*-21-*/ "topomatic first daughter";
156  varNames[iterator++] = /*-22-*/ "topomatic second daughter";
157  varNames[iterator++] = /*-23-*/ "topomatic max";
158  varNames[iterator++] = /*-24-*/ "topomatic min";
159 
160  varNames[iterator++] = /*-25-*/ "pointing angle XY [Cos(theta)]";
161  varNames[iterator++] = /*-26-*/ "vertex distance XY [cm]";
162  varNames[iterator++] = /*-27-*/ "normDecayLength XY";
163  varNames[iterator++] = /*-28-*/ "Chi2 per NDF vertex";
164 
165  varNames[iterator++] = /*-29-*/ "pointingAngleToBPlus";
166  varNames[iterator++] = /*-30-*/ "d0MotherToBPlus";
167  varNames[iterator++] = /*-31-*/ "d0FirstDaughterToBPlus";
168  varNames[iterator++] = /*-32-*/ "d0SecondDaughterToBPlus";
169  varNames[iterator++] = /*-33-*/ "impactProductToBPlus";
170  varNames[iterator++] = /*-34-*/ "impactProductXYToBPlus";
171  varNames[iterator++] = /*-35-*/ "normDecayLengthToBPlus";
172  varNames[iterator++] = /*-36-*/ "pseudoProperDecayTimeToBPlus";
173  varNames[iterator++] = /*-37-*/ "DecayTimeToBPlus";
174  varNames[iterator++] = /*-38-*/ "normalizedDecayTimeToBPlus";
175 
176  //BPlus cut variables
177  varNames[iterator++] = /*-39-*/ "inv. mass width[GeV]";
178  varNames[iterator++] = /*-40-*/ "delta mass width [GeV]";
179  varNames[iterator++] = /*-41-*/ "pointing angle [Cos(theta)]";
180  varNames[iterator++] = /*-42-*/ "dca [cm]";
181  varNames[iterator++] = /*-43-*/ "Pt BPlus [GeV/c]";
182  varNames[iterator++] = /*-44-*/ "Pt D0 [GeV/c]";
183  varNames[iterator++] = /*-45-*/ "Pt Pion [GeV/c]";
184  varNames[iterator++] = /*-46-*/ "d0 BPlus [cm]";
185  varNames[iterator++] = /*-47-*/ "d0 D0 [cm]";
186  varNames[iterator++] = /*-48-*/ "d0 Pion [cm]";
187  varNames[iterator++] = /*-49-*/ "d0d0 [cm^2]";
188  varNames[iterator++] = /*-50-*/ "d0d0 XY [cm^2]";
189 
190  varNames[iterator++] = /*-51-*/ "angle between both daughters";
191  varNames[iterator++] = /*-52-*/ "angle mother with first daughter";
192  varNames[iterator++] = /*-53-*/ "angle mother with second daughter";
193  varNames[iterator++] = /*-54-*/ "cosThetaStar";
194  varNames[iterator++] = /*-55-*/ "vertexDistance";
195  varNames[iterator++] = /*-56-*/ "pseudoProperDecayTime";
196  varNames[iterator++] = /*-57-*/ "DecayTime";
197  varNames[iterator++] = /*-58-*/ "normalizedDecayTime";
198  varNames[iterator++] = /*-59-*/ "normDecayLength";
199  varNames[iterator++] = /*-60-*/ "topomatic first daughter";
200  varNames[iterator++] = /*-61-*/ "topomatic second daughter";
201  varNames[iterator++] = /*-62-*/ "topomatic max";
202  varNames[iterator++] = /*-63-*/ "topomatic min";
203 
204  varNames[iterator++] = /*-64-*/ "pointing angle XY [Cos(theta)]";
205  varNames[iterator++] = /*-65-*/ "vertex distance XY [cm]";
206  varNames[iterator++] = /*-66-*/ "normDecayLength XY";
207  varNames[iterator++] = /*-67-*/ "Chi2 per NDF vertex";
208 
209  Bool_t isUpperCut[nvars] = {0};
210 
211  SetVarNames(nvars, varNames, isUpperCut);
212 
213  Float_t limits[2] = {0, 999999999.};
214  SetPtBins(2, limits);
215 
216 
217  //
218  // Initialization of D0 cuts for D0 pt bins
219  //
220 
221  const Int_t nvarsD0forD0ptbin = 29;
222  SetNVarsD0forD0ptbin(nvarsD0forD0ptbin);
223 
224  TString varNamesD0forD0ptbin[nvarsD0forD0ptbin];
225  iterator = 0;
226 
227  //D0 cut variables
228  varNamesD0forD0ptbin[iterator++] = /*-00-*/ "inv. mass width[GeV]";
229  varNamesD0forD0ptbin[iterator++] = /*-01-*/ "delta mass width [GeV]"; //not used for D0
230  varNamesD0forD0ptbin[iterator++] = /*-02-*/ "pointing angle [Cos(theta)]";
231  varNamesD0forD0ptbin[iterator++] = /*-03-*/ "dca [cm]";
232  varNamesD0forD0ptbin[iterator++] = /*-04-*/ "Pt D0 [GeV/c]";
233  varNamesD0forD0ptbin[iterator++] = /*-05-*/ "Pt first daughter [GeV/c]";
234  varNamesD0forD0ptbin[iterator++] = /*-06-*/ "Pt second daughter [GeV/c]";
235  varNamesD0forD0ptbin[iterator++] = /*-07-*/ "d0 D0 [cm]";
236  varNamesD0forD0ptbin[iterator++] = /*-08-*/ "d0 first daughter [cm]";
237  varNamesD0forD0ptbin[iterator++] = /*-09-*/ "d0 second daughter [cm]";
238  varNamesD0forD0ptbin[iterator++] = /*-10-*/ "d0d0 [cm^2]";
239  varNamesD0forD0ptbin[iterator++] = /*-11-*/ "d0d0 XY [cm^2]";
240 
241  varNamesD0forD0ptbin[iterator++] = /*-12-*/ "angle between both daughters";
242  varNamesD0forD0ptbin[iterator++] = /*-13-*/ "angle mother with first daughter";
243  varNamesD0forD0ptbin[iterator++] = /*-14-*/ "angle mother with second daughter";
244  varNamesD0forD0ptbin[iterator++] = /*-15-*/ "cosThetaStar";
245  varNamesD0forD0ptbin[iterator++] = /*-16-*/ "vertexDistance";
246  varNamesD0forD0ptbin[iterator++] = /*-17-*/ "pseudoProperDecayTime";
247  varNamesD0forD0ptbin[iterator++] = /*-18-*/ "DecayTime";
248  varNamesD0forD0ptbin[iterator++] = /*-19-*/ "normalizedDecayTime";
249  varNamesD0forD0ptbin[iterator++] = /*-20-*/ "normDecayLength";
250  varNamesD0forD0ptbin[iterator++] = /*-21-*/ "topomatic first daughter";
251  varNamesD0forD0ptbin[iterator++] = /*-22-*/ "topomatic second daughter";
252  varNamesD0forD0ptbin[iterator++] = /*-23-*/ "topomatic max";
253  varNamesD0forD0ptbin[iterator++] = /*-24-*/ "topomatic min";
254 
255  varNamesD0forD0ptbin[iterator++] = /*-25-*/ "pointing angle XY [Cos(theta)]";
256  varNamesD0forD0ptbin[iterator++] = /*-26-*/ "vertex distance XY [cm]";
257  varNamesD0forD0ptbin[iterator++] = /*-27-*/ "normDecayLength XY";
258  varNamesD0forD0ptbin[iterator++] = /*-28-*/ "Chi2 per NDF vertex";
259 
260  Bool_t isUpperCutD0forD0ptbin[nvarsD0forD0ptbin] = {0};
261 
262  SetVarNamesD0forD0ptbin(nvarsD0forD0ptbin, varNamesD0forD0ptbin, isUpperCutD0forD0ptbin);
263 
264  Float_t limitsD0forD0ptbin[2] = {0, 999999999.};
265  SetPtBinsD0forD0ptbin(2, limitsD0forD0ptbin);
266 
267  Bool_t forOpt[16] = {0}; //not yet used for BPlus analysis
268  SetVarsForOpt(16, forOpt);
269 
270 }
271 //--------------------------------------------------------------------------
273  AliRDHFCuts(source),
274  fTrackCutsSoftPi(0),
275  fMaxPtPid(source.fMaxPtPid),
276  fTPCflag(source.fTPCflag),
277  fCircRadius(source.fCircRadius),
278  fGetCutInfo(source.fGetCutInfo),
279 
280  fIsCutUsed(0x0),
281 
285  fCutsRDD0forD0ptbin(0x0),
291 
303 
315 
327 
332 
340 {
341  //
342  // Copy constructor
343  //
344 
346 
349  if (source.fIsCutUsed)
350  {
351  if (fIsCutUsed) {
352  delete [] fIsCutUsed;
353  fIsCutUsed = nullptr;
354  }
355  fIsCutUsed = new Bool_t[(source.GetNPtBins()) * (source.GetNVars())];
356 
357  for (Int_t i = 0; i < source.fnVars; ++i)
358  {
359  for (Int_t j = 0; j < source.fnPtBins; j++)
360  {
361  Bool_t bUse = source.GetIsCutUsed(i, j);
362  SetIsCutUsed(i, j, bUse);
363  }
364  }
365  }
366  if (source.fIsCutUsedD0forD0ptbin)
367  {
369  delete [] fIsCutUsedD0forD0ptbin;
370  fIsCutUsedD0forD0ptbin = nullptr;
371  }
373  for (Int_t i = 0; i < source.fnVarsD0forD0ptbin; ++i)
374  {
375  for (Int_t j = 0; j < source.fnPtBinsD0forD0ptbin; j++)
376  {
377  Bool_t bUse = source.GetIsCutUsedD0forD0ptbin(i, j);
378  SetIsCutUsedD0forD0ptbin(i, j, bUse);
379  }
380  }
381  }
384  {
388  }
390  for (Int_t i = 0; i < source.fnVariablesForCutOptimization; i++)
391  {
392  Bool_t bUpperCut = source.GetIsUpperCutForCutOptimization(i);
393  SetIsUpperCutForCutOptimization(i, bUpperCut);
394  }
395  }
396  if (source.fCutIndexForCutOptimization)
397  {
399  delete [] fCutIndexForCutOptimization;
400  fCutIndexForCutOptimization = nullptr;
401  }
403  for (Int_t i = 0; i < source.fnVariablesForCutOptimization; i++)
404  {
405  Int_t nCutIndex = source.GetCutIndexForCutOptimization(i);
406  SetCutIndexForCutOptimization(i, nCutIndex);
407  }
408  }
409  if (source.fSigmaForCutOptimization)
410  {
412  delete [] fSigmaForCutOptimization;
413  fSigmaForCutOptimization = nullptr;
414  }
416  for (Int_t i = 0; i < source.fnPtBins; i++)
417  {
418  Double_t binSigma = source.GetSigmaForCutOptimization(i);
419  SetSigmaForCutOptimization(binSigma, i);
420  }
421  }
424 
431 }
432 //--------------------------------------------------------------------------
434  //
435  // Default Destructor
436  //
437  if (fTrackCutsSoftPi) { delete fTrackCutsSoftPi; fTrackCutsSoftPi = nullptr;}
438  if (fIsCutUsed) { delete [] fIsCutUsed; fIsCutUsed = nullptr; }
439  if (fCutsRDD0forD0ptbin) { delete [] fCutsRDD0forD0ptbin; fCutsRDD0forD0ptbin = nullptr; }
448 }
449 //--------------------------------------------------------------------------
451 {
452  //
453  // assignment operator
454  //
455 
456  if (&source == this) return *this;
457 
458  AliRDHFCuts::operator=(source);
459 
460  if(source.GetTrackCutsSoftPi()) {
461  delete fTrackCutsSoftPi;
462  fTrackCutsSoftPi = new AliESDtrackCuts(*(source.GetTrackCutsSoftPi()));
463  }
464  fMaxPtPid = source.fMaxPtPid;
465  fTPCflag = source.fTPCflag;
466  fCircRadius = source.fCircRadius;
467  fGetCutInfo = source.fGetCutInfo;
502  fFiducialYCut = source.fFiducialYCut;
506 
509  if (source.fIsCutUsed)
510  {
511  if (fIsCutUsed) {
512  delete [] fIsCutUsed;
513  fIsCutUsed = nullptr;
514  }
515  fIsCutUsed = new Bool_t[(source.GetNPtBins()) * (source.GetNVars())];
516 
517  for (Int_t i = 0; i < source.fnVars; ++i)
518  {
519  for (Int_t j = 0; j < source.fnPtBins; j++)
520  {
521  Bool_t bUse = source.GetIsCutUsed(i, j);
522  SetIsCutUsed(i, j, bUse);
523  }
524  }
525  }
526  if (source.fIsCutUsedD0forD0ptbin)
527  {
529  delete [] fIsCutUsedD0forD0ptbin;
530  fIsCutUsedD0forD0ptbin = nullptr;
531  }
533  for (Int_t i = 0; i < source.fnVarsD0forD0ptbin; ++i)
534  {
535  for (Int_t j = 0; j < source.fnPtBinsD0forD0ptbin; j++)
536  {
537  Bool_t bUse = source.GetIsCutUsedD0forD0ptbin(i, j);
538  SetIsCutUsedD0forD0ptbin(i, j, bUse);
539  }
540  }
541  }
544  {
548  }
550  for (Int_t i = 0; i < source.fnVariablesForCutOptimization; i++)
551  {
552  Bool_t bUpperCut = source.GetIsUpperCutForCutOptimization(i);
553  SetIsUpperCutForCutOptimization(i, bUpperCut);
554  }
555  }
556  if (source.fCutIndexForCutOptimization)
557  {
559  delete [] fCutIndexForCutOptimization;
560  fCutIndexForCutOptimization = nullptr;
561  }
563  for (Int_t i = 0; i < source.fnVariablesForCutOptimization; i++)
564  {
565  Int_t nCutIndex = source.GetCutIndexForCutOptimization(i);
566  SetCutIndexForCutOptimization(i, nCutIndex);
567  }
568  }
569  if (source.fSigmaForCutOptimization)
570  {
572  delete [] fSigmaForCutOptimization;
573  fSigmaForCutOptimization = nullptr;
574  }
576  for (Int_t i = 0; i < source.fnPtBins; i++)
577  {
578  Double_t binSigma = source.GetSigmaForCutOptimization(i);
579  SetSigmaForCutOptimization(binSigma, i);
580  }
581  }
584 
591 
592  return *this;
593 }
594 //--------------------------------------------------------------------------
595 void AliRDHFCutsBPlustoD0Pi::GetCutVarsForOpt(AliAODRecoDecayHF* /*d*/, Float_t* /*vars*/, Int_t /*nvars*/, Int_t* /*pdgdaughters*/) {
596  // not yet used
597 
598  return;
599 }
600 //--------------------------------------------------------------------------
601 Int_t AliRDHFCutsBPlustoD0Pi::IsSelected(TObject* obj, Int_t selectionLevel, AliAODEvent* aod, Bool_t bCutArray[68]) {
602  //
603  // In this function we apply the selection cuts on the BPlus candidate and its daughters.
604  // The function returns 0 if the candidate is cut and is able to return information on which cuts the candidate passes.
605  //
606 
607  fIsSelectedCuts = 0;
608  fIsSelectedPID = 0;
609 
610  // The cuts used in this class have to be set in the maketfile.
611  if (!fCutsRD) {
612  cout << "Cut matrice not inizialized. Exit..." << endl;
613  return 0;
614  }
615 
616  AliAODRecoDecayHF2Prong* candidateBPlus = (AliAODRecoDecayHF2Prong*)obj;
617  if (!candidateBPlus) {
618  cout << "candidateBPlus null" << endl;
619  return 0;
620  }
621 
622  AliAODRecoDecayHF2Prong* candidateD0 = (AliAODRecoDecayHF2Prong*)candidateBPlus->GetDaughter(1);
623  if (!candidateD0) {
624  cout << "candidateD0 null" << endl;
625  return 0;
626  }
627 
628  AliAODTrack *candidatePion = (AliAODTrack*)candidateBPlus->GetDaughter(0);
629  if (!candidatePion) {
630  cout << "candidatePion null 1" << endl;
631  return 0;
632  }
633 
634  AliAODVertex * vertexBPlus = candidateBPlus->GetSecondaryVtx();
635  if (!vertexBPlus) {
636  cout << "vertexBPlus null" << endl;
637  return 0;
638  }
639 
640  AliAODVertex * primaryVertex = aod->GetPrimaryVertex();
641  if (!primaryVertex) {
642  cout << "primaryVertex null" << endl;
643  return 0;
644  }
645 
646  Int_t returnvalue = 1;
647  Bool_t bPassedCut = kFALSE;
648 
649  //get the magnetic field
650  Double_t bz = (Double_t)aod->GetMagneticField();
651 
652  // selection on candidate
653  if (selectionLevel == AliRDHFCuts::kAll ||
654  selectionLevel == AliRDHFCuts::kCandidate) {
655 
656  // We check to which pt bin the candidate belongs
657  Int_t ptbin = PtBin(candidateBPlus->Pt());
658  if (ptbin == -1) return -1;
659 
660  // We obtain the variable values in the section below
661  // D0Mass and BPlusmass
662  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
663  Double_t mBPlusPDG = TDatabasePDG::Instance()->GetParticle(521)->Mass();
664 
665  // delta mass PDG
666  Double_t deltaPDG = mBPlusPDG - mD0PDG;
667 
668  // Half width BPlus mass
669  UInt_t prongs[2];
670  prongs[0] = 211; prongs[1] = 421;
671  Double_t invMassBPlus = candidateBPlus->InvMass(2, prongs);
672  Double_t invMassDifference = TMath::Abs(mBPlusPDG - invMassBPlus);
673  Double_t invMassDelta = TMath::Abs(deltaPDG - (DeltaInvMassBPlusKpipi(candidateBPlus)));
674 
675  Double_t pointingAngle = candidateBPlus->CosPointingAngle();
676  Double_t dcaMother = candidateBPlus->GetDCA();
677  Double_t ptMother = candidateBPlus->Pt();
678  Double_t momentumMother = candidateBPlus->P();
679  Double_t ptD0 = candidateD0->Pt();
680  Double_t ptPion = candidatePion->Pt();
681 
682  AliExternalTrackParam motherTrack;
683  motherTrack.CopyFromVTrack(candidateBPlus);
684  Double_t d0z0[2], covd0z0[3], d0[2];
685  motherTrack.PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
686  d0[0] = d0z0[0];
687  Double_t d0Mother = TMath::Abs(d0[0]);
688  Double_t d0firstTrack = TMath::Abs(candidateBPlus->Getd0Prong(0));
689  Double_t d0secondTrack = TMath::Abs(candidateBPlus->Getd0Prong(1));
690 
691  Double_t impactProduct = candidateBPlus->Prodd0d0();
692  Double_t impactProductXY = TMath::Abs(candidateBPlus->ImpParXY());
693 
694  Double_t angleBetweenBothDaughters = (candidateD0->Px() * candidatePion->Px() + candidateD0->Py() * candidatePion->Py() + candidateD0->Pz() * candidatePion->Pz()) / (candidateD0->P() * candidatePion->P());
695  Double_t angleMotherFirstDaughter = (candidateBPlus->Px() * candidatePion->Px() + candidateBPlus->Py() * candidatePion->Py() + candidateBPlus->Pz() * candidatePion->Pz()) / (candidateBPlus->P() * candidatePion->P());
696  Double_t angleMotherSecondDaughter = (candidateBPlus->Px() * candidateD0->Px() + candidateBPlus->Py() * candidateD0->Py() + candidateBPlus->Pz() * candidateD0->Pz()) / (candidateBPlus->P() * candidateD0->P());
697 
698  Double_t cosThetaStar = candidateBPlus->CosThetaStar(0, 521, 211, 421);
699  Double_t vertexDistance = vertexBPlus->DistanceToVertex(primaryVertex);
700  Double_t normDecayLength = candidateBPlus->NormalizedDecayLength();
701  Double_t pdgMassMother = TDatabasePDG::Instance()->GetParticle(521)->Mass();
702  Double_t pseudoProperDecayLength = ((vertexBPlus->GetX() - primaryVertex->GetX()) * candidateBPlus->Px() / TMath::Abs(candidateBPlus->Pt())) + ((vertexBPlus->GetY() - primaryVertex->GetY()) * candidateBPlus->Py() / TMath::Abs(candidateBPlus->Pt()));
703  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother / ptMother;
704  Double_t decayTime = vertexDistance / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother / (momentumMother * momentumMother)) + 1)));
705 
706  Double_t phi = candidateBPlus->Phi();
707  Double_t theta = candidateBPlus->Theta();
708  Double_t covMatrix[21];
709  candidateBPlus->GetCovarianceXYZPxPyPz(covMatrix);
710 
711  Double_t cp = TMath::Cos(phi);
712  Double_t sp = TMath::Sin(phi);
713  Double_t ct = TMath::Cos(theta);
714  Double_t st = TMath::Sin(theta);
715 
716  Double_t errorMomentum = covMatrix[9] * cp * cp * ct * ct // GetCovPxPx
717  + covMatrix[13] * 2.*cp * sp * ct * ct // GetCovPxPy
718  + covMatrix[18] * 2.*cp * ct * st // GetCovPxPz
719  + covMatrix[14] * sp * sp * ct * ct // GetCovPyPy
720  + covMatrix[19] * 2.*sp * ct * st // GetCovPyPz
721  + covMatrix[20] * st * st; // GetCovPzPz
722  Double_t normalizedDecayTime = candidateBPlus->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother * errorMomentum * errorMomentum / (momentumMother * momentumMother)) + 1)));
723 
724  Double_t cosPointingAngleXY = candidateBPlus->CosPointingAngleXY();
725  Double_t distanceXYToVertex = vertexBPlus->DistanceXYToVertex(primaryVertex);
726  Double_t normalizedDecayLengthXY = candidateBPlus->NormalizedDecayLengthXY();
727  Double_t chi2Vertex = vertexBPlus->GetChi2perNDF();
728 
729  //Topomatic
730  Double_t dd0pr1 = 0.;
731  Double_t dd0pr2 = 0.;
732  Double_t dd0max = 0.;
733  Double_t dd0min = 0.;
734  for (Int_t ipr = 0; ipr < 2; ipr++)
735  {
736  Double_t diffIP, errdiffIP;
737  candidateBPlus->Getd0MeasMinusExpProng(ipr, bz, diffIP, errdiffIP);
738  Double_t normdd0 = 0.;
739  if (errdiffIP > 0.) normdd0 = diffIP / errdiffIP;
740  if (ipr == 0) dd0pr1 = normdd0;
741  if (ipr == 1) dd0pr2 = normdd0;
742  }
743  if (TMath::Abs(dd0pr1) > TMath::Abs(dd0pr2)) {dd0max = dd0pr1; dd0min = dd0pr2;}
744  else {dd0max = dd0pr2; dd0min = dd0pr1;}
745 
746 
747  // We apply the cuts
748  Int_t nCutIndex = 0;
749  Double_t cutVariableValue = 0.0;
750 
751  // "inv. mass width [GeV]" --------------------------------------------
752  nCutIndex = 39;
753  cutVariableValue = invMassDifference;
754  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
755  if (!bPassedCut && !fGetCutInfo) return 0;
756  //---------------------------------------------------------------------
757 
758  // "delta mass width [GeV]" -------------------------------------------
759  nCutIndex = 40;
760  cutVariableValue = invMassDelta;
761  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
762  if (!bPassedCut && !fGetCutInfo) return 0;
763  //---------------------------------------------------------------------
764 
765  // "pointing angle [Cos(theta)]" --------------------------------------
766  nCutIndex = 41;
767  cutVariableValue = pointingAngle;
768  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
769  if (!bPassedCut && !fGetCutInfo) return 0;
770  //---------------------------------------------------------------------
771 
772  // "dca [cm]" ---------------------------------------------------------
773  nCutIndex = 42;
774  cutVariableValue = dcaMother;
775  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
776  if (!bPassedCut && !fGetCutInfo) return 0;
777  //---------------------------------------------------------------------
778 
779  // "Pt BPlus [GeV/c]" ----------------------------------------------------
780  nCutIndex = 43;
781  cutVariableValue = ptMother;
782  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
783  if (!bPassedCut && !fGetCutInfo) return 0;
784  //---------------------------------------------------------------------
785 
786  // "Pt D0 [GeV/c]" -------------------------------------------------
787  nCutIndex = 44;
788  cutVariableValue = ptD0;
789  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
790  if (!bPassedCut && !fGetCutInfo) return 0;
791  //---------------------------------------------------------------------
792 
793  // "Pt Pion [GeV/c]" --------------------------------------------------
794  nCutIndex = 45;
795  cutVariableValue = ptPion;
796  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
797  if (!bPassedCut && !fGetCutInfo) return 0;
798  //---------------------------------------------------------------------
799 
800  // "d0 BPlus [cm]" -------------------------------------------------------
801  nCutIndex = 46;
802  cutVariableValue = d0Mother;
803  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
804  if (!bPassedCut && !fGetCutInfo) return 0;
805  //---------------------------------------------------------------------
806 
807  // "d0 D0 [cm]"-----------------------------------------------------
808  nCutIndex = 47;
809  cutVariableValue = d0firstTrack;
810  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
811  if (!bPassedCut && !fGetCutInfo) return 0;
812  //---------------------------------------------------------------------
813 
814  // "d0 Pion [cm]" -----------------------------------------------------
815  nCutIndex = 48;
816  cutVariableValue = d0secondTrack;
817  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
818  if (!bPassedCut && !fGetCutInfo) return 0;
819  //---------------------------------------------------------------------
820 
821  // "d0d0 [cm^2]" ------------------------------------------------------
822  nCutIndex = 49;
823  cutVariableValue = impactProduct;
824  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
825  if (!bPassedCut && !fGetCutInfo) return 0;
826  //---------------------------------------------------------------------
827 
828  // "d0d0 XY [cm^2]" ---------------------------------------------------
829  nCutIndex = 50;
830  cutVariableValue = impactProductXY;
831  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
832  if (!bPassedCut && !fGetCutInfo) return 0;
833  //---------------------------------------------------------------------
834 
835  // "angle between both daughters" -------------------------------------
836  nCutIndex = 51;
837  cutVariableValue = angleBetweenBothDaughters;
838  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
839  if (!bPassedCut && !fGetCutInfo) return 0;
840  //---------------------------------------------------------------------
841 
842  // "angle mother with first daughter" ---------------------------------
843  nCutIndex = 52;
844  cutVariableValue = angleMotherFirstDaughter;
845  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
846  if (!bPassedCut && !fGetCutInfo) return 0;
847  //---------------------------------------------------------------------
848 
849  // "angle mother with second daughter" --------------------------------
850  nCutIndex = 53;
851  cutVariableValue = angleMotherSecondDaughter;
852  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
853  if (!bPassedCut && !fGetCutInfo) return 0;
854  //---------------------------------------------------------------------
855 
856  // "cosThetaStar" -----------------------------------------------------
857  nCutIndex = 54;
858  cutVariableValue = cosThetaStar;
859  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
860  if (!bPassedCut && !fGetCutInfo) return 0;
861  //---------------------------------------------------------------------
862 
863  // "vertexDistance" ---------------------------------------------------
864  nCutIndex = 55;
865  cutVariableValue = vertexDistance;
866  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
867  if (!bPassedCut && !fGetCutInfo) return 0;
868  //---------------------------------------------------------------------
869 
870  // "pseudoProperDecayTime" --------------------------------------------
871  nCutIndex = 56;
872  cutVariableValue = pseudoProperDecayTime;
873  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
874  if (!bPassedCut && !fGetCutInfo) return 0;
875  //---------------------------------------------------------------------
876 
877  // "DecayTime" --------------------------------------------------------
878  nCutIndex = 57;
879  cutVariableValue = decayTime;
880  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
881  if (!bPassedCut && !fGetCutInfo) return 0;
882  //---------------------------------------------------------------------
883 
884  // "normalizedDecayTime" ----------------------------------------------------
885  nCutIndex = 58;
886  cutVariableValue = normalizedDecayTime;
887  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
888  if (!bPassedCut && !fGetCutInfo) return 0;
889  //---------------------------------------------------------------------
890 
891  // "normDecayLength" --------------------------------------------------
892  nCutIndex = 59;
893  cutVariableValue = normDecayLength;
894  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
895  if (!bPassedCut && !fGetCutInfo) return 0;
896  //---------------------------------------------------------------------
897 
898  // "topomatic first daughter" -----------------------------------------
899  nCutIndex = 60;
900  cutVariableValue = dd0pr1;
901  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
902  if (!bPassedCut && !fGetCutInfo) return 0;
903  //---------------------------------------------------------------------
904 
905  // "topomatic second daughter" ----------------------------------------
906  nCutIndex = 61;
907  cutVariableValue = dd0pr2;
908  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
909  if (!bPassedCut && !fGetCutInfo) return 0;
910  //---------------------------------------------------------------------
911 
912  // "topomatic max" ----------------------------------------------------
913  nCutIndex = 62;
914  cutVariableValue = dd0max;
915  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
916  if (!bPassedCut && !fGetCutInfo) return 0;
917  //---------------------------------------------------------------------
918 
919  // "topomatic min" ----------------------------------------------------
920  nCutIndex = 63;
921  cutVariableValue = dd0min;
922  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
923  if (!bPassedCut && !fGetCutInfo) return 0;
924  //---------------------------------------------------------------------
925 
926  // "pointing angle XY" ----------------------------------------------------
927  nCutIndex = 64;
928  cutVariableValue = cosPointingAngleXY;
929  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
930  if (!bPassedCut && !fGetCutInfo) return 0;
931  //---------------------------------------------------------------------
932 
933  // "vertex distance XY" ----------------------------------------------------
934  nCutIndex = 65;
935  cutVariableValue = distanceXYToVertex;
936  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
937  if (!bPassedCut && !fGetCutInfo) return 0;
938  //---------------------------------------------------------------------
939 
940  // "normalized decay length XY" ----------------------------------------------------
941  nCutIndex = 66;
942  cutVariableValue = normalizedDecayLengthXY;
943  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
944  if (!bPassedCut && !fGetCutInfo) return 0;
945  //---------------------------------------------------------------------
946 
947  // "chi squared per NDF" ----------------------------------------------------
948  nCutIndex = 67;
949  cutVariableValue = chi2Vertex;
950  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
951  if (!bPassedCut && !fGetCutInfo) return 0;
952  //---------------------------------------------------------------------
953 
954  // select D0
955  bPassedCut = IsD0FromBPlusSelected(ptMother, candidateBPlus, selectionLevel, aod, bCutArray);
956  if (!bPassedCut && !fGetCutInfo) return 0;
957  }
958 
959  if (bPassedCut == kFALSE)
960  {
961  returnvalue = 0;
962  } else
963  {
964  for (Int_t i = 39; i < 68; ++i)
965  {
966  if (bCutArray[i] == kTRUE) {
967  returnvalue = 0;
968  break;
969  }
970  }
971  }
972 
973  fIsSelectedCuts = returnvalue;
974 
975  return returnvalue;
976 }
977 //_________________________________________________________________________________________________
978 Int_t AliRDHFCutsBPlustoD0Pi::IsD0FromBPlusSelected(Double_t ptBPlus, TObject* obj, Int_t selectionLevel, AliAODEvent* aod, Bool_t bCutArray[68]) {
979  //
980  // Apply selection on D0 candidate from BPlus candidate. We have to pass the BPlus candidate to this function to get variables w.r.t. BPlus vertex.
981  //
982 
983  if (!fCutsRD) {
984  cout << "Cut matrice not inizialized. Exit..." << endl;
985  return 0;
986  }
987 
988  AliAODRecoDecayHF2Prong* candidateBPlus = (AliAODRecoDecayHF2Prong*)obj;
989  if (!candidateBPlus) {
990  cout << "candidateBPlus null" << endl;
991  return 0;
992  }
993 
994  AliAODRecoDecayHF2Prong* candidateD0 = (AliAODRecoDecayHF2Prong*)candidateBPlus->GetDaughter(1);
995  if (!candidateD0) {
996  cout << "candidateD0 null" << endl;
997  return 0;
998  }
999 
1000  AliAODTrack *candidatePion = (AliAODTrack*)candidateD0->GetDaughter(0);
1001  if (!candidatePion) {
1002  cout << "candidatePion null 2" << endl;
1003  return 0;
1004  }
1005 
1006  AliAODTrack *candidateKaon = (AliAODTrack*)candidateD0->GetDaughter(1);
1007  if (!candidateKaon) {
1008  cout << "candidateKaon null" << endl;
1009  return 0;
1010  }
1011 
1012  AliAODVertex * vertexBPlus = candidateBPlus->GetSecondaryVtx();
1013  if (!vertexBPlus) {
1014  cout << "vertexBPlus null" << endl;
1015  return 0;
1016  }
1017 
1018  AliAODVertex * vertexD0 = candidateD0->GetSecondaryVtx();
1019  if (!vertexD0) {
1020  cout << "vertexD0 null" << endl;
1021  return 0;
1022  }
1023 
1024  AliAODVertex * primaryVertex = aod->GetPrimaryVertex();
1025  if (!primaryVertex) {
1026  cout << "primaryVertex null" << endl;
1027  return 0;
1028  }
1029 
1030  Int_t returnvalue = 1;
1031  Bool_t bPassedCut = kFALSE;
1032 
1033  //get the magnetic field
1034  Double_t bz = (Double_t)aod->GetMagneticField();
1035 
1036 
1037  // selection on candidate
1038  if (selectionLevel == AliRDHFCuts::kAll ||
1039  selectionLevel == AliRDHFCuts::kCandidate) {
1040 
1041  Int_t ptbin = PtBin(ptBPlus);
1042 
1043  // D0mass
1044  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1045 
1046  // D0 window - invariant mass
1047  Int_t chargeBPlus = candidateBPlus->Charge();
1048  UInt_t prongs[2];
1049  if (chargeBPlus == 1)
1050  {
1051  prongs[0] = 211;
1052  prongs[1] = 321;
1053  }
1054  else if (chargeBPlus == -1)
1055  {
1056  prongs[1] = 211;
1057  prongs[0] = 321;
1058  }
1059  else
1060  {
1061  cout << "Wrong charge BPlus." << endl;
1062  return 0;
1063  }
1064  Double_t invMassD0 = candidateD0->InvMass(2, prongs);
1065  Double_t invMassDifference = TMath::Abs(mD0PDG - invMassD0);
1066 
1067  Double_t pointingAngle = candidateD0->CosPointingAngle();
1068  Double_t dcaMother = candidateD0->GetDCA();
1069  Double_t ptMother = candidateD0->Pt();
1070  Double_t momentumMother = candidateD0->P();
1071  Double_t ptPion = candidatePion->Pt();
1072  Double_t ptKaon = candidateKaon->Pt();
1073 
1074  AliExternalTrackParam motherTrack;
1075  motherTrack.CopyFromVTrack(candidateD0);
1076  Double_t d0z0[2], covd0z0[3], d0[2];
1077  motherTrack.PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
1078  d0[0] = d0z0[0];
1079  Double_t d0Mother = TMath::Abs(d0[0]);
1080  Double_t d0firstTrack = TMath::Abs(candidateD0->Getd0Prong(0));
1081  Double_t d0secondTrack = TMath::Abs(candidateD0->Getd0Prong(1));
1082 
1083  Double_t impactProduct = candidateD0->Prodd0d0();
1084  Double_t impactProductXY = TMath::Abs(candidateD0->ImpParXY());
1085 
1086  Double_t angleBetweenBothDaughters = (candidateKaon->Px() * candidatePion->Px() + candidateKaon->Py() * candidatePion->Py() + candidateKaon->Pz() * candidatePion->Pz()) / (candidateKaon->P() * candidatePion->P());
1087  Double_t angleMotherFirstDaughter = (candidateD0->Px() * candidatePion->Px() + candidateD0->Py() * candidatePion->Py() + candidateD0->Pz() * candidatePion->Pz()) / (candidateD0->P() * candidatePion->P());
1088  Double_t angleMotherSecondDaughter = (candidateD0->Px() * candidateKaon->Px() + candidateD0->Py() * candidateKaon->Py() + candidateD0->Pz() * candidateKaon->Pz()) / (candidateD0->P() * candidateKaon->P());
1089 
1090  Double_t cosThetaStar = candidateD0->CosThetaStar(0, 421, 211, 321);
1091  Double_t vertexDistance = vertexD0->DistanceToVertex(primaryVertex);
1092  Double_t normDecayLength = candidateD0->NormalizedDecayLength();
1093  Double_t pdgMassMother = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1094  Double_t pseudoProperDecayLength = ((vertexD0->GetX() - primaryVertex->GetX()) * candidateD0->Px() / TMath::Abs(candidateD0->Pt())) + ((vertexD0->GetY() - primaryVertex->GetY()) * candidateD0->Py() / TMath::Abs(candidateD0->Pt()));
1095  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother / ptMother;
1096  Double_t decayTime = vertexDistance / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother / (momentumMother * momentumMother)) + 1)));
1097 
1098  Double_t phi = candidateD0->Phi();
1099  Double_t theta = candidateD0->Theta();
1100  Double_t covMatrix[21];
1101  candidateD0->GetCovarianceXYZPxPyPz(covMatrix);
1102 
1103  Double_t cp = TMath::Cos(phi);
1104  Double_t sp = TMath::Sin(phi);
1105  Double_t ct = TMath::Cos(theta);
1106  Double_t st = TMath::Sin(theta);
1107 
1108  Double_t errorMomentum = covMatrix[9] * cp * cp * ct * ct // GetCovPxPx
1109  + covMatrix[13] * 2.*cp * sp * ct * ct // GetCovPxPy
1110  + covMatrix[18] * 2.*cp * ct * st // GetCovPxPz
1111  + covMatrix[14] * sp * sp * ct * ct // GetCovPyPy
1112  + covMatrix[19] * 2.*sp * ct * st // GetCovPyPz
1113  + covMatrix[20] * st * st; // GetCovPzPz
1114  Double_t normalizedDecayTime = candidateD0->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother * errorMomentum * errorMomentum / (momentumMother * momentumMother)) + 1)));
1115 
1116  Double_t cosPointingAngleXY = candidateD0->CosPointingAngleXY();
1117  Double_t distanceXYToVertex = vertexD0->DistanceXYToVertex(primaryVertex);
1118  Double_t normalizedDecayLengthXY = candidateD0->NormalizedDecayLengthXY();
1119  Double_t chi2Vertex = vertexD0->GetChi2perNDF();
1120 
1121 
1122  //Topomatic
1123  Double_t dd0pr1 = 0.;
1124  Double_t dd0pr2 = 0.;
1125  Double_t dd0max = 0.;
1126  Double_t dd0min = 0.;
1127  for (Int_t ipr = 0; ipr < 2; ipr++)
1128  {
1129  Double_t diffIP, errdiffIP;
1130  candidateD0->Getd0MeasMinusExpProng(ipr, bz, diffIP, errdiffIP);
1131  Double_t normdd0 = 0.;
1132  if (errdiffIP > 0.) normdd0 = diffIP / errdiffIP;
1133  if (ipr == 0) dd0pr1 = normdd0;
1134  if (ipr == 1) dd0pr2 = normdd0;
1135  }
1136  if (TMath::Abs(dd0pr1) > TMath::Abs(dd0pr2)) {dd0max = dd0pr1; dd0min = dd0pr2;}
1137  else {dd0max = dd0pr2; dd0min = dd0pr1;}
1138 
1139 
1140  // We apply the cuts
1141  Int_t nCutIndex = 0;
1142  Double_t cutVariableValue = 0.0;
1143 
1144  // "inv. mass width [GeV]" --------------------------------------------
1145  nCutIndex = 0;
1146  cutVariableValue = invMassDifference;
1147  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1148  if (!bPassedCut && !fGetCutInfo) return 0;
1149  //---------------------------------------------------------------------
1150 
1151  // "delta mass width [GeV]" -------------------------------------------
1152  // nCutIndex = 1; // not used for D0
1153  // cutVariableValue = invMassDelta;
1154  // bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1155  // if(!bPassedCut && !fGetCutInfo) return 0;
1156  //---------------------------------------------------------------------
1157 
1158  // "pointing angle [Cos(theta)]" --------------------------------------
1159  nCutIndex = 2;
1160  cutVariableValue = pointingAngle;
1161  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1162  if (!bPassedCut && !fGetCutInfo) return 0;
1163  //---------------------------------------------------------------------
1164 
1165  // "dca [cm]" ---------------------------------------------------------
1166  nCutIndex = 3;
1167  cutVariableValue = dcaMother;
1168  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1169  if (!bPassedCut && !fGetCutInfo) return 0;
1170  //---------------------------------------------------------------------
1171 
1172  // "Pt D0 [GeV/c]" ----------------------------------------------------
1173  nCutIndex = 4;
1174  cutVariableValue = ptMother;
1175  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1176  if (!bPassedCut && !fGetCutInfo) return 0;
1177  //---------------------------------------------------------------------
1178 
1179  // "Pt Kaon [GeV/c]" -------------------------------------------------
1180  nCutIndex = 5;
1181  cutVariableValue = ptKaon;
1182  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1183  if (!bPassedCut && !fGetCutInfo) return 0;
1184  //---------------------------------------------------------------------
1185 
1186  // "Pt Pion [GeV/c]" --------------------------------------------------
1187  nCutIndex = 6;
1188  cutVariableValue = ptPion;
1189  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1190  if (!bPassedCut && !fGetCutInfo) return 0;
1191  //---------------------------------------------------------------------
1192 
1193  // "d0 D0 [cm]" -------------------------------------------------------
1194  nCutIndex = 7;
1195  cutVariableValue = d0Mother;
1196  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1197  if (!bPassedCut && !fGetCutInfo) return 0;
1198  //---------------------------------------------------------------------
1199 
1200  // "d0 Kaon [cm]"-----------------------------------------------------
1201  nCutIndex = 8;
1202  cutVariableValue = d0firstTrack;
1203  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1204  if (!bPassedCut && !fGetCutInfo) return 0;
1205  //---------------------------------------------------------------------
1206 
1207  // "d0 Pion [cm]" -----------------------------------------------------
1208  nCutIndex = 9;
1209  cutVariableValue = d0secondTrack;
1210  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1211  if (!bPassedCut && !fGetCutInfo) return 0;
1212  //---------------------------------------------------------------------
1213 
1214  // "d0d0 [cm^2]" ------------------------------------------------------
1215  nCutIndex = 10;
1216  cutVariableValue = impactProduct;
1217  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1218  if (!bPassedCut && !fGetCutInfo) return 0;
1219  //---------------------------------------------------------------------
1220 
1221  // "d0d0 XY [cm^2]" ---------------------------------------------------
1222  nCutIndex = 11;
1223  cutVariableValue = impactProductXY;
1224  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1225  if (!bPassedCut && !fGetCutInfo) return 0;
1226  //---------------------------------------------------------------------
1227 
1228  // "angle between both daughters" -------------------------------------
1229  nCutIndex = 12;
1230  cutVariableValue = angleBetweenBothDaughters;
1231  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1232  if (!bPassedCut && !fGetCutInfo) return 0;
1233  //---------------------------------------------------------------------
1234 
1235  // "angle mother with first daughter" ---------------------------------
1236  nCutIndex = 13;
1237  cutVariableValue = angleMotherFirstDaughter;
1238  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1239  if (!bPassedCut && !fGetCutInfo) return 0;
1240  //---------------------------------------------------------------------
1241 
1242  // "angle mother with second daughter" --------------------------------
1243  nCutIndex = 14;
1244  cutVariableValue = angleMotherSecondDaughter;
1245  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1246  if (!bPassedCut && !fGetCutInfo) return 0;
1247  //---------------------------------------------------------------------
1248 
1249  // "cosThetaStar" -----------------------------------------------------
1250  nCutIndex = 15;
1251  cutVariableValue = cosThetaStar;
1252  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1253  if (!bPassedCut && !fGetCutInfo) return 0;
1254  //---------------------------------------------------------------------
1255 
1256  // "vertexDistance" ---------------------------------------------------
1257  nCutIndex = 16;
1258  cutVariableValue = vertexDistance;
1259  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1260  if (!bPassedCut && !fGetCutInfo) return 0;
1261  //---------------------------------------------------------------------
1262 
1263  // "pseudoProperDecayTime" --------------------------------------------
1264  nCutIndex = 17;
1265  cutVariableValue = pseudoProperDecayTime;
1266  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1267  if (!bPassedCut && !fGetCutInfo) return 0;
1268  //---------------------------------------------------------------------
1269 
1270  // "DecayTime" --------------------------------------------------------
1271  nCutIndex = 18;
1272  cutVariableValue = decayTime;
1273  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1274  if (!bPassedCut && !fGetCutInfo) return 0;
1275  //---------------------------------------------------------------------
1276 
1277  // "normalizedDecayTime" ----------------------------------------------------
1278  nCutIndex = 19;
1279  cutVariableValue = normalizedDecayTime;
1280  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1281  if (!bPassedCut && !fGetCutInfo) return 0;
1282  //---------------------------------------------------------------------
1283 
1284  // "normDecayLength" --------------------------------------------------
1285  nCutIndex = 20;
1286  cutVariableValue = normDecayLength;
1287  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1288  if (!bPassedCut && !fGetCutInfo) return 0;
1289  //---------------------------------------------------------------------
1290 
1291  // "topomatic first daughter" -----------------------------------------
1292  nCutIndex = 21;
1293  cutVariableValue = dd0pr1;
1294  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1295  if (!bPassedCut && !fGetCutInfo) return 0;
1296  //---------------------------------------------------------------------
1297 
1298  // "topomatic second daughter" ----------------------------------------
1299  nCutIndex = 22;
1300  cutVariableValue = dd0pr2;
1301  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1302  if (!bPassedCut && !fGetCutInfo) return 0;
1303  //---------------------------------------------------------------------
1304 
1305  // "topomatic max" ----------------------------------------------------
1306  nCutIndex = 23;
1307  cutVariableValue = dd0max;
1308  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1309  if (!bPassedCut && !fGetCutInfo) return 0;
1310  //---------------------------------------------------------------------
1311 
1312  // "topomatic min" ----------------------------------------------------
1313  nCutIndex = 24;
1314  cutVariableValue = dd0min;
1315  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1316  if (!bPassedCut && !fGetCutInfo) return 0;
1317  //---------------------------------------------------------------------
1318 
1319  // "pointing angle XY" ----------------------------------------------------
1320  nCutIndex = 25;
1321  cutVariableValue = cosPointingAngleXY;
1322  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1323  if (!bPassedCut && !fGetCutInfo) return 0;
1324  //---------------------------------------------------------------------
1325 
1326  // "vertex distance XY" ----------------------------------------------------
1327  nCutIndex = 26;
1328  cutVariableValue = distanceXYToVertex;
1329  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1330  if (!bPassedCut && !fGetCutInfo) return 0;
1331  //---------------------------------------------------------------------
1332 
1333  // "normalized decay length XY" ----------------------------------------------------
1334  nCutIndex = 27;
1335  cutVariableValue = normalizedDecayLengthXY;
1336  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1337  if (!bPassedCut && !fGetCutInfo) return 0;
1338  //---------------------------------------------------------------------
1339 
1340  // "chi squared per NDF" ----------------------------------------------------
1341  nCutIndex = 28;
1342  cutVariableValue = chi2Vertex;
1343  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1344  if (!bPassedCut && !fGetCutInfo) return 0;
1345  //---------------------------------------------------------------------
1346 
1347 
1348 
1349  AliAODRecoDecay* candidateD0toBPlus = (AliAODRecoDecay*)candidateD0;
1350  AliExternalTrackParam pionD0Track;
1351  AliExternalTrackParam kaonD0Track;
1352 
1353  Double_t d0z0DSVert[2], covd0z0DSVert[3], d0DSVert[2];
1354 
1355  pionD0Track.CopyFromVTrack(candidatePion);
1356  pionD0Track.PropagateToDCA(vertexBPlus, bz, 100., d0z0DSVert, covd0z0DSVert);
1357  d0DSVert[0] = d0z0DSVert[0];
1358 
1359  kaonD0Track.CopyFromVTrack(candidateKaon);
1360  kaonD0Track.PropagateToDCA(vertexBPlus, bz, 100., d0z0DSVert, covd0z0DSVert);
1361  d0DSVert[1] = d0z0DSVert[0];
1362 
1363  AliExternalTrackParam D0Track;
1364  D0Track.CopyFromVTrack(candidateD0);
1365  Double_t d0z0D0DSVert[2], covd0z0D0DSVert[3], d0D0DSVert;
1366  motherTrack.PropagateToDCA(vertexBPlus, bz, 100., d0z0D0DSVert, covd0z0D0DSVert);
1367  d0D0DSVert = TMath::Abs(d0z0D0DSVert[0]);
1368 
1369  Double_t impactProductToBPlus = d0DSVert[0] * d0DSVert[1];
1370  Double_t impactProductXYToBPlus = candidateD0toBPlus->ImpParXY(vertexBPlus);
1371 
1372  Double_t pointingAngleToBPlus = candidateD0toBPlus->CosPointingAngle(vertexBPlus);
1373  Double_t d0FirstDaughterToBPlus = TMath::Abs(d0DSVert[0]);
1374  Double_t d0SecondDaughterToBPlus = TMath::Abs(d0DSVert[1]);
1375  Double_t normDecayLengthToBPlus = candidateD0toBPlus->NormalizedDecayLength(vertexBPlus);
1376 
1377  Double_t pseudoProperDecayLengthDSVert = ((vertexD0->GetX() - vertexBPlus->GetX()) * candidateD0->Px() / TMath::Abs(candidateD0->Pt())) + ((vertexD0->GetY() - vertexBPlus->GetY()) * candidateD0->Py() / TMath::Abs(candidateD0->Pt()));
1378  Double_t pseudoProperDecayTimeToBPlus = pseudoProperDecayLengthDSVert * pdgMassMother / ptMother;
1379  Double_t DecayTimeToBPlus = vertexDistance / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother / (momentumMother * momentumMother)) + 1)));
1380 
1381  Double_t phiDSVert = candidateD0->Phi();
1382  Double_t thetaDSVert = candidateD0->Theta();
1383  Double_t covMatrixDSVert[21];
1384  candidateD0->GetCovarianceXYZPxPyPz(covMatrixDSVert);
1385 
1386  cp = TMath::Cos(phiDSVert);
1387  sp = TMath::Sin(phiDSVert);
1388  ct = TMath::Cos(thetaDSVert);
1389  st = TMath::Sin(thetaDSVert);
1390 
1391  errorMomentum = covMatrix[9] * cp * cp * ct * ct // GetCovPxPx
1392  + covMatrix[13] * 2.*cp * sp * ct * ct // GetCovPxPy
1393  + covMatrix[18] * 2.*cp * ct * st // GetCovPxPz
1394  + covMatrix[14] * sp * sp * ct * ct // GetCovPyPy
1395  + covMatrix[19] * 2.*sp * ct * st // GetCovPyPz
1396  + covMatrix[20] * st * st; // GetCovPzPz
1397  Double_t normalizedDecayTimeToBPlus = candidateD0toBPlus->NormalizedDecayLength(vertexBPlus) / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother * errorMomentum * errorMomentum / (momentumMother * momentumMother)) + 1)));
1398 
1399  // "pointingAngleToBPlus" ---------------------------------------------
1400  nCutIndex = 29;
1401  cutVariableValue = pointingAngleToBPlus;
1402  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1403  if (!bPassedCut && !fGetCutInfo) return 0;
1404  //---------------------------------------------------------------------
1405 
1406  // "d0MotherToBPlus" --------------------------------------------------
1407  nCutIndex = 30;
1408  cutVariableValue = d0D0DSVert;
1409  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1410  if (!bPassedCut && !fGetCutInfo) return 0;
1411  //---------------------------------------------------------------------
1412 
1413  // "d0FirstDaughterToBPlus" -------------------------------------------
1414  nCutIndex = 31;
1415  cutVariableValue = d0FirstDaughterToBPlus;
1416  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1417  if (!bPassedCut && !fGetCutInfo) return 0;
1418  //---------------------------------------------------------------------
1419 
1420  // "d0SecondDaughterToBPlus" ------------------------------------------
1421  nCutIndex = 32;
1422  cutVariableValue = d0SecondDaughterToBPlus;
1423  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1424  if (!bPassedCut && !fGetCutInfo) return 0;
1425  //---------------------------------------------------------------------
1426 
1427  // "impactProductToBPlus" ---------------------------------------------
1428  nCutIndex = 33;
1429  cutVariableValue = impactProductToBPlus;
1430  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1431  if (!bPassedCut && !fGetCutInfo) return 0;
1432  //---------------------------------------------------------------------
1433 
1434  // "impactProductXYToBPlus" -------------------------------------------
1435  nCutIndex = 34;
1436  cutVariableValue = impactProductXYToBPlus;
1437  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1438  if (!bPassedCut && !fGetCutInfo) return 0;
1439  //---------------------------------------------------------------------
1440 
1441  // "normDecayLengthToBPlus" -------------------------------------------
1442  nCutIndex = 35;
1443  cutVariableValue = normDecayLengthToBPlus;
1444  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1445  if (!bPassedCut && !fGetCutInfo) return 0;
1446  //---------------------------------------------------------------------
1447 
1448  // "pseudoProperDecayTimeToBPlus" -------------------------------------
1449  nCutIndex = 36;
1450  cutVariableValue = pseudoProperDecayTimeToBPlus;
1451  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1452  if (!bPassedCut && !fGetCutInfo) return 0;
1453  //---------------------------------------------------------------------
1454 
1455  // "DecayTimeToBPlus" -------------------------------------------------
1456  nCutIndex = 37;
1457  cutVariableValue = DecayTimeToBPlus;
1458  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1459  if (!bPassedCut && !fGetCutInfo) return 0;
1460  //---------------------------------------------------------------------
1461 
1462  // "normalizedDecayTimeToBPlus" ---------------------------------------------
1463  nCutIndex = 38;
1464  cutVariableValue = normalizedDecayTimeToBPlus;
1465  bPassedCut = ApplyCutOnVariable(nCutIndex, ptbin, cutVariableValue, bCutArray);
1466  if (!bPassedCut && !fGetCutInfo) return 0;
1467  //---------------------------------------------------------------------
1468  }
1469 
1470  for (Int_t i = 0; i < 39; ++i)
1471  {
1472  if (bCutArray[i] == kTRUE) {
1473  returnvalue = 0;
1474  break;
1475  }
1476  }
1477 
1478  return returnvalue;
1479 }
1480 //----------------------------------------------------------------------------------
1482 //
1483  // Apply selection on D0 candidate.
1484  //
1485 
1486  if (!fCutsRDD0forD0ptbin) {
1487  cout << "Cut matrice not inizialized. Exit..." << endl;
1488  return 0;
1489  }
1490 
1491  AliAODRecoDecayHF2Prong* candidateD0 = (AliAODRecoDecayHF2Prong*)obj;
1492  if (!candidateD0) {
1493  cout << "candidateD0 null" << endl;
1494  return 0;
1495  }
1496 
1497  AliAODTrack *candidatePion = (AliAODTrack*)candidateD0->GetDaughter(0);
1498  if (!candidatePion) {
1499  cout << "candidatePion null 3" << endl;
1500  return 0;
1501  }
1502 
1503  AliAODTrack *candidateKaon = (AliAODTrack*)candidateD0->GetDaughter(1);
1504  if (!candidateKaon) {
1505  cout << "candidateKaon null" << endl;
1506  return 0;
1507  }
1508 
1509  AliAODVertex * vertexD0 = candidateD0->GetSecondaryVtx();
1510  if (!vertexD0) {
1511  cout << "vertexD0 null" << endl;
1512  return 0;
1513  }
1514 
1515  AliAODVertex * primaryVertex = aod->GetPrimaryVertex();
1516  if (!primaryVertex) {
1517  cout << "primaryVertex null" << endl;
1518  return 0;
1519  }
1520 
1521  Int_t returnvalue = 1;
1522  Bool_t bPassedCut = kFALSE;
1523 
1524  //get the magnetic field
1525  Double_t bz = (Double_t)aod->GetMagneticField();
1526 
1527 
1528  // selection on candidate
1529  if (selectionLevel == AliRDHFCuts::kAll ||
1530  selectionLevel == AliRDHFCuts::kCandidate) {
1531 
1532  Int_t ptbin = PtBinD0forD0ptbin(candidateD0->Pt());
1533  if (ptbin == -1) return -1;
1534 
1535  // D0mass
1536  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1537 
1538  // D0 window - invariant mass
1539  UInt_t prongs[2];
1540  prongs[0] = 211; prongs[1] = 321;
1541  Double_t invMassD0 = candidateD0->InvMass(2, prongs);
1542  Double_t invMassDifference = TMath::Abs(mD0PDG - invMassD0);
1543 
1544  UInt_t prongs2[2];
1545  prongs2[1] = 211; prongs2[0] = 321;
1546  Double_t invMassD02 = candidateD0->InvMass(2, prongs2);
1547  Double_t invMassDifference2 = TMath::Abs(mD0PDG - invMassD02);
1548 
1549  if (invMassDifference > invMassDifference2) invMassDifference = invMassDifference2;
1550 
1551  Double_t pointingAngle = candidateD0->CosPointingAngle();
1552  Double_t dcaMother = candidateD0->GetDCA();
1553  Double_t ptMother = candidateD0->Pt();
1554  Double_t momentumMother = candidateD0->P();
1555  Double_t ptPion = candidatePion->Pt();
1556  Double_t ptKaon = candidateKaon->Pt();
1557 
1558  AliExternalTrackParam motherTrack;
1559  motherTrack.CopyFromVTrack(candidateD0);
1560  Double_t d0z0[2], covd0z0[3], d0[2];
1561  motherTrack.PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
1562  d0[0] = d0z0[0];
1563  Double_t d0Mother = TMath::Abs(d0[0]);
1564  Double_t d0firstTrack = TMath::Abs(candidateD0->Getd0Prong(0));
1565  Double_t d0secondTrack = TMath::Abs(candidateD0->Getd0Prong(1));
1566 
1567  Double_t impactProduct = candidateD0->Prodd0d0();
1568  Double_t impactProductXY = TMath::Abs(candidateD0->ImpParXY());
1569 
1570  Double_t angleBetweenBothDaughters = (candidateKaon->Px() * candidatePion->Px() + candidateKaon->Py() * candidatePion->Py() + candidateKaon->Pz() * candidatePion->Pz()) / (candidateKaon->P() * candidatePion->P());
1571  Double_t angleMotherFirstDaughter = (candidateD0->Px() * candidatePion->Px() + candidateD0->Py() * candidatePion->Py() + candidateD0->Pz() * candidatePion->Pz()) / (candidateD0->P() * candidatePion->P());
1572  Double_t angleMotherSecondDaughter = (candidateD0->Px() * candidateKaon->Px() + candidateD0->Py() * candidateKaon->Py() + candidateD0->Pz() * candidateKaon->Pz()) / (candidateD0->P() * candidateKaon->P());
1573 
1574  Double_t cosThetaStar = candidateD0->CosThetaStar(0, 421, 211, 321);
1575  Double_t vertexDistance = vertexD0->DistanceToVertex(primaryVertex);
1576  Double_t normDecayLength = candidateD0->NormalizedDecayLength();
1577  Double_t pdgMassMother = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1578  Double_t pseudoProperDecayLength = ((vertexD0->GetX() - primaryVertex->GetX()) * candidateD0->Px() / TMath::Abs(candidateD0->Pt())) + ((vertexD0->GetY() - primaryVertex->GetY()) * candidateD0->Py() / TMath::Abs(candidateD0->Pt()));
1579  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother / ptMother;
1580  Double_t decayTime = vertexDistance / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother / (momentumMother * momentumMother)) + 1)));
1581 
1582  Double_t phi = candidateD0->Phi();
1583  Double_t theta = candidateD0->Theta();
1584  Double_t covMatrix[21];
1585  candidateD0->GetCovarianceXYZPxPyPz(covMatrix);
1586 
1587  Double_t cp = TMath::Cos(phi);
1588  Double_t sp = TMath::Sin(phi);
1589  Double_t ct = TMath::Cos(theta);
1590  Double_t st = TMath::Sin(theta);
1591 
1592  Double_t errorMomentum = covMatrix[9] * cp * cp * ct * ct // GetCovPxPx
1593  + covMatrix[13] * 2.*cp * sp * ct * ct // GetCovPxPy
1594  + covMatrix[18] * 2.*cp * ct * st // GetCovPxPz
1595  + covMatrix[14] * sp * sp * ct * ct // GetCovPyPy
1596  + covMatrix[19] * 2.*sp * ct * st // GetCovPyPz
1597  + covMatrix[20] * st * st; // GetCovPzPz
1598  Double_t normalizedDecayTime = candidateD0->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother * errorMomentum * errorMomentum / (momentumMother * momentumMother)) + 1)));
1599 
1600  Double_t cosPointingAngleXY = candidateD0->CosPointingAngleXY();
1601  Double_t distanceXYToVertex = vertexD0->DistanceXYToVertex(primaryVertex);
1602  Double_t normalizedDecayLengthXY = candidateD0->NormalizedDecayLengthXY();
1603  Double_t chi2Vertex = vertexD0->GetChi2perNDF();
1604 
1605  //Topomatic
1606  Double_t dd0pr1 = 0.;
1607  Double_t dd0pr2 = 0.;
1608  Double_t dd0max = 0.;
1609  Double_t dd0min = 0.;
1610  for (Int_t ipr = 0; ipr < 2; ipr++)
1611  {
1612  Double_t diffIP, errdiffIP;
1613  candidateD0->Getd0MeasMinusExpProng(ipr, bz, diffIP, errdiffIP);
1614  Double_t normdd0 = 0.;
1615  if (errdiffIP > 0.) normdd0 = diffIP / errdiffIP;
1616  if (ipr == 0) dd0pr1 = normdd0;
1617  if (ipr == 1) dd0pr2 = normdd0;
1618  }
1619  if (TMath::Abs(dd0pr1) > TMath::Abs(dd0pr2)) {dd0max = dd0pr1; dd0min = dd0pr2;}
1620  else {dd0max = dd0pr2; dd0min = dd0pr1;}
1621 
1622 
1623  // We apply the cuts
1624  Int_t nCutIndex = 0;
1625  Double_t cutVariableValue = 0.0;
1626 
1627  // "inv. mass width [GeV]" --------------------------------------------
1628  nCutIndex = 0;
1629  cutVariableValue = invMassDifference;
1630  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1631  if (!bPassedCut && !fGetCutInfo) return 0;
1632  //---------------------------------------------------------------------
1633 
1634  // "delta mass width [GeV]" -------------------------------------------
1635  // nCutIndex = 1; // not used for D0
1636  // cutVariableValue = invMassDelta;
1637  // bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1638  // if(!bPassedCut && !fGetCutInfo) return 0;
1639  //---------------------------------------------------------------------
1640 
1641  // "pointing angle [Cos(theta)]" --------------------------------------
1642  nCutIndex = 2;
1643  cutVariableValue = pointingAngle;
1644  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1645  if (!bPassedCut && !fGetCutInfo) return 0;
1646  //---------------------------------------------------------------------
1647 
1648  // "dca [cm]" ---------------------------------------------------------
1649  nCutIndex = 3;
1650  cutVariableValue = dcaMother;
1651  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1652  if (!bPassedCut && !fGetCutInfo) return 0;
1653  //---------------------------------------------------------------------
1654 
1655  // "Pt D0 [GeV/c]" ----------------------------------------------------
1656  nCutIndex = 4;
1657  cutVariableValue = ptMother;
1658  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1659  if (!bPassedCut && !fGetCutInfo) return 0;
1660  //---------------------------------------------------------------------
1661 
1662  // "Pt Kaon [GeV/c]" -------------------------------------------------
1663  nCutIndex = 5;
1664  cutVariableValue = ptKaon;
1665  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1666  if (!bPassedCut && !fGetCutInfo) return 0;
1667  //---------------------------------------------------------------------
1668 
1669  // "Pt Pion [GeV/c]" --------------------------------------------------
1670  nCutIndex = 6;
1671  cutVariableValue = ptPion;
1672  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1673  if (!bPassedCut && !fGetCutInfo) return 0;
1674  //---------------------------------------------------------------------
1675 
1676  // "d0 D0 [cm]" -------------------------------------------------------
1677  nCutIndex = 7;
1678  cutVariableValue = d0Mother;
1679  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1680  if (!bPassedCut && !fGetCutInfo) return 0;
1681  //---------------------------------------------------------------------
1682 
1683  // "d0 Kaon [cm]"-----------------------------------------------------
1684  nCutIndex = 8;
1685  cutVariableValue = d0firstTrack;
1686  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1687  if (!bPassedCut && !fGetCutInfo) return 0;
1688  //---------------------------------------------------------------------
1689 
1690  // "d0 Pion [cm]" -----------------------------------------------------
1691  nCutIndex = 9;
1692  cutVariableValue = d0secondTrack;
1693  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1694  if (!bPassedCut && !fGetCutInfo) return 0;
1695  //---------------------------------------------------------------------
1696 
1697  // "d0d0 [cm^2]" ------------------------------------------------------
1698  nCutIndex = 10;
1699  cutVariableValue = impactProduct;
1700  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1701  if (!bPassedCut && !fGetCutInfo) return 0;
1702  //---------------------------------------------------------------------
1703 
1704  // "d0d0 XY [cm^2]" ---------------------------------------------------
1705  nCutIndex = 11;
1706  cutVariableValue = impactProductXY;
1707  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1708  if (!bPassedCut && !fGetCutInfo) return 0;
1709  //---------------------------------------------------------------------
1710 
1711  // "angle between both daughters" -------------------------------------
1712  nCutIndex = 12;
1713  cutVariableValue = angleBetweenBothDaughters;
1714  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1715  if (!bPassedCut && !fGetCutInfo) return 0;
1716  //---------------------------------------------------------------------
1717 
1718  // "angle mother with first daughter" ---------------------------------
1719  nCutIndex = 13;
1720  cutVariableValue = angleMotherFirstDaughter;
1721  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1722  if (!bPassedCut && !fGetCutInfo) return 0;
1723  //---------------------------------------------------------------------
1724 
1725  // "angle mother with second daughter" --------------------------------
1726  nCutIndex = 14;
1727  cutVariableValue = angleMotherSecondDaughter;
1728  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1729  if (!bPassedCut && !fGetCutInfo) return 0;
1730  //---------------------------------------------------------------------
1731 
1732  // "cosThetaStar" -----------------------------------------------------
1733  nCutIndex = 15;
1734  cutVariableValue = cosThetaStar;
1735  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1736  if (!bPassedCut && !fGetCutInfo) return 0;
1737  //---------------------------------------------------------------------
1738 
1739  // "vertexDistance" ---------------------------------------------------
1740  nCutIndex = 16;
1741  cutVariableValue = vertexDistance;
1742  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1743  if (!bPassedCut && !fGetCutInfo) return 0;
1744  //---------------------------------------------------------------------
1745 
1746  // "pseudoProperDecayTime" --------------------------------------------
1747  nCutIndex = 17;
1748  cutVariableValue = pseudoProperDecayTime;
1749  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1750  if (!bPassedCut && !fGetCutInfo) return 0;
1751  //---------------------------------------------------------------------
1752 
1753  // "DecayTime" --------------------------------------------------------
1754  nCutIndex = 18;
1755  cutVariableValue = decayTime;
1756  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1757  if (!bPassedCut && !fGetCutInfo) return 0;
1758  //---------------------------------------------------------------------
1759 
1760  // "normalizedDecayTime" ----------------------------------------------------
1761  nCutIndex = 19;
1762  cutVariableValue = normalizedDecayTime;
1763  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1764  if (!bPassedCut && !fGetCutInfo) return 0;
1765  //---------------------------------------------------------------------
1766 
1767  // "normDecayLength" --------------------------------------------------
1768  nCutIndex = 20;
1769  cutVariableValue = normDecayLength;
1770  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1771  if (!bPassedCut && !fGetCutInfo) return 0;
1772  //---------------------------------------------------------------------
1773 
1774  // "topomatic first daughter" -----------------------------------------
1775  nCutIndex = 21;
1776  cutVariableValue = dd0pr1;
1777  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1778  if (!bPassedCut && !fGetCutInfo) return 0;
1779  //---------------------------------------------------------------------
1780 
1781  // "topomatic second daughter" ----------------------------------------
1782  nCutIndex = 22;
1783  cutVariableValue = dd0pr2;
1784  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1785  if (!bPassedCut && !fGetCutInfo) return 0;
1786  //---------------------------------------------------------------------
1787 
1788  // "topomatic max" ----------------------------------------------------
1789  nCutIndex = 23;
1790  cutVariableValue = dd0max;
1791  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1792  if (!bPassedCut && !fGetCutInfo) return 0;
1793  //---------------------------------------------------------------------
1794 
1795  // "topomatic min" ----------------------------------------------------
1796  nCutIndex = 24;
1797  cutVariableValue = dd0min;
1798  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1799  if (!bPassedCut && !fGetCutInfo) return 0;
1800  //---------------------------------------------------------------------
1801 
1802  // "pointing angle XY" ----------------------------------------------------
1803  nCutIndex = 25;
1804  cutVariableValue = cosPointingAngleXY;
1805  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1806  if (!bPassedCut && !fGetCutInfo) return 0;
1807  //---------------------------------------------------------------------
1808 
1809  // "vertex distance XY" ----------------------------------------------------
1810  nCutIndex = 26;
1811  cutVariableValue = distanceXYToVertex;
1812  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1813  if (!bPassedCut && !fGetCutInfo) return 0;
1814  //---------------------------------------------------------------------
1815 
1816  // "normalized decay length XY" ----------------------------------------------------
1817  nCutIndex = 27;
1818  cutVariableValue = normalizedDecayLengthXY;
1819  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1820  if (!bPassedCut && !fGetCutInfo) return 0;
1821  //---------------------------------------------------------------------
1822 
1823  // "chi squared per NDF" ----------------------------------------------------
1824  nCutIndex = 28;
1825  cutVariableValue = chi2Vertex;
1826  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex, ptbin, cutVariableValue, bCutArray);
1827  if (!bPassedCut && !fGetCutInfo) return 0;
1828  //---------------------------------------------------------------------
1829 
1830  }
1831 
1832  for (Int_t i = 0; i < 29; ++i)
1833  {
1834  if (bCutArray[i] == kTRUE) {
1835  returnvalue = 0;
1836  break;
1837  }
1838  }
1839 
1840  return returnvalue;
1841 }
1842 //--------------------------------------------------------------------------
1844  //
1845  // In this function we apply the selection cuts on the BPlus candidate and its daughters.
1846  // The function returns 0 if the candidate is cut.
1847  //
1848 
1849  // The cuts used in this class have to be set in the maketfile.
1850  if (!fCutsRD) {
1851  cout << "Cut matrice not inizialized. Exit..." << endl;
1852  return 0;
1853  }
1854 
1855  AliAODRecoDecayHF2Prong* candidateBPlus = (AliAODRecoDecayHF2Prong*)obj;
1856  if (!candidateBPlus) {
1857  cout << "candidateBPlus null" << endl;
1858  return 0;
1859  }
1860 
1861  AliAODRecoDecayHF2Prong* candidateD0 = (AliAODRecoDecayHF2Prong*)candidateBPlus->GetDaughter(1);
1862  if (!candidateD0) {
1863  cout << "candidateD0 null" << endl;
1864  return 0;
1865  }
1866 
1867  AliAODTrack *candidatePion = (AliAODTrack*)candidateBPlus->GetDaughter(0);
1868  if (!candidatePion) {
1869  cout << "candidatePion null 1" << endl;
1870  return 0;
1871  }
1872 
1873  AliAODVertex * vertexBPlus = candidateBPlus->GetSecondaryVtx();
1874  if (!vertexBPlus) {
1875  cout << "vertexBPlus null" << endl;
1876  return 0;
1877  }
1878 
1879  AliAODVertex * primaryVertex = aod->GetPrimaryVertex();
1880  if (!primaryVertex) {
1881  cout << "primaryVertex null" << endl;
1882  return 0;
1883  }
1884 
1885  //get the magnetic field
1886  Double_t bz = (Double_t)aod->GetMagneticField();
1887 
1888  // selection on candidate
1889  if (selectionLevel == AliRDHFCuts::kAll ||
1890  selectionLevel == AliRDHFCuts::kCandidate) {
1891 
1892  // We check to which pt bin the candidate belongs
1893  Int_t ptbin = PtBin(candidateBPlus->Pt());
1894  if (ptbin == -1) return -1;
1895 
1896  // We obtain the variable values in the section below
1897  // D0Mass and BPlusmass
1898  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1899  Double_t mBPlusPDG = TDatabasePDG::Instance()->GetParticle(521)->Mass();
1900 
1901  // delta mass PDG
1902  Double_t deltaPDG = mBPlusPDG - mD0PDG;
1903 
1904  // Half width BPlus mass
1905  UInt_t prongs[2];
1906  prongs[0] = 211; prongs[1] = 421;
1907  Double_t invMassBPlus = candidateBPlus->InvMass(2, prongs);
1908  Double_t invMassDifference = TMath::Abs(mBPlusPDG - invMassBPlus);
1909  Double_t invMassDelta = TMath::Abs(deltaPDG - (DeltaInvMassBPlusKpipi(candidateBPlus)));
1910 
1911  Double_t pointingAngle = candidateBPlus->CosPointingAngle();
1912  Double_t dcaMother = candidateBPlus->GetDCA();
1913  Double_t ptMother = candidateBPlus->Pt();
1914  Double_t momentumMother = candidateBPlus->P();
1915  Double_t ptD0 = candidateD0->Pt();
1916  Double_t ptPion = candidatePion->Pt();
1917 
1918  AliExternalTrackParam motherTrack;
1919  motherTrack.CopyFromVTrack(candidateBPlus);
1920  Double_t d0z0[2], covd0z0[3], d0[2];
1921  motherTrack.PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
1922  d0[0] = d0z0[0];
1923  Double_t d0Mother = TMath::Abs(d0[0]);
1924  Double_t d0firstTrack = TMath::Abs(candidateBPlus->Getd0Prong(0));
1925  Double_t d0secondTrack = TMath::Abs(candidateBPlus->Getd0Prong(1));
1926 
1927  Double_t impactProduct = candidateBPlus->Prodd0d0();
1928  Double_t impactProductXY = TMath::Abs(candidateBPlus->ImpParXY());
1929 
1930  Double_t angleBetweenBothDaughters = (candidateD0->Px() * candidatePion->Px() + candidateD0->Py() * candidatePion->Py() + candidateD0->Pz() * candidatePion->Pz()) / (candidateD0->P() * candidatePion->P());
1931  Double_t angleMotherFirstDaughter = (candidateBPlus->Px() * candidatePion->Px() + candidateBPlus->Py() * candidatePion->Py() + candidateBPlus->Pz() * candidatePion->Pz()) / (candidateBPlus->P() * candidatePion->P());
1932  Double_t angleMotherSecondDaughter = (candidateBPlus->Px() * candidateD0->Px() + candidateBPlus->Py() * candidateD0->Py() + candidateBPlus->Pz() * candidateD0->Pz()) / (candidateBPlus->P() * candidateD0->P());
1933 
1934  Double_t cosThetaStar = candidateBPlus->CosThetaStar(0, 521, 211, 421);
1935  Double_t vertexDistance = vertexBPlus->DistanceToVertex(primaryVertex);
1936  Double_t normDecayLength = candidateBPlus->NormalizedDecayLength();
1937  Double_t pdgMassMother = TDatabasePDG::Instance()->GetParticle(521)->Mass();
1938  Double_t pseudoProperDecayLength = ((vertexBPlus->GetX() - primaryVertex->GetX()) * candidateBPlus->Px() / TMath::Abs(candidateBPlus->Pt())) + ((vertexBPlus->GetY() - primaryVertex->GetY()) * candidateBPlus->Py() / TMath::Abs(candidateBPlus->Pt()));
1939  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother / ptMother;
1940  Double_t decayTime = vertexDistance / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother / (momentumMother * momentumMother)) + 1)));
1941 
1942  Double_t phi = candidateBPlus->Phi();
1943  Double_t theta = candidateBPlus->Theta();
1944  Double_t covMatrix[21];
1945  candidateBPlus->GetCovarianceXYZPxPyPz(covMatrix);
1946 
1947  Double_t cp = TMath::Cos(phi);
1948  Double_t sp = TMath::Sin(phi);
1949  Double_t ct = TMath::Cos(theta);
1950  Double_t st = TMath::Sin(theta);
1951 
1952  Double_t errorMomentum = covMatrix[9] * cp * cp * ct * ct // GetCovPxPx
1953  + covMatrix[13] * 2.*cp * sp * ct * ct // GetCovPxPy
1954  + covMatrix[18] * 2.*cp * ct * st // GetCovPxPz
1955  + covMatrix[14] * sp * sp * ct * ct // GetCovPyPy
1956  + covMatrix[19] * 2.*sp * ct * st // GetCovPyPz
1957  + covMatrix[20] * st * st; // GetCovPzPz
1958  Double_t normalizedDecayTime = candidateBPlus->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother * errorMomentum * errorMomentum / (momentumMother * momentumMother)) + 1)));
1959 
1960  Double_t cosPointingAngleXY = candidateBPlus->CosPointingAngleXY();
1961  Double_t distanceXYToVertex = vertexBPlus->DistanceXYToVertex(primaryVertex);
1962  Double_t normalizedDecayLengthXY = candidateBPlus->NormalizedDecayLengthXY();
1963  Double_t chi2Vertex = vertexBPlus->GetChi2perNDF();
1964 
1965  //Topomatic
1966  Double_t dd0pr1 = 0.;
1967  Double_t dd0pr2 = 0.;
1968  Double_t dd0max = 0.;
1969  Double_t dd0min = 0.;
1970  for (Int_t ipr = 0; ipr < 2; ipr++)
1971  {
1972  Double_t diffIP, errdiffIP;
1973  candidateBPlus->Getd0MeasMinusExpProng(ipr, bz, diffIP, errdiffIP);
1974  Double_t normdd0 = 0.;
1975  if (errdiffIP > 0.) normdd0 = diffIP / errdiffIP;
1976  if (ipr == 0) dd0pr1 = normdd0;
1977  if (ipr == 1) dd0pr2 = normdd0;
1978  }
1979  if (TMath::Abs(dd0pr1) > TMath::Abs(dd0pr2)) {dd0max = dd0pr1; dd0min = dd0pr2;}
1980  else {dd0max = dd0pr2; dd0min = dd0pr1;}
1981 
1982 
1983  // We apply the cuts
1984  Int_t nCutIndex = 0;
1985  Double_t cutVariableValue = 0.0;
1986 
1987  // "inv. mass width [GeV]" --------------------------------------------
1988  nCutIndex = 39;
1989  cutVariableValue = invMassDifference;
1990  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
1991  //---------------------------------------------------------------------
1992 
1993  // "delta mass width [GeV]" -------------------------------------------
1994  nCutIndex = 40;
1995  cutVariableValue = invMassDelta;
1996  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
1997  //---------------------------------------------------------------------
1998 
1999  // "pointing angle [Cos(theta)]" --------------------------------------
2000  nCutIndex = 41;
2001  cutVariableValue = pointingAngle;
2002  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2003  //---------------------------------------------------------------------
2004 
2005  // "dca [cm]" ---------------------------------------------------------
2006  nCutIndex = 42;
2007  cutVariableValue = dcaMother;
2008  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2009  //---------------------------------------------------------------------
2010 
2011  // "Pt BPlus [GeV/c]" ----------------------------------------------------
2012  nCutIndex = 43;
2013  cutVariableValue = ptMother;
2014  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2015  //---------------------------------------------------------------------
2016 
2017  // "Pt D0 [GeV/c]" -------------------------------------------------
2018  nCutIndex = 44;
2019  cutVariableValue = ptD0;
2020  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2021  //---------------------------------------------------------------------
2022 
2023  // "Pt Pion [GeV/c]" --------------------------------------------------
2024  nCutIndex = 45;
2025  cutVariableValue = ptPion;
2026  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2027  //---------------------------------------------------------------------
2028 
2029  // "d0 BPlus [cm]" -------------------------------------------------------
2030  nCutIndex = 46;
2031  cutVariableValue = d0Mother;
2032  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2033  //---------------------------------------------------------------------
2034 
2035  // "d0 D0 [cm]"-----------------------------------------------------
2036  nCutIndex = 47;
2037  cutVariableValue = d0firstTrack;
2038  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2039  //---------------------------------------------------------------------
2040 
2041  // "d0 Pion [cm]" -----------------------------------------------------
2042  nCutIndex = 48;
2043  cutVariableValue = d0secondTrack;
2044  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2045  //---------------------------------------------------------------------
2046 
2047  // "d0d0 [cm^2]" ------------------------------------------------------
2048  nCutIndex = 49;
2049  cutVariableValue = impactProduct;
2050  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2051  //---------------------------------------------------------------------
2052 
2053  // "d0d0 XY [cm^2]" ---------------------------------------------------
2054  nCutIndex = 50;
2055  cutVariableValue = impactProductXY;
2056  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2057  //---------------------------------------------------------------------
2058 
2059  // "angle between both daughters" -------------------------------------
2060  nCutIndex = 51;
2061  cutVariableValue = angleBetweenBothDaughters;
2062  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2063  //---------------------------------------------------------------------
2064 
2065  // "angle mother with first daughter" ---------------------------------
2066  nCutIndex = 52;
2067  cutVariableValue = angleMotherFirstDaughter;
2068  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2069  //---------------------------------------------------------------------
2070 
2071  // "angle mother with second daughter" --------------------------------
2072  nCutIndex = 53;
2073  cutVariableValue = angleMotherSecondDaughter;
2074  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2075  //---------------------------------------------------------------------
2076 
2077  // "cosThetaStar" -----------------------------------------------------
2078  nCutIndex = 54;
2079  cutVariableValue = cosThetaStar;
2080  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2081  //---------------------------------------------------------------------
2082 
2083  // "vertexDistance" ---------------------------------------------------
2084  nCutIndex = 55;
2085  cutVariableValue = vertexDistance;
2086  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2087  //---------------------------------------------------------------------
2088 
2089  // "pseudoProperDecayTime" --------------------------------------------
2090  nCutIndex = 56;
2091  cutVariableValue = pseudoProperDecayTime;
2092  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2093  //---------------------------------------------------------------------
2094 
2095  // "DecayTime" --------------------------------------------------------
2096  nCutIndex = 57;
2097  cutVariableValue = decayTime;
2098  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2099  //---------------------------------------------------------------------
2100 
2101  // "normalizedDecayTime" ----------------------------------------------------
2102  nCutIndex = 58;
2103  cutVariableValue = normalizedDecayTime;
2104  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2105  //---------------------------------------------------------------------
2106 
2107  // "normDecayLength" --------------------------------------------------
2108  nCutIndex = 59;
2109  cutVariableValue = normDecayLength;
2110  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2111  //---------------------------------------------------------------------
2112 
2113  // "topomatic first daughter" -----------------------------------------
2114  nCutIndex = 60;
2115  cutVariableValue = dd0pr1;
2116  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2117  //---------------------------------------------------------------------
2118 
2119  // "topomatic second daughter" ----------------------------------------
2120  nCutIndex = 61;
2121  cutVariableValue = dd0pr2;
2122  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2123  //---------------------------------------------------------------------
2124 
2125  // "topomatic max" ----------------------------------------------------
2126  nCutIndex = 62;
2127  cutVariableValue = dd0max;
2128  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2129  //---------------------------------------------------------------------
2130 
2131  // "topomatic min" ----------------------------------------------------
2132  nCutIndex = 63;
2133  cutVariableValue = dd0min;
2134  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2135  //---------------------------------------------------------------------
2136 
2137  // "pointing angle XY" ----------------------------------------------------
2138  nCutIndex = 64;
2139  cutVariableValue = cosPointingAngleXY;
2140  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2141  //---------------------------------------------------------------------
2142 
2143  // "vertex distance XY" ----------------------------------------------------
2144  nCutIndex = 65;
2145  cutVariableValue = distanceXYToVertex;
2146  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2147  //---------------------------------------------------------------------
2148 
2149  // "normalized decay length XY" ----------------------------------------------------
2150  nCutIndex = 66;
2151  cutVariableValue = normalizedDecayLengthXY;
2152  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2153  //---------------------------------------------------------------------
2154 
2155  // "chi squared per NDF" ----------------------------------------------------
2156  nCutIndex = 67;
2157  cutVariableValue = chi2Vertex;
2158  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2159  //---------------------------------------------------------------------
2160 
2161  // select D0
2162  if (!IsD0FromBPlusSelectedMVA(ptMother, candidateBPlus, selectionLevel, aod, primaryVertex, bz)) return 0;
2163  }
2164 
2165  return 1;
2166 }
2167 //_________________________________________________________________________________________________
2168 Int_t AliRDHFCutsBPlustoD0Pi::IsD0FromBPlusSelectedMVA(Double_t ptBPlus, TObject* obj, Int_t selectionLevel, AliAODEvent* /*aod*/, AliAODVertex *primaryVertex, Double_t bz) {
2169  //
2170  // Apply selection on D0 candidate from BPlus candidate. We have to pass the BPlus candidate to this function to get variables w.r.t. BPlus vertex.
2171  //
2172 
2173  if (!fCutsRD) {
2174  cout << "Cut matrice not inizialized. Exit..." << endl;
2175  return 0;
2176  }
2177 
2178  AliAODRecoDecayHF2Prong* candidateBPlus = (AliAODRecoDecayHF2Prong*)obj;
2179  if (!candidateBPlus) {
2180  cout << "candidateBPlus null" << endl;
2181  return 0;
2182  }
2183 
2184  AliAODRecoDecayHF2Prong* candidateD0 = (AliAODRecoDecayHF2Prong*)candidateBPlus->GetDaughter(1);
2185  if (!candidateD0) {
2186  cout << "candidateD0 null" << endl;
2187  return 0;
2188  }
2189 
2190  AliAODTrack *candidatePion = (AliAODTrack*)candidateD0->GetDaughter(0);
2191  if (!candidatePion) {
2192  cout << "candidatePion null 2" << endl;
2193  return 0;
2194  }
2195 
2196  AliAODTrack *candidateKaon = (AliAODTrack*)candidateD0->GetDaughter(1);
2197  if (!candidateKaon) {
2198  cout << "candidateKaon null" << endl;
2199  return 0;
2200  }
2201 
2202  AliAODVertex * vertexBPlus = candidateBPlus->GetSecondaryVtx();
2203  if (!vertexBPlus) {
2204  cout << "vertexBPlus null" << endl;
2205  return 0;
2206  }
2207 
2208  AliAODVertex * vertexD0 = candidateD0->GetSecondaryVtx();
2209  if (!vertexD0) {
2210  cout << "vertexD0 null" << endl;
2211  return 0;
2212  }
2213 
2214  if (!primaryVertex) {
2215  cout << "primaryVertex null" << endl;
2216  return 0;
2217  }
2218 
2219  // selection on candidate
2220  if (selectionLevel == AliRDHFCuts::kAll ||
2221  selectionLevel == AliRDHFCuts::kCandidate) {
2222 
2223  Int_t ptbin = PtBin(ptBPlus);
2224 
2225  // D0mass
2226  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2227 
2228  // D0 window - invariant mass
2229  Int_t chargeBPlus = candidateBPlus->Charge();
2230  UInt_t prongs[2];
2231  if (chargeBPlus == 1)
2232  {
2233  prongs[0] = 211;
2234  prongs[1] = 321;
2235  }
2236  else if (chargeBPlus == -1)
2237  {
2238  prongs[1] = 211;
2239  prongs[0] = 321;
2240  }
2241  else
2242  {
2243  cout << "Wrong charge BPlus." << endl;
2244  return 0;
2245  }
2246  Double_t invMassD0 = candidateD0->InvMass(2, prongs);
2247  Double_t invMassDifference = TMath::Abs(mD0PDG - invMassD0);
2248 
2249  Double_t pointingAngle = candidateD0->CosPointingAngle();
2250  Double_t dcaMother = candidateD0->GetDCA();
2251  Double_t ptMother = candidateD0->Pt();
2252  Double_t momentumMother = candidateD0->P();
2253  Double_t ptPion = candidatePion->Pt();
2254  Double_t ptKaon = candidateKaon->Pt();
2255 
2256  AliExternalTrackParam motherTrack;
2257  motherTrack.CopyFromVTrack(candidateD0);
2258  Double_t d0z0[2], covd0z0[3], d0[2];
2259  motherTrack.PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
2260  d0[0] = d0z0[0];
2261  Double_t d0Mother = TMath::Abs(d0[0]);
2262  Double_t d0firstTrack = TMath::Abs(candidateD0->Getd0Prong(0));
2263  Double_t d0secondTrack = TMath::Abs(candidateD0->Getd0Prong(1));
2264 
2265  Double_t impactProduct = candidateD0->Prodd0d0();
2266  Double_t impactProductXY = TMath::Abs(candidateD0->ImpParXY());
2267 
2268  Double_t angleBetweenBothDaughters = (candidateKaon->Px() * candidatePion->Px() + candidateKaon->Py() * candidatePion->Py() + candidateKaon->Pz() * candidatePion->Pz()) / (candidateKaon->P() * candidatePion->P());
2269  Double_t angleMotherFirstDaughter = (candidateD0->Px() * candidatePion->Px() + candidateD0->Py() * candidatePion->Py() + candidateD0->Pz() * candidatePion->Pz()) / (candidateD0->P() * candidatePion->P());
2270  Double_t angleMotherSecondDaughter = (candidateD0->Px() * candidateKaon->Px() + candidateD0->Py() * candidateKaon->Py() + candidateD0->Pz() * candidateKaon->Pz()) / (candidateD0->P() * candidateKaon->P());
2271 
2272  Double_t cosThetaStar = candidateD0->CosThetaStar(0, 421, 211, 321);
2273  Double_t vertexDistance = vertexD0->DistanceToVertex(primaryVertex);
2274  Double_t normDecayLength = candidateD0->NormalizedDecayLength();
2275  Double_t pdgMassMother = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2276  Double_t pseudoProperDecayLength = ((vertexD0->GetX() - primaryVertex->GetX()) * candidateD0->Px() / TMath::Abs(candidateD0->Pt())) + ((vertexD0->GetY() - primaryVertex->GetY()) * candidateD0->Py() / TMath::Abs(candidateD0->Pt()));
2277  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother / ptMother;
2278  Double_t decayTime = vertexDistance / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother / (momentumMother * momentumMother)) + 1)));
2279 
2280  Double_t phi = candidateD0->Phi();
2281  Double_t theta = candidateD0->Theta();
2282  Double_t covMatrix[21];
2283  candidateD0->GetCovarianceXYZPxPyPz(covMatrix);
2284 
2285  Double_t cp = TMath::Cos(phi);
2286  Double_t sp = TMath::Sin(phi);
2287  Double_t ct = TMath::Cos(theta);
2288  Double_t st = TMath::Sin(theta);
2289 
2290  Double_t errorMomentum = covMatrix[9] * cp * cp * ct * ct // GetCovPxPx
2291  + covMatrix[13] * 2.*cp * sp * ct * ct // GetCovPxPy
2292  + covMatrix[18] * 2.*cp * ct * st // GetCovPxPz
2293  + covMatrix[14] * sp * sp * ct * ct // GetCovPyPy
2294  + covMatrix[19] * 2.*sp * ct * st // GetCovPyPz
2295  + covMatrix[20] * st * st; // GetCovPzPz
2296  Double_t normalizedDecayTime = candidateD0->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother * errorMomentum * errorMomentum / (momentumMother * momentumMother)) + 1)));
2297 
2298  Double_t cosPointingAngleXY = candidateD0->CosPointingAngleXY();
2299  Double_t distanceXYToVertex = vertexD0->DistanceXYToVertex(primaryVertex);
2300  Double_t normalizedDecayLengthXY = candidateD0->NormalizedDecayLengthXY();
2301  Double_t chi2Vertex = vertexD0->GetChi2perNDF();
2302 
2303 
2304  //Topomatic
2305  Double_t dd0pr1 = 0.;
2306  Double_t dd0pr2 = 0.;
2307  Double_t dd0max = 0.;
2308  Double_t dd0min = 0.;
2309  for (Int_t ipr = 0; ipr < 2; ipr++)
2310  {
2311  Double_t diffIP, errdiffIP;
2312  candidateD0->Getd0MeasMinusExpProng(ipr, bz, diffIP, errdiffIP);
2313  Double_t normdd0 = 0.;
2314  if (errdiffIP > 0.) normdd0 = diffIP / errdiffIP;
2315  if (ipr == 0) dd0pr1 = normdd0;
2316  if (ipr == 1) dd0pr2 = normdd0;
2317  }
2318  if (TMath::Abs(dd0pr1) > TMath::Abs(dd0pr2)) {dd0max = dd0pr1; dd0min = dd0pr2;}
2319  else {dd0max = dd0pr2; dd0min = dd0pr1;}
2320 
2321 
2322  // We apply the cuts
2323  Int_t nCutIndex = 0;
2324  Double_t cutVariableValue = 0.0;
2325 
2326  // "inv. mass width [GeV]" --------------------------------------------
2327  nCutIndex = 0;
2328  cutVariableValue = invMassDifference;
2329  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2330  //---------------------------------------------------------------------
2331 
2332  // "delta mass width [GeV]" -------------------------------------------
2333  // nCutIndex = 1; // not used for D0
2334  // cutVariableValue = invMassDelta;
2335  // if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2336  //---------------------------------------------------------------------
2337 
2338  // "pointing angle [Cos(theta)]" --------------------------------------
2339  nCutIndex = 2;
2340  cutVariableValue = pointingAngle;
2341  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2342  //---------------------------------------------------------------------
2343 
2344  // "dca [cm]" ---------------------------------------------------------
2345  nCutIndex = 3;
2346  cutVariableValue = dcaMother;
2347  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2348  //---------------------------------------------------------------------
2349 
2350  // "Pt D0 [GeV/c]" ----------------------------------------------------
2351  nCutIndex = 4;
2352  cutVariableValue = ptMother;
2353  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2354  //---------------------------------------------------------------------
2355 
2356  // "Pt Kaon [GeV/c]" -------------------------------------------------
2357  nCutIndex = 5;
2358  cutVariableValue = ptKaon;
2359  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2360  //---------------------------------------------------------------------
2361 
2362  // "Pt Pion [GeV/c]" --------------------------------------------------
2363  nCutIndex = 6;
2364  cutVariableValue = ptPion;
2365  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2366  //---------------------------------------------------------------------
2367 
2368  // "d0 D0 [cm]" -------------------------------------------------------
2369  nCutIndex = 7;
2370  cutVariableValue = d0Mother;
2371  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2372  //---------------------------------------------------------------------
2373 
2374  // "d0 Kaon [cm]"-----------------------------------------------------
2375  nCutIndex = 8;
2376  cutVariableValue = d0firstTrack;
2377  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2378  //---------------------------------------------------------------------
2379 
2380  // "d0 Pion [cm]" -----------------------------------------------------
2381  nCutIndex = 9;
2382  cutVariableValue = d0secondTrack;
2383  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2384  //---------------------------------------------------------------------
2385 
2386  // "d0d0 [cm^2]" ------------------------------------------------------
2387  nCutIndex = 10;
2388  cutVariableValue = impactProduct;
2389  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2390  //---------------------------------------------------------------------
2391 
2392  // "d0d0 XY [cm^2]" ---------------------------------------------------
2393  nCutIndex = 11;
2394  cutVariableValue = impactProductXY;
2395  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2396  //---------------------------------------------------------------------
2397 
2398  // "angle between both daughters" -------------------------------------
2399  nCutIndex = 12;
2400  cutVariableValue = angleBetweenBothDaughters;
2401  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2402  //---------------------------------------------------------------------
2403 
2404  // "angle mother with first daughter" ---------------------------------
2405  nCutIndex = 13;
2406  cutVariableValue = angleMotherFirstDaughter;
2407  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2408  //---------------------------------------------------------------------
2409 
2410  // "angle mother with second daughter" --------------------------------
2411  nCutIndex = 14;
2412  cutVariableValue = angleMotherSecondDaughter;
2413  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2414  //---------------------------------------------------------------------
2415 
2416  // "cosThetaStar" -----------------------------------------------------
2417  nCutIndex = 15;
2418  cutVariableValue = cosThetaStar;
2419  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2420  //---------------------------------------------------------------------
2421 
2422  // "vertexDistance" ---------------------------------------------------
2423  nCutIndex = 16;
2424  cutVariableValue = vertexDistance;
2425  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2426  //---------------------------------------------------------------------
2427 
2428  // "pseudoProperDecayTime" --------------------------------------------
2429  nCutIndex = 17;
2430  cutVariableValue = pseudoProperDecayTime;
2431  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2432  //---------------------------------------------------------------------
2433 
2434  // "DecayTime" --------------------------------------------------------
2435  nCutIndex = 18;
2436  cutVariableValue = decayTime;
2437  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2438  //---------------------------------------------------------------------
2439 
2440  // "normalizedDecayTime" ----------------------------------------------------
2441  nCutIndex = 19;
2442  cutVariableValue = normalizedDecayTime;
2443  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2444  //---------------------------------------------------------------------
2445 
2446  // "normDecayLength" --------------------------------------------------
2447  nCutIndex = 20;
2448  cutVariableValue = normDecayLength;
2449  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2450  //---------------------------------------------------------------------
2451 
2452  // "topomatic first daughter" -----------------------------------------
2453  nCutIndex = 21;
2454  cutVariableValue = dd0pr1;
2455  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2456  //---------------------------------------------------------------------
2457 
2458  // "topomatic second daughter" ----------------------------------------
2459  nCutIndex = 22;
2460  cutVariableValue = dd0pr2;
2461  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2462  //---------------------------------------------------------------------
2463 
2464  // "topomatic max" ----------------------------------------------------
2465  nCutIndex = 23;
2466  cutVariableValue = dd0max;
2467  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2468  //---------------------------------------------------------------------
2469 
2470  // "topomatic min" ----------------------------------------------------
2471  nCutIndex = 24;
2472  cutVariableValue = dd0min;
2473  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2474  //---------------------------------------------------------------------
2475 
2476  // "pointing angle XY" ----------------------------------------------------
2477  nCutIndex = 25;
2478  cutVariableValue = cosPointingAngleXY;
2479  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2480  //---------------------------------------------------------------------
2481 
2482  // "vertex distance XY" ----------------------------------------------------
2483  nCutIndex = 26;
2484  cutVariableValue = distanceXYToVertex;
2485  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2486  //---------------------------------------------------------------------
2487 
2488  // "normalized decay length XY" ----------------------------------------------------
2489  nCutIndex = 27;
2490  cutVariableValue = normalizedDecayLengthXY;
2491  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2492  //---------------------------------------------------------------------
2493 
2494  // "chi squared per NDF" ----------------------------------------------------
2495  nCutIndex = 28;
2496  cutVariableValue = chi2Vertex;
2497  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2498  //---------------------------------------------------------------------
2499 
2500 
2501 
2502  AliAODRecoDecay* candidateD0toBPlus = (AliAODRecoDecay*)candidateD0;
2503  AliExternalTrackParam pionD0Track;
2504  AliExternalTrackParam kaonD0Track;
2505 
2506  Double_t d0z0DSVert[2], covd0z0DSVert[3], d0DSVert[2];
2507 
2508  pionD0Track.CopyFromVTrack(candidatePion);
2509  pionD0Track.PropagateToDCA(vertexBPlus, bz, 100., d0z0DSVert, covd0z0DSVert);
2510  d0DSVert[0] = d0z0DSVert[0];
2511 
2512  kaonD0Track.CopyFromVTrack(candidateKaon);
2513  kaonD0Track.PropagateToDCA(vertexBPlus, bz, 100., d0z0DSVert, covd0z0DSVert);
2514  d0DSVert[1] = d0z0DSVert[0];
2515 
2516  AliExternalTrackParam D0Track;
2517  D0Track.CopyFromVTrack(candidateD0);
2518  Double_t d0z0D0DSVert[2], covd0z0D0DSVert[3], d0D0DSVert;
2519  motherTrack.PropagateToDCA(vertexBPlus, bz, 100., d0z0D0DSVert, covd0z0D0DSVert);
2520  d0D0DSVert = TMath::Abs(d0z0D0DSVert[0]);
2521 
2522  Double_t impactProductToBPlus = d0DSVert[0] * d0DSVert[1];
2523  Double_t impactProductXYToBPlus = candidateD0toBPlus->ImpParXY(vertexBPlus);
2524 
2525  Double_t pointingAngleToBPlus = candidateD0toBPlus->CosPointingAngle(vertexBPlus);
2526  Double_t d0FirstDaughterToBPlus = TMath::Abs(d0DSVert[0]);
2527  Double_t d0SecondDaughterToBPlus = TMath::Abs(d0DSVert[1]);
2528  Double_t normDecayLengthToBPlus = candidateD0toBPlus->NormalizedDecayLength(vertexBPlus);
2529 
2530  Double_t pseudoProperDecayLengthDSVert = ((vertexD0->GetX() - vertexBPlus->GetX()) * candidateD0->Px() / TMath::Abs(candidateD0->Pt())) + ((vertexD0->GetY() - vertexBPlus->GetY()) * candidateD0->Py() / TMath::Abs(candidateD0->Pt()));
2531  Double_t pseudoProperDecayTimeToBPlus = pseudoProperDecayLengthDSVert * pdgMassMother / ptMother;
2532  Double_t DecayTimeToBPlus = vertexDistance / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother / (momentumMother * momentumMother)) + 1)));
2533 
2534  Double_t phiDSVert = candidateD0->Phi();
2535  Double_t thetaDSVert = candidateD0->Theta();
2536  Double_t covMatrixDSVert[21];
2537  candidateD0->GetCovarianceXYZPxPyPz(covMatrixDSVert);
2538 
2539  cp = TMath::Cos(phiDSVert);
2540  sp = TMath::Sin(phiDSVert);
2541  ct = TMath::Cos(thetaDSVert);
2542  st = TMath::Sin(thetaDSVert);
2543 
2544  errorMomentum = covMatrix[9] * cp * cp * ct * ct // GetCovPxPx
2545  + covMatrix[13] * 2.*cp * sp * ct * ct // GetCovPxPy
2546  + covMatrix[18] * 2.*cp * ct * st // GetCovPxPz
2547  + covMatrix[14] * sp * sp * ct * ct // GetCovPyPy
2548  + covMatrix[19] * 2.*sp * ct * st // GetCovPyPz
2549  + covMatrix[20] * st * st; // GetCovPzPz
2550  Double_t normalizedDecayTimeToBPlus = candidateD0toBPlus->NormalizedDecayLength(vertexBPlus) / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother * errorMomentum * errorMomentum / (momentumMother * momentumMother)) + 1)));
2551 
2552  // "pointingAngleToBPlus" ---------------------------------------------
2553  nCutIndex = 29;
2554  cutVariableValue = pointingAngleToBPlus;
2555  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2556  //---------------------------------------------------------------------
2557 
2558  // "d0MotherToBPlus" --------------------------------------------------
2559  nCutIndex = 30;
2560  cutVariableValue = d0D0DSVert;
2561  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2562  //---------------------------------------------------------------------
2563 
2564  // "d0FirstDaughterToBPlus" -------------------------------------------
2565  nCutIndex = 31;
2566  cutVariableValue = d0FirstDaughterToBPlus;
2567  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2568  //---------------------------------------------------------------------
2569 
2570  // "d0SecondDaughterToBPlus" ------------------------------------------
2571  nCutIndex = 32;
2572  cutVariableValue = d0SecondDaughterToBPlus;
2573  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2574  //---------------------------------------------------------------------
2575 
2576  // "impactProductToBPlus" ---------------------------------------------
2577  nCutIndex = 33;
2578  cutVariableValue = impactProductToBPlus;
2579  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2580  //---------------------------------------------------------------------
2581 
2582  // "impactProductXYToBPlus" -------------------------------------------
2583  nCutIndex = 34;
2584  cutVariableValue = impactProductXYToBPlus;
2585  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2586  //---------------------------------------------------------------------
2587 
2588  // "normDecayLengthToBPlus" -------------------------------------------
2589  nCutIndex = 35;
2590  cutVariableValue = normDecayLengthToBPlus;
2591  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2592  //---------------------------------------------------------------------
2593 
2594  // "pseudoProperDecayTimeToBPlus" -------------------------------------
2595  nCutIndex = 36;
2596  cutVariableValue = pseudoProperDecayTimeToBPlus;
2597  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2598  //---------------------------------------------------------------------
2599 
2600  // "DecayTimeToBPlus" -------------------------------------------------
2601  nCutIndex = 37;
2602  cutVariableValue = DecayTimeToBPlus;
2603  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2604  //---------------------------------------------------------------------
2605 
2606  // "normalizedDecayTimeToBPlus" ---------------------------------------------
2607  nCutIndex = 38;
2608  cutVariableValue = normalizedDecayTimeToBPlus;
2609  if(!ApplyCutOnVariableMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2610  //---------------------------------------------------------------------
2611  }
2612 
2613  return 1;
2614 }
2615 //----------------------------------------------------------------------------------
2616 Int_t AliRDHFCutsBPlustoD0Pi::IsD0forD0ptbinSelectedMVA(TObject* obj, Int_t selectionLevel, AliAODEvent* aod, AliAODVertex *primaryVertex, Double_t bz) {
2617  //
2618  // Apply selection on D0 candidate (for MVA, need to be merged with function for standard analysis)
2619  //
2620 
2621  if (!fCutsRDD0forD0ptbin) {
2622  cout << "Cut matrice not inizialized. Exit..." << endl;
2623  return 0;
2624  }
2625 
2626  AliAODRecoDecayHF2Prong* candidateD0 = (AliAODRecoDecayHF2Prong*)obj;
2627  if (!candidateD0) {
2628  cout << "candidateD0 null" << endl;
2629  return 0;
2630  }
2631 
2632  AliAODTrack *candidatePion = (AliAODTrack*)candidateD0->GetDaughter(0);
2633  if (!candidatePion) {
2634  cout << "candidatePion null 3" << endl;
2635  return 0;
2636  }
2637 
2638  AliAODTrack *candidateKaon = (AliAODTrack*)candidateD0->GetDaughter(1);
2639  if (!candidateKaon) {
2640  cout << "candidateKaon null" << endl;
2641  return 0;
2642  }
2643 
2644  AliAODVertex * vertexD0 = candidateD0->GetSecondaryVtx();
2645  if (!vertexD0) {
2646  cout << "vertexD0 null" << endl;
2647  return 0;
2648  }
2649 
2650  if (!primaryVertex) {
2651  cout << "primaryVertex null" << endl;
2652  return 0;
2653  }
2654 
2655  Bool_t bPassedCut = kFALSE;
2656 
2657  // selection on candidate
2658  if (selectionLevel == AliRDHFCuts::kAll ||
2659  selectionLevel == AliRDHFCuts::kCandidate) {
2660 
2661  Int_t ptbin = PtBinD0forD0ptbin(candidateD0->Pt());
2662  if (ptbin == -1) return -1;
2663 
2664  // D0mass
2665  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2666 
2667  // D0 window - invariant mass
2668  UInt_t prongs[2];
2669  prongs[0] = 211; prongs[1] = 321;
2670  Double_t invMassD0 = candidateD0->InvMass(2, prongs);
2671  Double_t invMassDifference = TMath::Abs(mD0PDG - invMassD0);
2672 
2673  UInt_t prongs2[2];
2674  prongs2[1] = 211; prongs2[0] = 321;
2675  Double_t invMassD02 = candidateD0->InvMass(2, prongs2);
2676  Double_t invMassDifference2 = TMath::Abs(mD0PDG - invMassD02);
2677 
2678  if (invMassDifference > invMassDifference2) invMassDifference = invMassDifference2;
2679 
2680  Double_t pointingAngle = candidateD0->CosPointingAngle();
2681  Double_t dcaMother = candidateD0->GetDCA();
2682  Double_t ptMother = candidateD0->Pt();
2683  Double_t momentumMother = candidateD0->P();
2684  Double_t ptPion = candidatePion->Pt();
2685  Double_t ptKaon = candidateKaon->Pt();
2686 
2687  AliExternalTrackParam motherTrack;
2688  motherTrack.CopyFromVTrack(candidateD0);
2689  Double_t d0z0[2], covd0z0[3], d0[2];
2690  motherTrack.PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
2691  d0[0] = d0z0[0];
2692  Double_t d0Mother = TMath::Abs(d0[0]);
2693  Double_t d0firstTrack = TMath::Abs(candidateD0->Getd0Prong(0));
2694  Double_t d0secondTrack = TMath::Abs(candidateD0->Getd0Prong(1));
2695 
2696  Double_t impactProduct = candidateD0->Prodd0d0();
2697  Double_t impactProductXY = TMath::Abs(candidateD0->ImpParXY());
2698 
2699  Double_t angleBetweenBothDaughters = (candidateKaon->Px() * candidatePion->Px() + candidateKaon->Py() * candidatePion->Py() + candidateKaon->Pz() * candidatePion->Pz()) / (candidateKaon->P() * candidatePion->P());
2700  Double_t angleMotherFirstDaughter = (candidateD0->Px() * candidatePion->Px() + candidateD0->Py() * candidatePion->Py() + candidateD0->Pz() * candidatePion->Pz()) / (candidateD0->P() * candidatePion->P());
2701  Double_t angleMotherSecondDaughter = (candidateD0->Px() * candidateKaon->Px() + candidateD0->Py() * candidateKaon->Py() + candidateD0->Pz() * candidateKaon->Pz()) / (candidateD0->P() * candidateKaon->P());
2702 
2703  Double_t cosThetaStar = candidateD0->CosThetaStar(0, 421, 211, 321);
2704  Double_t vertexDistance = vertexD0->DistanceToVertex(primaryVertex);
2705  Double_t normDecayLength = candidateD0->NormalizedDecayLength();
2706  Double_t pdgMassMother = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2707  Double_t pseudoProperDecayLength = ((vertexD0->GetX() - primaryVertex->GetX()) * candidateD0->Px() / TMath::Abs(candidateD0->Pt())) + ((vertexD0->GetY() - primaryVertex->GetY()) * candidateD0->Py() / TMath::Abs(candidateD0->Pt()));
2708  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother / ptMother;
2709  Double_t decayTime = vertexDistance / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother / (momentumMother * momentumMother)) + 1)));
2710 
2711  Double_t phi = candidateD0->Phi();
2712  Double_t theta = candidateD0->Theta();
2713  Double_t covMatrix[21];
2714  candidateD0->GetCovarianceXYZPxPyPz(covMatrix);
2715 
2716  Double_t cp = TMath::Cos(phi);
2717  Double_t sp = TMath::Sin(phi);
2718  Double_t ct = TMath::Cos(theta);
2719  Double_t st = TMath::Sin(theta);
2720 
2721  Double_t errorMomentum = covMatrix[9] * cp * cp * ct * ct // GetCovPxPx
2722  + covMatrix[13] * 2.*cp * sp * ct * ct // GetCovPxPy
2723  + covMatrix[18] * 2.*cp * ct * st // GetCovPxPz
2724  + covMatrix[14] * sp * sp * ct * ct // GetCovPyPy
2725  + covMatrix[19] * 2.*sp * ct * st // GetCovPyPz
2726  + covMatrix[20] * st * st; // GetCovPzPz
2727  Double_t normalizedDecayTime = candidateD0->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother * errorMomentum * errorMomentum / (momentumMother * momentumMother)) + 1)));
2728 
2729  Double_t cosPointingAngleXY = candidateD0->CosPointingAngleXY();
2730  Double_t distanceXYToVertex = vertexD0->DistanceXYToVertex(primaryVertex);
2731  Double_t normalizedDecayLengthXY = candidateD0->NormalizedDecayLengthXY();
2732  Double_t chi2Vertex = vertexD0->GetChi2perNDF();
2733 
2734  //Topomatic
2735  Double_t dd0pr1 = 0.;
2736  Double_t dd0pr2 = 0.;
2737  Double_t dd0max = 0.;
2738  Double_t dd0min = 0.;
2739  for (Int_t ipr = 0; ipr < 2; ipr++)
2740  {
2741  Double_t diffIP, errdiffIP;
2742  candidateD0->Getd0MeasMinusExpProng(ipr, bz, diffIP, errdiffIP);
2743  Double_t normdd0 = 0.;
2744  if (errdiffIP > 0.) normdd0 = diffIP / errdiffIP;
2745  if (ipr == 0) dd0pr1 = normdd0;
2746  if (ipr == 1) dd0pr2 = normdd0;
2747  }
2748  if (TMath::Abs(dd0pr1) > TMath::Abs(dd0pr2)) {dd0max = dd0pr1; dd0min = dd0pr2;}
2749  else {dd0max = dd0pr2; dd0min = dd0pr1;}
2750 
2751 
2752  // We apply the cuts
2753  Int_t nCutIndex = 0;
2754  Double_t cutVariableValue = 0.0;
2755 
2756  // "inv. mass width [GeV]" --------------------------------------------
2757  nCutIndex = 0;
2758  cutVariableValue = invMassDifference;
2759  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2760  //---------------------------------------------------------------------
2761 
2762  // "delta mass width [GeV]" -------------------------------------------
2763  // nCutIndex = 1; // not used for D0
2764  // cutVariableValue = invMassDelta;
2765  // if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2766  //---------------------------------------------------------------------
2767 
2768  // "pointing angle [Cos(theta)]" --------------------------------------
2769  nCutIndex = 2;
2770  cutVariableValue = pointingAngle;
2771  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2772  //---------------------------------------------------------------------
2773 
2774  // "dca [cm]" ---------------------------------------------------------
2775  nCutIndex = 3;
2776  cutVariableValue = dcaMother;
2777  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2778  //---------------------------------------------------------------------
2779 
2780  // "Pt D0 [GeV/c]" ----------------------------------------------------
2781  nCutIndex = 4;
2782  cutVariableValue = ptMother;
2783  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2784  //---------------------------------------------------------------------
2785 
2786  // "Pt Kaon [GeV/c]" -------------------------------------------------
2787  nCutIndex = 5;
2788  cutVariableValue = ptKaon;
2789  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2790  //---------------------------------------------------------------------
2791 
2792  // "Pt Pion [GeV/c]" --------------------------------------------------
2793  nCutIndex = 6;
2794  cutVariableValue = ptPion;
2795  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2796  //---------------------------------------------------------------------
2797 
2798  // "d0 D0 [cm]" -------------------------------------------------------
2799  nCutIndex = 7;
2800  cutVariableValue = d0Mother;
2801  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2802  //---------------------------------------------------------------------
2803 
2804  // "d0 Kaon [cm]"-----------------------------------------------------
2805  nCutIndex = 8;
2806  cutVariableValue = d0firstTrack;
2807  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2808  //---------------------------------------------------------------------
2809 
2810  // "d0 Pion [cm]" -----------------------------------------------------
2811  nCutIndex = 9;
2812  cutVariableValue = d0secondTrack;
2813  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2814  //---------------------------------------------------------------------
2815 
2816  // "d0d0 [cm^2]" ------------------------------------------------------
2817  nCutIndex = 10;
2818  cutVariableValue = impactProduct;
2819  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2820  //---------------------------------------------------------------------
2821 
2822  // "d0d0 XY [cm^2]" ---------------------------------------------------
2823  nCutIndex = 11;
2824  cutVariableValue = impactProductXY;
2825  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2826  //---------------------------------------------------------------------
2827 
2828  // "angle between both daughters" -------------------------------------
2829  nCutIndex = 12;
2830  cutVariableValue = angleBetweenBothDaughters;
2831  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2832  //---------------------------------------------------------------------
2833 
2834  // "angle mother with first daughter" ---------------------------------
2835  nCutIndex = 13;
2836  cutVariableValue = angleMotherFirstDaughter;
2837  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2838  //---------------------------------------------------------------------
2839 
2840  // "angle mother with second daughter" --------------------------------
2841  nCutIndex = 14;
2842  cutVariableValue = angleMotherSecondDaughter;
2843  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2844  //---------------------------------------------------------------------
2845 
2846  // "cosThetaStar" -----------------------------------------------------
2847  nCutIndex = 15;
2848  cutVariableValue = cosThetaStar;
2849  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2850  //---------------------------------------------------------------------
2851 
2852  // "vertexDistance" ---------------------------------------------------
2853  nCutIndex = 16;
2854  cutVariableValue = vertexDistance;
2855  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2856  //---------------------------------------------------------------------
2857 
2858  // "pseudoProperDecayTime" --------------------------------------------
2859  nCutIndex = 17;
2860  cutVariableValue = pseudoProperDecayTime;
2861  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2862  //---------------------------------------------------------------------
2863 
2864  // "DecayTime" --------------------------------------------------------
2865  nCutIndex = 18;
2866  cutVariableValue = decayTime;
2867  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2868  //---------------------------------------------------------------------
2869 
2870  // "normalizedDecayTime" ----------------------------------------------------
2871  nCutIndex = 19;
2872  cutVariableValue = normalizedDecayTime;
2873  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2874  //---------------------------------------------------------------------
2875 
2876  // "normDecayLength" --------------------------------------------------
2877  nCutIndex = 20;
2878  cutVariableValue = normDecayLength;
2879  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2880  //---------------------------------------------------------------------
2881 
2882  // "topomatic first daughter" -----------------------------------------
2883  nCutIndex = 21;
2884  cutVariableValue = dd0pr1;
2885  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2886  //---------------------------------------------------------------------
2887 
2888  // "topomatic second daughter" ----------------------------------------
2889  nCutIndex = 22;
2890  cutVariableValue = dd0pr2;
2891  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2892  //---------------------------------------------------------------------
2893 
2894  // "topomatic max" ----------------------------------------------------
2895  nCutIndex = 23;
2896  cutVariableValue = dd0max;
2897  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2898  //---------------------------------------------------------------------
2899 
2900  // "topomatic min" ----------------------------------------------------
2901  nCutIndex = 24;
2902  cutVariableValue = dd0min;
2903  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2904  //---------------------------------------------------------------------
2905 
2906  // "pointing angle XY" ----------------------------------------------------
2907  nCutIndex = 25;
2908  cutVariableValue = cosPointingAngleXY;
2909  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2910  //---------------------------------------------------------------------
2911 
2912  // "vertex distance XY" ----------------------------------------------------
2913  nCutIndex = 26;
2914  cutVariableValue = distanceXYToVertex;
2915  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2916  //---------------------------------------------------------------------
2917 
2918  // "normalized decay length XY" ----------------------------------------------------
2919  nCutIndex = 27;
2920  cutVariableValue = normalizedDecayLengthXY;
2921  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2922  //---------------------------------------------------------------------
2923 
2924  // "chi squared per NDF" ----------------------------------------------------
2925  nCutIndex = 28;
2926  cutVariableValue = chi2Vertex;
2927  if(!ApplyCutOnVariableD0forD0ptbinMVA(nCutIndex, ptbin, cutVariableValue)) return 0;
2928  //---------------------------------------------------------------------
2929  }
2930 
2931  if(selectionLevel == AliRDHFCuts::kAll || selectionLevel == AliRDHFCuts::kTracks){
2932 
2933  if(!AreDaughtersSelected(candidateD0,aod)) return 0;
2934 
2935  AliExternalTrackParam particleTrack;
2936  particleTrack.CopyFromVTrack(candidatePion);
2937  Double_t d0[2], covd0[3];
2938  particleTrack.PropagateToDCA(primaryVertex, bz, 100., d0, covd0);
2939 
2940  if (TMath::Abs(d0[0]) < GetMind0D0FirstDaughter()) return 0;
2941 
2942  AliExternalTrackParam particleTrack2;
2943  particleTrack2.CopyFromVTrack(candidateKaon);
2944  Double_t d02[2], covd02[3];
2945  particleTrack2.PropagateToDCA(primaryVertex, bz, 100., d02, covd02);
2946 
2947  if (TMath::Abs(d0[0]) < GetMind0D0SecondDaughter()) return 0;
2948 
2949  if (UseFilterBitD0FirstDaughter() == kTRUE) {
2950  if (!(candidatePion->TestFilterMask(BIT(GetFilterBitD0FirstDaughter())))) return 0;
2951  if (!(candidateKaon->TestFilterMask(BIT(GetFilterBitD0SecondDaughter())))) return 0;
2952  }
2953  }
2954 
2955  return 1;
2956 }
2957 //----------------------------------------------------------------------------------
2958 Int_t AliRDHFCutsBPlustoD0Pi::IsBplusPionSelectedMVA(TObject* obj,Int_t selectionLevel, AliAODEvent* aod, AliAODVertex *primaryVertex, Double_t bz) {
2959  //
2960  // Apply selection on D0 candidate.
2961  //
2962 
2963  AliAODTrack* candidatePion = (AliAODTrack*)obj;
2964  if (!candidatePion) {
2965  cout << "candidatePion null" << endl;
2966  return 0;
2967  }
2968 
2969  if(selectionLevel == AliRDHFCuts::kAll || selectionLevel == AliRDHFCuts::kTracks){
2970 
2971  //quick track quality pre selection to speed things up.
2972  //The first 4 cuts are also in AliESDtrackCuts object
2973  if (candidatePion->GetITSNcls() < 1) return 0;
2974  if (candidatePion->GetTPCNcls() < 1) return 0;
2975  if (candidatePion->GetStatus()&AliESDtrack::kITSpureSA) return 0;
2976  if (!(candidatePion->GetStatus()&AliESDtrack::kITSin)) return 0;
2977  if (candidatePion->GetID() < 0) return 0;
2978  Double_t covtest[21];
2979  if (!candidatePion->GetCovarianceXYZPxPyPz(covtest)) return 0;
2980 
2981  Double_t pos[3],cov[6];
2982  primaryVertex->GetXYZ(pos);
2983  primaryVertex->GetCovarianceMatrix(cov);
2984  const AliESDVertex vESD(pos,cov,100.,100);
2985 
2986  if(!IsDaughterSelected(candidatePion,&vESD,fTrackCutsSoftPi,aod)) return 0;
2987 
2988  AliExternalTrackParam particleTrack;
2989  particleTrack.CopyFromVTrack(candidatePion);
2990  Double_t d0[2], covd0[3];
2991  particleTrack.PropagateToDCA(primaryVertex, bz, 100., d0, covd0);
2992 
2993  if (TMath::Abs(d0[0]) < GetMind0BPlusPion()) return 0;
2994 
2995  if (UseFilterBitBPlusPion() == kTRUE) {
2996  if (!(candidatePion->TestFilterMask(BIT(GetFilterBitBPlusPion())))) return 0;
2997  }
2998 
2999  }
3000 
3001  if(selectionLevel == AliRDHFCuts::kAll || selectionLevel == AliRDHFCuts::kPID){
3002  if (!(SelectPID(candidatePion, 2))) return 0;
3003  }
3004 
3005  return 1;
3006 
3007 }
3008 //---------------------------------------------------------------------------
3009 Int_t AliRDHFCutsBPlustoD0Pi::IsD0SelectedPreRecVtxMVA(AliAODRecoDecayHF2Prong* d, AliAODTrack* pion, AliAODVertex *primaryVertex, Double_t bz, Int_t selLevel = 0)
3010 {
3011  //
3012  // Preselection before reconstruction Bplus vertex to save time.
3013  //
3014  // selLevel =
3015  // 0: Everything
3016  // 1: PID of D0 daughters (using charge pion track. Always true when PID is off)
3017  // 2: Invariant mass cut D0 (Update earlier D0forD0ptbin cut with correct D0/D0bar mass)
3018  // 3: Bplus invariant mass cut without reconstruction vertex
3019  // 4: Everything except Bplus invariant mass cut
3020  //
3021  // returnvalue:
3022  // 1: selected
3023  // 0: selected (but Wrong-Sign pair)
3024  // -1: all other (rejected)
3025 
3026  if (!fCutsRD) {
3027  cout << "Cut matrice not inizialized. Exit..." << endl;
3028  return -1;
3029  }
3030 
3031  if (!d) {
3032  cout << "candidate D0 null" << endl;
3033  return -1;
3034  }
3035 
3036  if (!pion) {
3037  cout << "candidate Pion null" << endl;
3038  return -1;
3039  }
3040 
3041  if (!primaryVertex) {
3042  cout << "primaryVertex null" << endl;
3043  return -1;
3044  }
3045 
3046  if(selLevel > 4 || selLevel < 0){
3047  AliWarning(Form("selLevel %d is not supported. return -1.",selLevel));
3048  return -1;
3049  }
3050 
3051  Int_t returnvaluePID = 1;
3052  Int_t returnvalueD0 = 1;
3053  Int_t returnvalueBplus = 1;
3054 
3055  // PID D0 daughters (check if the pions have the opposite charge)
3056  if(selLevel == 0 || selLevel == 1 || selLevel == 4){
3057  AliAODTrack* track0 = (AliAODTrack*)d->GetDaughter(0);
3058  AliAODTrack* track1 = (AliAODTrack*)d->GetDaughter(1);
3059 
3060  if (pion->Charge() == -1) {
3061  if ((SelectPID(track0, 2)) && (SelectPID(track1, 3))) returnvaluePID = 1;
3062  else if ((SelectPID(track0, 3)) && (SelectPID(track1, 2))) returnvaluePID = 0;
3063  else returnvaluePID = -1;
3064  } else if (pion->Charge() == 1) {
3065  if ((SelectPID(track0, 3)) && (SelectPID(track1, 2))) returnvaluePID = 1;
3066  else if ((SelectPID(track0, 2)) && (SelectPID(track1, 3))) returnvaluePID = 0;
3067  else returnvaluePID = -1;
3068  } else {
3069  returnvaluePID = -1;
3070  }
3071  }
3072 
3073  // D0 window - invariant mass cut
3074  if(selLevel == 0 || selLevel == 2 || selLevel == 4){
3075  Double_t pdgMassD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
3076  Int_t ptbin = PtBinD0forD0ptbin(d->Pt());
3077  if(ptbin == -1) returnvalueD0 = -1;
3078 
3079  if (pion->Charge() == -1) {
3080  Double_t invMassDifference = TMath::Abs(d->InvMassD0() - pdgMassD0);
3081  if(!ApplyCutOnVariableD0forD0ptbinMVA(0, ptbin, invMassDifference)) returnvalueD0 = -1;
3082  } else if (pion->Charge() == 1) {
3083  Double_t invMassDifference = TMath::Abs(d->InvMassD0bar() - pdgMassD0);
3084  if(!ApplyCutOnVariableD0forD0ptbinMVA(0, ptbin, invMassDifference)) returnvalueD0 = -1;
3085  } else {
3086  returnvalueD0 = -1;
3087  }
3088  }
3089 
3090  // Bplus invariant mass window (using PV instead of SV, to speed thing up)
3091  if(selLevel == 0 || selLevel == 3){
3092  AliExternalTrackParam firstTrack;
3093  firstTrack.CopyFromVTrack(pion);
3094  AliExternalTrackParam secondTrack;
3095  secondTrack.CopyFromVTrack(d);
3096 
3097  Double_t d0z0[2], covd0z0[3], d0[2], d0err[2];
3098  firstTrack.PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
3099  d0[0] = d0z0[0]; d0err[0] = TMath::Sqrt(covd0z0[0]);
3100  secondTrack.PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
3101  d0[1] = d0z0[0]; d0err[1] = TMath::Sqrt(covd0z0[0]);
3102 
3103  Double_t px[2],py[2],pz[2],momentum[3];
3104  firstTrack.GetPxPyPz(momentum);
3105  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
3106  secondTrack.GetPxPyPz(momentum);
3107  px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
3108 
3109  Double_t xdummy = 0., ydummy = 0.;
3110  Double_t dca = secondTrack.GetDCA(&firstTrack, bz, xdummy, ydummy);
3111 
3112  AliAODRecoDecayHF2Prong trackBPlus(primaryVertex, px, py, pz, d0, d0err, dca);
3113 
3114  UInt_t pdg2[2] = {211,421};
3115  Double_t invMassBPlus = trackBPlus.InvMass(2,pdg2);
3116  Double_t mBPlusPDG = TDatabasePDG::Instance()->GetParticle(521)->Mass();
3117  Double_t invMassDifference = TMath::Abs(mBPlusPDG - invMassBPlus) - 0.1; //-0.1 to be on the safe side as invmass is not completely correct here, will be recomputed in IsSelected()
3118 
3119  Int_t ptbin = PtBin(trackBPlus.Pt());
3120  if (ptbin == -1) returnvalueBplus = -1;
3121  else{
3122  if(!ApplyCutOnVariableMVA(39, ptbin, invMassDifference)) returnvalueBplus = -1;
3123  }
3124  }
3125 
3126  if(returnvaluePID == 1 && returnvalueD0 == 1 && returnvalueBplus == 1) return 1;
3127  else if(returnvaluePID == 0 && returnvalueD0 == 1 && returnvalueBplus == 1) return 0;
3128  else return -1;
3129 
3130 }
3131 //----------------------------------------------------------------------------------
3133 {
3134  //
3135  // Should be improved (for example pT dependence)
3136  //
3137 
3138  if (y > GetFiducialYCut()) return kFALSE;
3139 
3140  return kTRUE;
3141 }
3142 //_______________________________________________________________________________-
3144 {
3145  //
3146  // PID method, n sigma approach default // not used for BPlus, done seperately for each daughter
3147  //
3148 
3149  // AliAODRecoDecayHF2Prong* dstar = (AliAODRecoDecayHF2Prong*)obj;
3150  // if(!dstar){
3151  // cout<<"AliAODRecoDecayHF2Prong null"<<endl;
3152  // return 0;
3153  // }
3154 
3155  // if(!fUsePID || dstar->Pt() > fMaxPtPid) return 3;
3156 
3157  // AliAODRecoDecayHF2Prong* d0 = (AliAODRecoDecayHF2Prong*)dstar->Get2Prong();
3158  // if(!d0){
3159  // cout<<"AliAODRecoDecayHF2Prong null"<<endl;
3160  // return 0;
3161  // }
3162 
3163  // // here the PID
3164  // AliAODTrack *pos = (AliAODTrack*)dstar->Get2Prong()->GetDaughter(0);
3165  // AliAODTrack *neg = (AliAODTrack*)dstar->Get2Prong()->GetDaughter(1);
3166 
3167  // if (dstar->Charge()>0){
3168  // if(!SelectPID(pos,2)) return 0;//pion+
3169  // if(!SelectPID(neg,3)) return 0;//kaon-
3170  // }else{
3171  // if(!SelectPID(pos,3)) return 0;//kaon+
3172  // if(!SelectPID(neg,2)) return 0;//pion-
3173  // }
3174 
3175  // if ((fPidHF->GetMatch() == 10 || fPidHF->GetMatch() == 11) && fPidHF->GetITS()) { //ITS n sigma band
3176  // AliAODTrack *softPion = (AliAODTrack*) dstar->GetBachelor();
3177 
3178  // if (fPidHF->CheckBands(AliPID::kPion, AliPIDResponse::kITS, softPion) == -1) {
3179  // return 0;
3180  // }
3181  // }
3182 
3183  return 3;
3184 }
3185 //_______________________________________________________________________________-
3187 {
3188  //
3189  // here the PID
3190 
3191  Bool_t isParticle = kTRUE;
3192  if (!fUsePID) return isParticle;
3193  Int_t match = fPidHF->GetMatch();
3194 
3195  if (match == 1) { //n-sigma
3196  Bool_t TPCon = TMath::Abs(2) > 1e-4 ? kTRUE : kFALSE;
3197  Bool_t TOFon = TMath::Abs(3) > 1e-4 ? kTRUE : kFALSE;
3198 
3199  Bool_t isTPC = kTRUE;
3200  Bool_t isTOF = kTRUE;
3201 
3202  if (TPCon) { //TPC
3203  if (fPidHF->CheckStatus(track, "TPC")) {
3204  if (type == 2) isTPC = fPidHF->IsPionRaw(track, "TPC");
3205  if (type == 3) isTPC = fPidHF->IsKaonRaw(track, "TPC");
3206  }
3207  }
3208  if (TOFon) { //TOF
3209  if (fPidHF->CheckStatus(track, "TOF")) {
3210  if (type == 2) isTOF = fPidHF->IsPionRaw(track, "TOF");
3211  if (type == 3) isTOF = fPidHF->IsKaonRaw(track, "TOF");
3212  }
3213  }
3214 
3215  //--------------------------------
3216  // cut on high momentum in the TPC
3217  //--------------------------------
3218  Double_t pPIDcut = track->P();
3219  if (pPIDcut > fTPCflag) isTPC = 1;
3220 
3221  isParticle = isTPC && isTOF;
3222  }
3223 
3224  if (match == 2) { //bayesian
3225  //Double_t priors[5]={0.01,0.001,0.3,0.3,0.3};
3226  Double_t prob[5] = {1., 1., 1., 1., 1.};
3227 
3228  //fPidHF->SetPriors(priors,5);
3229  // fPidHF->BayesianProbability(track,prob);
3230 
3231  Double_t max = 0.;
3232  Int_t k = -1;
3233  for (Int_t i = 0; i < 5; i++) {
3234  if (prob[i] > max) {k = i; max = prob[i];}
3235  }
3236  isParticle = Bool_t(k == type);
3237  }
3238 
3239  if (match == 10 || match == 11) { //Assymetric PID using histograms
3240  Int_t checkTPC = fPidHF->CheckBands((AliPID::EParticleType) type, AliPIDResponse::kTPC, track);
3241  Int_t checkTOF = fPidHF->CheckBands((AliPID::EParticleType) type, AliPIDResponse::kTOF, track);
3242 
3243  isParticle = checkTPC >= 0 && checkTOF >= 0 ? kTRUE : kFALSE; //Standard: both at least compatible
3244  if (match == 11) { //Extra requirement: at least one identified
3245  isParticle = isParticle && checkTPC + checkTOF >= 1 ? kTRUE : kFALSE;
3246  }
3247  }
3248 
3249  if (match == 12) { //Circular cut
3250  Double_t nSigmaTOF = 0;
3251  Double_t nSigmaTPC = 0;
3252 
3253  Double_t radius = fCircRadius;
3254 
3255  isParticle = kTRUE;
3256  if (radius > 0) {
3257  Int_t TPCok = fPidHF->GetnSigmaTPC(track, type, nSigmaTPC);
3258  Int_t TOFok = fPidHF->GetnSigmaTOF(track, type, nSigmaTOF);
3259  if (TPCok != -1 && TOFok != -1) {
3260  //Both detectors gave PID information
3261  isParticle = TMath::Sqrt(nSigmaTPC * nSigmaTPC + nSigmaTOF * nSigmaTOF) <= radius ? kTRUE : kFALSE;
3262  }
3263  else {
3264  //Maximum one detector gave PID information
3265  if (TPCok != -1) {
3266  isParticle = nSigmaTPC <= radius ? kTRUE : kFALSE;
3267  }
3268  if (TOFok != -1) {
3269  isParticle = nSigmaTOF <= radius ? kTRUE : kFALSE;
3270  }
3271  }
3272  }
3273  }
3274 
3275  return isParticle;
3276 
3277 }
3278 //-------------------------------------------------------------------------------------
3280 {
3284 
3285  AliAODRecoDecayHF2Prong * D0 = (AliAODRecoDecayHF2Prong*)BPlus->GetDaughter(1);
3286 
3287  Double_t e[3] = {0};
3288  if (BPlus->Charge() == -1)
3289  {
3290  e[0] = D0->EProng(0, 211);
3291  e[1] = D0->EProng(1, 321);
3292  } else if (BPlus->Charge() == 1) {
3293  e[0] = D0->EProng(0, 321);
3294  e[1] = D0->EProng(1, 211);
3295  }
3296  e[2] = BPlus->EProng(0, 211);
3297 
3298  Double_t esum = e[0] + e[1] + e[2];
3299  Double_t invMassBPlus = TMath::Sqrt(esum * esum - BPlus->P2());
3300 
3301  Double_t invMassD0 = -1;
3302 
3303  if (BPlus->Charge() == -1) {invMassD0 = D0->InvMassD0();}
3304  else {invMassD0 = D0->InvMassD0bar();}
3305  if (invMassD0 == -1) {cout << "wrong invmass delta D0 BPlus" << endl;}
3306 
3307  return invMassBPlus - invMassD0;
3308 }
3309 
3310 
3311 //---------------------------------------------------------------------------
3312 //
3313 // DO for D0 pt bin functions
3314 //
3315 //---------------------------------------------------------------------------
3316 
3318  //
3319  // store the cuts
3320  //
3321 
3322  if (nVars != fnVarsD0forD0ptbin) {
3323  printf("Wrong number of variables: it has to be %d\n", fnVarsD0forD0ptbin);
3324  AliFatal("exiting");
3325  }
3326  if (nPtBins != fnPtBinsD0forD0ptbin) {
3327  printf("Wrong number of pt bins: it has to be %d\n", fnPtBinsD0forD0ptbin);
3328  AliFatal("exiting");
3329  }
3330 
3332 
3333 
3334  for (Int_t iv = 0; iv < fnVarsD0forD0ptbin; iv++)
3335  {
3336  for (Int_t ib = 0; ib < fnPtBinsD0forD0ptbin; ib++)
3337  {
3338  //check
3339 
3341  {
3342  cout << "Overflow, exit..." << endl;
3343  return;
3344  }
3345 
3346  fCutsRDD0forD0ptbin[GetGlobalIndexD0forD0ptbin(iv, ib)] = cutsRDD0forD0ptbin[iv][ib];
3347 
3348  }
3349  }
3350 
3351  return;
3352 }
3353 //---------------------------------------------------------------------------
3354 void AliRDHFCutsBPlustoD0Pi::SetCutsD0forD0ptbin(Int_t glIndex, Float_t *cutsRDD0forD0ptbin) {
3355  //
3356  // store the cuts
3357  //
3358  if (glIndex != fGlobalIndexD0forD0ptbin) {
3359  cout << "Wrong array size: it has to be " << fGlobalIndexD0forD0ptbin << endl;
3360  AliFatal("exiting");
3361  }
3363 
3364  for (Int_t iGl = 0; iGl < fGlobalIndexD0forD0ptbin; iGl++) {
3365  fCutsRDD0forD0ptbin[iGl] = cutsRDD0forD0ptbin[iGl];
3366  }
3367  return;
3368 }
3369 //---------------------------------------------------------------------------
3371  //
3372  //give the pt bin where the pt lies.
3373  //
3374  Int_t ptbin = -1;
3375  if (pt < fPtBinLimitsD0forD0ptbin[0])return ptbin;
3376  for (Int_t i = 0; i < fnPtBinsD0forD0ptbin; i++) {
3377  if (pt < fPtBinLimitsD0forD0ptbin[i + 1]) {
3378  ptbin = i;
3379  break;
3380  }
3381  }
3382  return ptbin;
3383 }
3384 //---------------------------------------------------------------------------
3386  // Set the pt bins
3387 
3389  delete [] fPtBinLimitsD0forD0ptbin;
3390  fPtBinLimitsD0forD0ptbin = NULL;
3391  printf("Changing the pt bins\n");
3392  }
3393 
3394  if (nPtBinLimits != fnPtBinsD0forD0ptbin + 1) {
3395  cout << "Warning: ptBinLimits dimension " << nPtBinLimits << " != nPtBins+1 (" << fnPtBinsD0forD0ptbin + 1 << ")\nSetting nPtBins to " << nPtBinLimits - 1 << endl;
3396  SetNPtBinsD0forD0ptbin(nPtBinLimits - 1);
3397  }
3398 
3399  fnPtBinLimitsD0forD0ptbin = nPtBinLimits;
3401  //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
3403  for (Int_t ib = 0; ib < nPtBinLimits; ib++) fPtBinLimitsD0forD0ptbin[ib] = ptBinLimits[ib];
3404  for (Int_t ib = 0; ib < nPtBinLimits; ib++) std::cout << "limit " << ib << " = " << fPtBinLimitsD0forD0ptbin[ib] << std::endl;
3405 
3406  return;
3407 }
3408 //---------------------------------------------------------------------------
3409 Int_t AliRDHFCutsBPlustoD0Pi::ApplyCutOnVariable(Int_t nCutIndex, Int_t ptbin, Float_t cutVariableValue, Bool_t bCutArray[68]) {
3410 
3411  if (GetIsCutUsed(nCutIndex, ptbin) == kTRUE)
3412  {
3413  Bool_t bCut = kFALSE;
3414  if (GetIsUpperCut(nCutIndex) == kTRUE)
3415  {
3416  if (cutVariableValue > fCutsRD[GetGlobalIndex(nCutIndex, ptbin)]) bCut = kTRUE;
3417  } else
3418  {
3419  if (cutVariableValue < fCutsRD[GetGlobalIndex(nCutIndex, ptbin)]) bCut = kTRUE;
3420  }
3421  if (bCut == kTRUE) {bCutArray[nCutIndex] = 1; return 0;}
3422  }
3423  return 1;
3424 }
3425 //---------------------------------------------------------------------------
3426 Int_t AliRDHFCutsBPlustoD0Pi::ApplyCutOnVariableD0forD0ptbin(Int_t nCutIndex, Int_t ptbin, Float_t cutVariableValue, Bool_t bCutArray[29]) {
3427 
3428  // std::cout << "index: " << nCutIndex << ", ptbin: " << ptbin << ", used: " << GetIsCutUsedD0forD0ptbin(nCutIndex,ptbin) << ", upper: " << GetIsUpperCutD0forD0ptbin(nCutIndex) << ", value: " << cutVariableValue << ", cutvalue: " << fCutsRDD0forD0ptbin[GetGlobalIndexD0forD0ptbin(nCutIndex,ptbin)] << std::endl;
3429  if (GetIsCutUsedD0forD0ptbin(nCutIndex, ptbin) == kTRUE)
3430  {
3431  Bool_t bCut = kFALSE;
3432  if (GetIsUpperCutD0forD0ptbin(nCutIndex) == kTRUE)
3433  {
3434  if (cutVariableValue > fCutsRDD0forD0ptbin[GetGlobalIndexD0forD0ptbin(nCutIndex, ptbin)]) bCut = kTRUE;
3435  } else
3436  {
3437  if (cutVariableValue < fCutsRDD0forD0ptbin[GetGlobalIndexD0forD0ptbin(nCutIndex, ptbin)]) bCut = kTRUE;
3438  }
3439  if (bCut == kTRUE) {bCutArray[nCutIndex] = 1; return 0;}
3440  }
3441  return 1;
3442 }
3443 //---------------------------------------------------------------------------
3445 
3446  if (GetIsCutUsed(nCutIndex, ptbin) == kTRUE)
3447  {
3448  if (GetIsUpperCut(nCutIndex) == kTRUE)
3449  {
3450  if (cutVariableValue > fCutsRD[GetGlobalIndex(nCutIndex, ptbin)]) return 0;
3451  } else
3452  {
3453  if (cutVariableValue < fCutsRD[GetGlobalIndex(nCutIndex, ptbin)]) return 0;
3454  }
3455  }
3456  return 1;
3457 }
3458 //---------------------------------------------------------------------------
3460 
3461  if (GetIsCutUsedD0forD0ptbin(nCutIndex, ptbin) == kTRUE)
3462  {
3463  if (GetIsUpperCutD0forD0ptbin(nCutIndex) == kTRUE)
3464  {
3465  if (cutVariableValue > fCutsRDD0forD0ptbin[GetGlobalIndexD0forD0ptbin(nCutIndex, ptbin)]) return 0;
3466  } else
3467  {
3468  if (cutVariableValue < fCutsRDD0forD0ptbin[GetGlobalIndexD0forD0ptbin(nCutIndex, ptbin)]) return 0;
3469  }
3470  }
3471  return 1;
3472 }
3473 //---------------------------------------------------------------------------
3475  // Set the variable names
3476 
3477  if (fVarNamesD0forD0ptbin) {
3478  delete [] fVarNamesD0forD0ptbin;
3479  fVarNamesD0forD0ptbin = NULL;
3480  //printf("Changing the variable names\n");
3481  }
3482  if (nVars != fnVarsD0forD0ptbin) {
3483  printf("Wrong number of variables: it has to be %d\n", fnVarsD0forD0ptbin);
3484  return;
3485  }
3486  //fnVars=nVars;
3487  fVarNamesD0forD0ptbin = new TString[nVars];
3488  fIsUpperCutD0forD0ptbin = new Bool_t[nVars];
3489  for (Int_t iv = 0; iv < nVars; iv++) {
3490  fVarNamesD0forD0ptbin[iv] = varNames[iv];
3491  fIsUpperCutD0forD0ptbin[iv] = isUpperCut[iv];
3492  }
3493 
3494  return;
3495 }
3496 //---------------------------------------------------------------------------
3498  // Set the use of all cuts to kFALSE and the cut value to zero. This function has to be called after setting the pt bins.
3499 
3500  if (fIsCutUsed) {
3501  delete [] fIsCutUsed;
3502  fIsCutUsed = NULL;
3503  }
3504 
3505  fIsCutUsed = new Bool_t[(GetNPtBins()) * (GetNVars())];
3506  for (Int_t iv = 0; iv < (GetNPtBins()) * (GetNVars()); iv++) {
3507  fIsCutUsed[iv] = kFALSE;
3508  }
3509 
3510  if (!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
3511 
3512 
3513  for (Int_t iv = 0; iv < fnVars; iv++) {
3514 
3515  for (Int_t ib = 0; ib < fnPtBins; ib++) {
3516 
3517  //check
3518  if (GetGlobalIndex(iv, ib) >= fGlobalIndex) {
3519  cout << "Overflow, exit..." << endl;
3520  return;
3521  }
3522 
3523  fCutsRD[GetGlobalIndex(iv, ib)] = 0;
3524 
3525  }
3526  }
3527 
3528  // D0 for D0 pt bin
3529 
3530  if (fIsCutUsedD0forD0ptbin) {
3531  delete [] fIsCutUsedD0forD0ptbin;
3532  fIsCutUsedD0forD0ptbin = NULL;
3533  }
3534 
3536  for (Int_t iv = 0; iv < (GetNPtBinsD0forD0ptbin()) * (GetNVarsD0forD0ptbin()); iv++) {
3537  fIsCutUsedD0forD0ptbin[iv] = kFALSE;
3538 
3539  }
3540 
3542 
3543  for (Int_t iv = 0; iv < fnVarsD0forD0ptbin; iv++)
3544  {
3545  for (Int_t ib = 0; ib < fnPtBinsD0forD0ptbin; ib++)
3546  {
3547  //check
3548 
3550  {
3551  cout << "Overflow, exit..." << endl;
3552  return;
3553  }
3554 
3556  }
3557  }
3558 
3559  return;
3560 }
3561 //---------------------------------------------------------------------------
3563  // Set the cut value and direction
3564 
3565  fIsCutUsed[GetGlobalIndex(nCutIndex, ptBin)] = kTRUE;
3566  fIsUpperCut[nCutIndex] = cutDirection;
3567  fCutsRD[GetGlobalIndex(nCutIndex, ptBin)] = cutValue;
3568 
3569  return;
3570 }
3571 //---------------------------------------------------------------------------
3573  // Set the cut value and direction
3574 
3575  fIsCutUsedD0forD0ptbin[GetGlobalIndexD0forD0ptbin(nCutIndex, ptBin)] = kTRUE;
3576  fIsUpperCutD0forD0ptbin[nCutIndex] = cutDirection;
3577  fCutsRDD0forD0ptbin[GetGlobalIndexD0forD0ptbin(nCutIndex, ptBin)] = cutValue;
3578 
3579  return;
3580 }
3581 //---------------------------------------------------------------------------
3583  //
3584  // Initialize the cuts for optimization
3585  //
3586 
3587  fnVariablesForCutOptimization = nVariables;
3588  fnCutsForOptimization = nCutsForOptimization;
3590 
3591 
3593  {
3596  }
3598 
3600  {
3601  delete [] fCutIndexForCutOptimization;
3602  fCutIndexForCutOptimization = nullptr;
3603  }
3605 
3607  {
3608  delete [] fSigmaForCutOptimization;
3609  fSigmaForCutOptimization = nullptr;
3610  }
3612 
3614 
3615 
3616  for (Int_t iv = 0; iv < fnVariablesForCutOptimization; iv++)
3617  {
3618  for (Int_t ib = 0; ib < fnPtBins; ib++)
3619  {
3620  for (Int_t ic = 0; ic < fnCutsForOptimization; ++ic)
3621  {
3623  {
3624  cout << "Overflow, exit..." << endl;
3625  return;
3626  }
3627 
3629  }
3630  }
3631  }
3632 
3633  return;
3634 }
3635 //---------------------------------------------------------------------------
3636 void AliRDHFCutsBPlustoD0Pi::SetCutsForCutOptimization(Int_t glIndex, Float_t *cutsRDForCutOptimization) {
3637  //
3638  // store the cuts
3639  //
3640  if (glIndex != fGlobalIndexCutOptimization) {
3641  cout << "Wrong array size: it has to be " << fGlobalIndexCutOptimization << endl;
3642  AliFatal("exiting");
3643  }
3645 
3646  for (Int_t iGl = 0; iGl < fGlobalIndexCutOptimization; iGl++) {
3647  fCutsRDForCutOptimization[iGl] = cutsRDForCutOptimization[iGl];
3648  }
3649  return;
3650 }
3651 //---------------------------------------------------------------------------
3653  // Set the cut value and direction
3654  if (nVariable > GetnVariablesForCutOptimization() || nVariable < 0) cout << "wrong variable number for GetNCutsForOptimization" << endl;
3655 
3656  fIsUpperCutForCutOptimization[nVariable] = cutDirection;
3657  fCutIndexForCutOptimization[nVariable] = nCutIndex;
3658 
3659  for (Int_t iCut = 0; iCut < fnCutsForOptimization; ++iCut)
3660  {
3661  fCutsRDForCutOptimization[GetGlobalIndexForCutOptimization(iCut, nVariable, ptBin)] = cutValues[iCut];
3662  }
3663 
3664  return;
3665 }
Bool_t GetIsCutUsed(Int_t nCutIndex, Int_t ptbin) const
Double_t NormalizedDecayLengthXY() const
Int_t fIsSelectedCuts
fix the daughter track references
Definition: AliRDHFCuts.h:476
Int_t IsD0FromBPlusSelected(Double_t ptBPlus, TObject *obj, Int_t selectionLevel, AliAODEvent *aod, Bool_t bCutArray[68])
Double_t NormalizedDecayLength() const
Int_t GetNVars() const
Definition: AliRDHFCuts.h:264
double Double_t
Definition: External.C:58
Bool_t GetIsCutUsedD0forD0ptbin(Int_t nCutIndex, Int_t ptbin) const
Bool_t GetIsUpperCutD0forD0ptbin(Int_t nCutIndex)
void SetHardSelectionArrayITSD0SecondDaughter(const Bool_t array[7]=0)
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
Int_t GetnSigmaTOF(AliAODTrack *track, Int_t species, Double_t &sigma) const
Int_t GetGlobalIndexD0forD0ptbin(Int_t iVar, Int_t iPtBin) const
Bool_t * fIsUpperCut
Definition: AliRDHFCuts.h:454
Bool_t IsKaonRaw(AliAODTrack *track, TString detector) const
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
AliESDtrackCuts * fTrackCutsSoftPi
Int_t GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &sigma) const
Double_t DeltaInvMassBPlusKpipi(AliAODRecoDecayHF2Prong *BPlus) const
Double_t bz
Double_t ImpParXY() const
void SetVarNamesD0forD0ptbin(Int_t nVars, TString *varNames, Bool_t *isUpperCut)
Int_t ApplyCutOnVariableMVA(Int_t nCutIndex, Int_t ptbin, Float_t cutVariableValue)
Bool_t fHardSelectionArrayITSD0SecondDaughter[7]
void SetNVars(Int_t nVars)
Definition: AliRDHFCuts.h:424
Double_t CosPointingAngleXY() const
void SetIsCutUsedD0forD0ptbin(Int_t nCutIndex, Int_t ptbin, Bool_t isCutUsed)
Int_t IsD0forD0ptbinSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod, Bool_t *bCutArray)
AliRDHFCuts & operator=(const AliRDHFCuts &source)
Bool_t fUsePID
Definition: AliRDHFCuts.h:455
virtual void GetCutVarsForOpt(AliAODRecoDecayHF *d, Float_t *vars, Int_t nvars, Int_t *pdgdaughters)
void AddTrackCutsSoftPi(const AliESDtrackCuts *cuts)
virtual AliESDtrackCuts * GetTrackCutsSoftPi() const
const Int_t nPtBins
Int_t ApplyCutOnVariableD0forD0ptbin(Int_t nCutIndex, Int_t ptbin, Float_t cutVariableValue, Bool_t bCutArray[29])
Int_t CheckBands(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, AliAODTrack *track)
void SetNPtBinsD0forD0ptbin(Int_t nptBins)
Int_t PtBinD0forD0ptbin(Double_t pt) const
Bool_t GetIsUpperCutForCutOptimization(Int_t nVariable) const
int Int_t
Definition: External.C:63
Bool_t * GetIsUpperCut() const
Definition: AliRDHFCuts.h:275
Int_t IsD0forD0ptbinSelectedMVA(TObject *obj, Int_t selectionLevel, AliAODEvent *aod, AliAODVertex *primaryVertex, Double_t bz)
unsigned int UInt_t
Definition: External.C:33
virtual Int_t IsSelectedPID(AliAODRecoDecayHF *rd)
void SetSoftSelectionArrayITSD0SecondDaughter(const Bool_t array[7]=0)
float Float_t
Definition: External.C:68
void SetCut(Int_t nCutIndex, Int_t ptBin, AliRDHFCutsBPlustoD0Pi::EUpperCut cutDirection, Float_t cutValue)
Bool_t CheckStatus(AliAODTrack *track, TString detectors) const
Float_t fMaxPtPid
cuts for pion from Bplus (AOD converted to ESD on the flight!)
void SetIsCutUsed(Int_t nCutIndex, Int_t ptbin, Bool_t isCutUsed)
void InitializeCutsForCutOptimization(Int_t nCutsForOptimization, Int_t nVariables)
AliRDHFCutsBPlustoD0Pi(const char *name="BPlustoD0PiCuts")
virtual Int_t SelectPID(AliAODTrack *track, Int_t type)
void SetSoftSelectionArrayITSBPlusPion(const Bool_t array[7]=0)
Int_t fIsSelectedPID
outcome of cuts selection
Definition: AliRDHFCuts.h:477
Int_t ApplyCutOnVariableD0forD0ptbinMVA(Int_t nCutIndex, Int_t ptbin, Float_t cutVariableValue)
Float_t * fCutsRD
fnVars*fnPtBins
Definition: AliRDHFCuts.h:453
Int_t IsBplusPionSelectedMVA(TObject *obj, Int_t selectionLevel, AliAODEvent *aod, AliAODVertex *primaryVertex, Double_t bz)
void SetHardSelectionArrayITSD0FirstDaughter(const Bool_t array[7]=0)
void SetVarsForOpt(Int_t nVars, Bool_t *forOpt)
void SetVarNames(Int_t nVars, TString *varNames, Bool_t *isUpperCut)
void SetCutForCutOptimization(Int_t nCutIndex, Int_t nVariable, Int_t ptBin, AliRDHFCutsBPlustoD0Pi::EUpperCut cutDirection, Float_t *cutValues)
Bool_t AreDaughtersSelected(AliAODRecoDecayHF *rd, const AliAODEvent *aod=0x0) const
void SetIsUpperCutForCutOptimization(Int_t nVariable, Bool_t isUpperCut)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod)
Bool_t IsPionRaw(AliAODTrack *track, TString detector) const
Int_t IsD0FromBPlusSelectedMVA(Double_t ptBPlus, TObject *obj, Int_t selectionLevel, AliAODEvent *aod, AliAODVertex *primaryVertex, Double_t bz)
Int_t GetGlobalIndexForCutOptimization(Int_t iCut, Int_t iVar, Int_t iPtBin)
Double_t GetFiducialYCut() const
void SetPtBins(Int_t nPtBinLimits, Float_t *ptBinLimits)
void SetHardSelectionArrayITSBPlusPion(const Bool_t array[7]=0)
void SetSoftSelectionArrayITSD0FirstDaughter(const Bool_t array[7]=0)
Int_t GetMatch() const
Definition: AliAODPidHF.h:145
Int_t IsD0SelectedPreRecVtxMVA(AliAODRecoDecayHF2Prong *d, AliAODTrack *pion, AliAODVertex *primaryVertex, Double_t bz, Int_t selLevel)
Double_t GetSigmaForCutOptimization(Int_t iPtBin) const
void SetCutsForCutOptimization(Int_t glIndex, Float_t *cutsRDForCutOptimization)
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:263
void SetNVarsD0forD0ptbin(Int_t nVars)
void SetPtBinsD0forD0ptbin(Int_t nPtBinLimits, Float_t *ptBinLimits)
AliRDHFCutsBPlustoD0Pi & operator=(const AliRDHFCutsBPlustoD0Pi &source)
Int_t GetCutIndexForCutOptimization(Int_t nVariable) const
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
Int_t fnPtBins
cuts on the candidate
Definition: AliRDHFCuts.h:445
AliAODPidHF * fPidHF
enable AOD049 centrality cleanup
Definition: AliRDHFCuts.h:457
Int_t GetNPtBinsD0forD0ptbin() const
Int_t fGlobalIndex
Definition: AliRDHFCuts.h:452
Int_t PtBin(Double_t pt) const
Bool_t fSoftSelectionArrayITSD0SecondDaughter[7]
void SetSigmaForCutOptimization(Double_t value, Int_t iPtBin)
void SetCutD0forD0ptbin(Int_t nCutIndex, Int_t ptBin, AliRDHFCutsBPlustoD0Pi::EUpperCut cutDirection, Float_t cutValue)
Int_t GetGlobalIndex(Int_t iVar, Int_t iPtBin) const
Int_t ApplyCutOnVariable(Int_t nCutIndex, Int_t ptbin, Float_t cutVariableValue, Bool_t bCutArray[68])
void SetCutsD0forD0ptbin(Int_t nVars, Int_t nPtBins, Float_t **cutsRDD0forD0ptbin)
Bool_t fGetCutInfo
Radius for circular PID nsigma cut.
Bool_t IsDaughterSelected(AliAODTrack *track, const AliESDVertex *primary, AliESDtrackCuts *cuts, const AliAODEvent *aod=0x0) const
void SetCutIndexForCutOptimization(Int_t nVariable, Int_t nCutIndex)
Int_t fnVars
Definition: AliRDHFCuts.h:448