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