AliPhysics  63e47e1 (63e47e1)
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  fMaxPtPid(9999.),
55  fTPCflag(999.),
56  fCircRadius(0.),
57  fGetCutInfo(0),
58  fIsCutUsed(0x0),
59  fnVarsD0forD0ptbin(0),
60  fnPtBinsD0forD0ptbin(1),
61  fGlobalIndexD0forD0ptbin(0),
62  fCutsRDD0forD0ptbin(0x0),
63  fnPtBinLimitsD0forD0ptbin(1),
64  fPtBinLimitsD0forD0ptbin(0x0),
65  fIsUpperCutD0forD0ptbin(0x0),
66  fIsCutUsedD0forD0ptbin(0x0),
67  fVarNamesD0forD0ptbin(0x0),
68  fMinITSNclsD0FirstDaughter(0),
69  fMinTPCNclsD0FirstDaughter(0),
70  fUseITSRefitD0FirstDaughter(0),
71  fUseTPCRefitD0FirstDaughter(0),
72  fUseFilterBitD0FirstDaughter(0),
73  fFilterBitD0FirstDaughter(0),
74  fMinPtD0FirstDaughter(0),
75  fMinITSNclsD0SecondDaughter(0),
76  fMinTPCNclsD0SecondDaughter(0),
77  fUseITSRefitD0SecondDaughter(0),
78  fUseTPCRefitD0SecondDaughter(0),
79  fUseFilterBitD0SecondDaughter(0),
80  fFilterBitD0SecondDaughter(0),
81  fMinPtD0SecondDaughter(0),
82  fMinITSNclsBPlusPion(0),
83  fMinTPCNclsBPlusPion(0),
84  fUseITSRefitBPlusPion(0),
85  fUseTPCRefitBPlusPion(0),
86  fUseFilterBitBPlusPion(0),
87  fFilterBitBPlusPion(0),
88  fMinPtBPlusPion(0),
89  fMaxAbsEtaD0FirstDaughter(999.),
90  fHardSelectionArrayITSD0FirstDaughter(),
91  fSoftSelectionArrayITSD0FirstDaughter(),
92  fNSoftITSCutD0FirstDaughter(0),
93  fMaxAbsEtaD0SecondDaughter(999.),
94  fHardSelectionArrayITSD0SecondDaughter(),
95  fSoftSelectionArrayITSD0SecondDaughter(),
96  fNSoftITSCutD0SecondDaughter(0),
97  fMaxAbsEtaBPlusPion(999.),
98  fHardSelectionArrayITSBPlusPion(),
99  fSoftSelectionArrayITSBPlusPion(),
100  fNSoftITSCutBPlusPion(0),
101  fMind0D0FirstDaughter(0),
102  fMind0D0SecondDaughter(0),
103  fMind0BPlusPion(0.)
104 {
105  //
106  // Default Constructor
107  //
108 
109  // Main cut setup as function of BPlus pt bins
110  const Int_t nvars=68;
111  SetNVars(nvars);
112 
113  TString varNames[nvars];
114  Int_t iterator = 0;
115 
116  //D0 cut variables
117  varNames[iterator++]= /*-00-*/ "inv. mass width[GeV]";
118  varNames[iterator++]= /*-01-*/ "delta mass width [GeV]"; //not used for D0
119  varNames[iterator++]= /*-02-*/ "pointing angle [Cos(theta)]";
120  varNames[iterator++]= /*-03-*/ "dca [cm]";
121  varNames[iterator++]= /*-04-*/ "Pt D0 [GeV/c]";
122  varNames[iterator++]= /*-05-*/ "Pt first daughter [GeV/c]";
123  varNames[iterator++]= /*-06-*/ "Pt second daughter [GeV/c]";
124  varNames[iterator++]= /*-07-*/ "d0 D0 [cm]";
125  varNames[iterator++]= /*-08-*/ "d0 first daughter [cm]";
126  varNames[iterator++]= /*-09-*/ "d0 second daughter [cm]";
127  varNames[iterator++]= /*-10-*/ "d0d0 [cm^2]";
128  varNames[iterator++]= /*-11-*/ "d0d0 XY [cm^2]";
129 
130  varNames[iterator++]= /*-12-*/ "angle between both daughters";
131  varNames[iterator++]= /*-13-*/ "angle mother with first daughter";
132  varNames[iterator++]= /*-14-*/ "angle mother with second daughter";
133  varNames[iterator++]= /*-15-*/ "cosThetaStar";
134  varNames[iterator++]= /*-16-*/ "vertexDistance";
135  varNames[iterator++]= /*-17-*/ "pseudoProperDecayTime";
136  varNames[iterator++]= /*-18-*/ "DecayTime";
137  varNames[iterator++]= /*-19-*/ "normalizedDecayTime";
138  varNames[iterator++]= /*-20-*/ "normDecayLength";
139  varNames[iterator++]= /*-21-*/ "topomatic first daughter";
140  varNames[iterator++]= /*-22-*/ "topomatic second daughter";
141  varNames[iterator++]= /*-23-*/ "topomatic max";
142  varNames[iterator++]= /*-24-*/ "topomatic min";
143 
144  varNames[iterator++]= /*-25-*/ "pointing angle XY [Cos(theta)]";
145  varNames[iterator++]= /*-26-*/ "vertex distance XY [cm]";
146  varNames[iterator++]= /*-27-*/ "normDecayLength XY";
147  varNames[iterator++]= /*-28-*/ "Chi2 per NDF vertex";
148 
149  varNames[iterator++]= /*-29-*/ "pointingAngleToBPlus";
150  varNames[iterator++]= /*-30-*/ "d0MotherToBPlus";
151  varNames[iterator++]= /*-31-*/ "d0FirstDaughterToBPlus";
152  varNames[iterator++]= /*-32-*/ "d0SecondDaughterToBPlus";
153  varNames[iterator++]= /*-33-*/ "impactProductToBPlus";
154  varNames[iterator++]= /*-34-*/ "impactProductXYToBPlus";
155  varNames[iterator++]= /*-35-*/ "normDecayLengthToBPlus";
156  varNames[iterator++]= /*-36-*/ "pseudoProperDecayTimeToBPlus";
157  varNames[iterator++]= /*-37-*/ "DecayTimeToBPlus";
158  varNames[iterator++]= /*-38-*/ "normalizedDecayTimeToBPlus";
159 
160  //BPlus cut variables
161  varNames[iterator++]= /*-39-*/ "inv. mass width[GeV]";
162  varNames[iterator++]= /*-40-*/ "delta mass width [GeV]";
163  varNames[iterator++]= /*-41-*/ "pointing angle [Cos(theta)]";
164  varNames[iterator++]= /*-42-*/ "dca [cm]";
165  varNames[iterator++]= /*-43-*/ "Pt BPlus [GeV/c]";
166  varNames[iterator++]= /*-44-*/ "Pt D0 [GeV/c]";
167  varNames[iterator++]= /*-45-*/ "Pt Pion [GeV/c]";
168  varNames[iterator++]= /*-46-*/ "d0 BPlus [cm]";
169  varNames[iterator++]= /*-47-*/ "d0 D0 [cm]";
170  varNames[iterator++]= /*-48-*/ "d0 Pion [cm]";
171  varNames[iterator++]= /*-49-*/ "d0d0 [cm^2]";
172  varNames[iterator++]= /*-50-*/ "d0d0 XY [cm^2]";
173 
174  varNames[iterator++]= /*-51-*/ "angle between both daughters";
175  varNames[iterator++]= /*-52-*/ "angle mother with first daughter";
176  varNames[iterator++]= /*-53-*/ "angle mother with second daughter";
177  varNames[iterator++]= /*-54-*/ "cosThetaStar";
178  varNames[iterator++]= /*-55-*/ "vertexDistance";
179  varNames[iterator++]= /*-56-*/ "pseudoProperDecayTime";
180  varNames[iterator++]= /*-57-*/ "DecayTime";
181  varNames[iterator++]= /*-58-*/ "normalizedDecayTime";
182  varNames[iterator++]= /*-59-*/ "normDecayLength";
183  varNames[iterator++]= /*-60-*/ "topomatic first daughter";
184  varNames[iterator++]= /*-61-*/ "topomatic second daughter";
185  varNames[iterator++]= /*-62-*/ "topomatic max";
186  varNames[iterator++]= /*-63-*/ "topomatic min";
187 
188  varNames[iterator++]= /*-64-*/ "pointing angle XY [Cos(theta)]";
189  varNames[iterator++]= /*-65-*/ "vertex distance XY [cm]";
190  varNames[iterator++]= /*-66-*/ "normDecayLength XY";
191  varNames[iterator++]= /*-67-*/ "Chi2 per NDF vertex";
192 
193  Bool_t isUpperCut[nvars]={0};
194 
195  SetVarNames(nvars,varNames,isUpperCut);
196 
197  Float_t limits[2]={0,999999999.};
198  SetPtBins(2,limits);
199 
200 
201  //
202  // Initialization of D0 cuts for D0 pt bins
203  //
204 
205  const Int_t nvarsD0forD0ptbin=29;
206  SetNVarsD0forD0ptbin(nvarsD0forD0ptbin);
207 
208  TString varNamesD0forD0ptbin[nvarsD0forD0ptbin];
209  iterator = 0;
210 
211  //D0 cut variables
212  varNamesD0forD0ptbin[iterator++]= /*-00-*/ "inv. mass width[GeV]";
213  varNamesD0forD0ptbin[iterator++]= /*-01-*/ "delta mass width [GeV]"; //not used for D0
214  varNamesD0forD0ptbin[iterator++]= /*-02-*/ "pointing angle [Cos(theta)]";
215  varNamesD0forD0ptbin[iterator++]= /*-03-*/ "dca [cm]";
216  varNamesD0forD0ptbin[iterator++]= /*-04-*/ "Pt D0 [GeV/c]";
217  varNamesD0forD0ptbin[iterator++]= /*-05-*/ "Pt first daughter [GeV/c]";
218  varNamesD0forD0ptbin[iterator++]= /*-06-*/ "Pt second daughter [GeV/c]";
219  varNamesD0forD0ptbin[iterator++]= /*-07-*/ "d0 D0 [cm]";
220  varNamesD0forD0ptbin[iterator++]= /*-08-*/ "d0 first daughter [cm]";
221  varNamesD0forD0ptbin[iterator++]= /*-09-*/ "d0 second daughter [cm]";
222  varNamesD0forD0ptbin[iterator++]= /*-10-*/ "d0d0 [cm^2]";
223  varNamesD0forD0ptbin[iterator++]= /*-11-*/ "d0d0 XY [cm^2]";
224 
225  varNamesD0forD0ptbin[iterator++]= /*-12-*/ "angle between both daughters";
226  varNamesD0forD0ptbin[iterator++]= /*-13-*/ "angle mother with first daughter";
227  varNamesD0forD0ptbin[iterator++]= /*-14-*/ "angle mother with second daughter";
228  varNamesD0forD0ptbin[iterator++]= /*-15-*/ "cosThetaStar";
229  varNamesD0forD0ptbin[iterator++]= /*-16-*/ "vertexDistance";
230  varNamesD0forD0ptbin[iterator++]= /*-17-*/ "pseudoProperDecayTime";
231  varNamesD0forD0ptbin[iterator++]= /*-18-*/ "DecayTime";
232  varNamesD0forD0ptbin[iterator++]= /*-19-*/ "normalizedDecayTime";
233  varNamesD0forD0ptbin[iterator++]= /*-20-*/ "normDecayLength";
234  varNamesD0forD0ptbin[iterator++]= /*-21-*/ "topomatic first daughter";
235  varNamesD0forD0ptbin[iterator++]= /*-22-*/ "topomatic second daughter";
236  varNamesD0forD0ptbin[iterator++]= /*-23-*/ "topomatic max";
237  varNamesD0forD0ptbin[iterator++]= /*-24-*/ "topomatic min";
238 
239  varNamesD0forD0ptbin[iterator++]= /*-25-*/ "pointing angle XY [Cos(theta)]";
240  varNamesD0forD0ptbin[iterator++]= /*-26-*/ "vertex distance XY [cm]";
241  varNamesD0forD0ptbin[iterator++]= /*-27-*/ "normDecayLength XY";
242  varNamesD0forD0ptbin[iterator++]= /*-28-*/ "Chi2 per NDF vertex";
243 
244  Bool_t isUpperCutD0forD0ptbin[nvarsD0forD0ptbin]={0};
245 
246  SetVarNamesD0forD0ptbin(nvarsD0forD0ptbin,varNamesD0forD0ptbin,isUpperCutD0forD0ptbin);
247 
248  Float_t limitsD0forD0ptbin[2]={0,999999999.};
249  SetPtBinsD0forD0ptbin(2,limitsD0forD0ptbin);
250 
251  Bool_t forOpt[16]={0}; //not yet used for BPlus analysis
252  SetVarsForOpt(16,forOpt);
253 
254 }
255 //--------------------------------------------------------------------------
257  AliRDHFCuts(source),
258  fMaxPtPid(source.fMaxPtPid),
259  fTPCflag(source.fTPCflag),
260  fCircRadius(source.fCircRadius),
261  fGetCutInfo(source.fGetCutInfo),
262  fIsCutUsed(0x0),
266  fCutsRDD0forD0ptbin(0x0),
308 {
309  //
310  // Copy constructor
311  //
312 
315  if(source.fIsCutUsed)
316  {
317  if(fIsCutUsed) {
318  delete [] fIsCutUsed;
319  fIsCutUsed = NULL;
320  }
321  fIsCutUsed = new Bool_t[(source.GetNPtBins())*(source.GetNVars())];
322 
323  for (Int_t i = 0; i < source.fnVars; ++i)
324  {
325  for(Int_t j = 0; j < source.fnPtBins; j++)
326  {
327  Bool_t bUse = source.GetIsCutUsed(i,j);
328  SetIsCutUsed(i,j,bUse);
329  }
330  }
331  }
332  if(source.fIsCutUsedD0forD0ptbin)
333  {
335  delete [] fIsCutUsedD0forD0ptbin;
336  fIsCutUsedD0forD0ptbin = NULL;
337  }
339  for (Int_t i = 0; i < source.fnVarsD0forD0ptbin; ++i)
340  {
341  for(Int_t j = 0; j < source.fnPtBinsD0forD0ptbin; j++)
342  {
343  Bool_t bUse = source.GetIsCutUsedD0forD0ptbin(i,j);
344  SetIsCutUsedD0forD0ptbin(i,j,bUse);
345  }
346  }
347  }
349 
356 }
357 //--------------------------------------------------------------------------
359  //
360  // Default Destructor
361  //
362  if(fIsCutUsed) { delete [] fIsCutUsed; fIsCutUsed=0; }
368 }
369 //--------------------------------------------------------------------------
371 {
372  //
373  // assignment operator
374  //
375 
376  if(&source == this) return *this;
377 
378  AliRDHFCuts::operator=(source);
379 
380  fMaxPtPid = source.fMaxPtPid;
381  fTPCflag = source.fTPCflag;
382  fCircRadius = source.fCircRadius;
383  fGetCutInfo = source.fGetCutInfo;
418 
421  if(source.fIsCutUsed)
422  {
423  if(fIsCutUsed) {
424  delete [] fIsCutUsed;
425  fIsCutUsed = NULL;
426  }
427  fIsCutUsed = new Bool_t[(source.GetNPtBins())*(source.GetNVars())];
428 
429  for (Int_t i = 0; i < source.fnVars; ++i)
430  {
431  for(Int_t j = 0; j < source.fnPtBins; j++)
432  {
433  Bool_t bUse = source.GetIsCutUsed(i,j);
434  SetIsCutUsed(i,j,bUse);
435  }
436  }
437  }
438  if(source.fIsCutUsedD0forD0ptbin)
439  {
441  delete [] fIsCutUsedD0forD0ptbin;
442  fIsCutUsedD0forD0ptbin = NULL;
443  }
445  for (Int_t i = 0; i < source.fnVarsD0forD0ptbin; ++i)
446  {
447  for(Int_t j = 0; j < source.fnPtBinsD0forD0ptbin; j++)
448  {
449  Bool_t bUse = source.GetIsCutUsedD0forD0ptbin(i,j);
450  SetIsCutUsedD0forD0ptbin(i,j,bUse);
451  }
452  }
453  }
455 
462 
463  return *this;
464 }
465 //--------------------------------------------------------------------------
467  // not yet used
468 
469  return;
470 }
471 //--------------------------------------------------------------------------
472 Int_t AliRDHFCutsBPlustoD0Pi::IsSelected(TObject* obj,Int_t selectionLevel, AliAODEvent* aod, Bool_t bCutArray[68]) {
473  //
474  // In this function we apply the selection cuts on the BPlus candidate and its daughters.
475  // The function returns 0 if the candidate is cut and is able to return information on which cuts the candidate passes.
476  //
477 
478  fIsSelectedCuts=0;
479  fIsSelectedPID=0;
480 
481  // The cuts used in this class have to be set in the maketfile.
482  if(!fCutsRD){
483  cout<<"Cut matrice not inizialized. Exit..."<<endl;
484  return 0;
485  }
486 
487  AliAODRecoDecayHF2Prong* candidateBPlus = (AliAODRecoDecayHF2Prong*)obj;
488  if(!candidateBPlus){
489  cout<<"candidateBPlus null"<<endl;
490  return 0;
491  }
492 
493  AliAODRecoDecayHF2Prong* candidateD0 = (AliAODRecoDecayHF2Prong*)candidateBPlus->GetDaughter(1);
494  if(!candidateD0){
495  cout<<"candidateD0 null"<<endl;
496  return 0;
497  }
498 
499  AliAODTrack *candidatePion = (AliAODTrack*)candidateBPlus->GetDaughter(0);
500  if(!candidatePion){
501  cout<<"candidatePion null 1"<<endl;
502  return 0;
503  }
504 
505  AliAODVertex * vertexBPlus = candidateBPlus->GetSecondaryVtx();
506  if(!vertexBPlus){
507  cout<<"vertexBPlus null"<<endl;
508  return 0;
509  }
510 
511  AliAODVertex * primaryVertex = aod->GetPrimaryVertex();
512  if(!primaryVertex){
513  cout<<"primaryVertex null"<<endl;
514  return 0;
515  }
516 
517  Int_t returnvalue=1;
518  Bool_t bPassedCut = kFALSE;
519 
520  //get the magnetic field
521  Double_t bz = (Double_t)aod->GetMagneticField();
522 
523  // selection on candidate
524  if(selectionLevel==AliRDHFCuts::kAll ||
525  selectionLevel==AliRDHFCuts::kCandidate) {
526 
527  // We check to which pt bin the candidate belongs
528  Int_t ptbin = PtBin(candidateBPlus->Pt());
529  if(ptbin==-1) return -1;
530 
531  // We obtain the variable values in the section below
532  // D0Mass and BPlusmass
533  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
534  Double_t mBPlusPDG = TDatabasePDG::Instance()->GetParticle(521)->Mass();
535 
536  // delta mass PDG
537  Double_t deltaPDG = mBPlusPDG - mD0PDG;
538 
539  // Half width BPlus mass
540  UInt_t prongs[2];
541  prongs[0] = 211; prongs[1] = 421;
542  Double_t invMassBPlus = candidateBPlus->InvMass(2,prongs);
543  Double_t invMassDifference = TMath::Abs(mBPlusPDG - invMassBPlus);
544  Double_t invMassDelta = TMath::Abs(deltaPDG-(DeltaInvMassBPlusKpipi(candidateBPlus)));
545 
546  Double_t pointingAngle = candidateBPlus->CosPointingAngle();
547  Double_t dcaMother = candidateBPlus->GetDCA();
548  Double_t ptMother = candidateBPlus->Pt();
549  Double_t momentumMother = candidateBPlus->P();
550  Double_t ptD0 = candidateD0->Pt();
551  Double_t ptPion = candidatePion->Pt();
552 
553  AliExternalTrackParam motherTrack;
554  motherTrack.CopyFromVTrack(candidateBPlus);
555  Double_t d0z0[2],covd0z0[3],d0[2];
556  motherTrack.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
557  d0[0] = d0z0[0];
558  Double_t d0Mother = TMath::Abs(d0[0]);
559  Double_t d0firstTrack = TMath::Abs(candidateBPlus->Getd0Prong(0));
560  Double_t d0secondTrack = TMath::Abs(candidateBPlus->Getd0Prong(1));
561 
562  Double_t impactProduct = candidateBPlus->Prodd0d0();
563  Double_t impactProductXY = TMath::Abs(candidateBPlus->ImpParXY());
564 
565  Double_t angleBetweenBothDaughters = (candidateD0->Px() * candidatePion->Px() + candidateD0->Py() * candidatePion->Py() + candidateD0->Pz() * candidatePion->Pz()) /(candidateD0->P() * candidatePion->P());
566  Double_t angleMotherFirstDaughter = (candidateBPlus->Px() * candidatePion->Px() + candidateBPlus->Py() * candidatePion->Py() + candidateBPlus->Pz() * candidatePion->Pz()) /(candidateBPlus->P() * candidatePion->P());
567  Double_t angleMotherSecondDaughter = (candidateBPlus->Px() * candidateD0->Px() + candidateBPlus->Py() * candidateD0->Py() + candidateBPlus->Pz() * candidateD0->Pz()) /(candidateBPlus->P() * candidateD0->P());
568 
569  Double_t cosThetaStar = candidateBPlus->CosThetaStar(0,521,211,421);
570  Double_t vertexDistance = vertexBPlus->DistanceToVertex(primaryVertex);
571  Double_t normDecayLength = candidateBPlus->NormalizedDecayLength();
572  Double_t pdgMassMother = TDatabasePDG::Instance()->GetParticle(521)->Mass();
573  Double_t pseudoProperDecayLength = ((vertexBPlus->GetX() - primaryVertex->GetX()) * candidateBPlus->Px() / TMath::Abs(candidateBPlus->Pt())) + ((vertexBPlus->GetY() - primaryVertex->GetY()) * candidateBPlus->Py() / TMath::Abs(candidateBPlus->Pt()));
574  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother/ptMother;
575  Double_t decayTime = vertexDistance / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother/(momentumMother*momentumMother)) + 1)));
576 
577  Double_t phi = candidateBPlus->Phi();
578  Double_t theta = candidateBPlus->Theta();
579  Double_t covMatrix[21];
580  candidateBPlus->GetCovarianceXYZPxPyPz(covMatrix);
581 
582  Double_t cp = TMath::Cos(phi);
583  Double_t sp = TMath::Sin(phi);
584  Double_t ct = TMath::Cos(theta);
585  Double_t st = TMath::Sin(theta);
586 
587  Double_t errorMomentum = covMatrix[9]*cp*cp*ct*ct // GetCovPxPx
588  +covMatrix[13]*2.*cp*sp*ct*ct // GetCovPxPy
589  +covMatrix[18]*2.*cp*ct*st // GetCovPxPz
590  +covMatrix[14]*sp*sp*ct*ct // GetCovPyPy
591  +covMatrix[19]*2.*sp*ct*st // GetCovPyPz
592  +covMatrix[20]*st*st; // GetCovPzPz
593  Double_t normalizedDecayTime = candidateBPlus->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother*errorMomentum*errorMomentum/(momentumMother*momentumMother)) + 1)));
594 
595  Double_t cosPointingAngleXY = candidateBPlus->CosPointingAngleXY();
596  Double_t distanceXYToVertex = vertexBPlus->DistanceXYToVertex(primaryVertex);
597  Double_t normalizedDecayLengthXY = candidateBPlus->NormalizedDecayLengthXY();
598  Double_t chi2Vertex = vertexBPlus->GetChi2perNDF();
599 
600  //Topomatic
601  Double_t dd0pr1=0.;
602  Double_t dd0pr2=0.;
603  Double_t dd0max=0.;
604  Double_t dd0min=0.;
605  for(Int_t ipr=0; ipr<2; ipr++)
606  {
607  Double_t diffIP, errdiffIP;
608  candidateBPlus->Getd0MeasMinusExpProng(ipr,bz,diffIP,errdiffIP);
609  Double_t normdd0=0.;
610  if(errdiffIP>0.) normdd0=diffIP/errdiffIP;
611  if(ipr==0) dd0pr1=normdd0;
612  if(ipr==1) dd0pr2=normdd0;
613  }
614  if(TMath::Abs(dd0pr1)>TMath::Abs(dd0pr2)) {dd0max=dd0pr1; dd0min=dd0pr2;}
615  else {dd0max=dd0pr2; dd0min=dd0pr1;}
616 
617 
618  // We apply the cuts
619  Int_t nCutIndex = 0;
620  Double_t cutVariableValue = 0.0;
621 
622  // "inv. mass width [GeV]" --------------------------------------------
623  nCutIndex = 39;
624  cutVariableValue = invMassDifference;
625  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
626  if(!bPassedCut && !fGetCutInfo) return 0;
627  //---------------------------------------------------------------------
628 
629  // "delta mass width [GeV]" -------------------------------------------
630  nCutIndex = 40;
631  cutVariableValue = invMassDelta;
632  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
633  if(!bPassedCut && !fGetCutInfo) return 0;
634  //---------------------------------------------------------------------
635 
636  // "pointing angle [Cos(theta)]" --------------------------------------
637  nCutIndex = 41;
638  cutVariableValue = pointingAngle;
639  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
640  if(!bPassedCut && !fGetCutInfo) return 0;
641  //---------------------------------------------------------------------
642 
643  // "dca [cm]" ---------------------------------------------------------
644  nCutIndex = 42;
645  cutVariableValue = dcaMother;
646  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
647  if(!bPassedCut && !fGetCutInfo) return 0;
648  //---------------------------------------------------------------------
649 
650  // "Pt BPlus [GeV/c]" ----------------------------------------------------
651  nCutIndex = 43;
652  cutVariableValue = ptMother;
653  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
654  if(!bPassedCut && !fGetCutInfo) return 0;
655  //---------------------------------------------------------------------
656 
657  // "Pt D0 [GeV/c]" -------------------------------------------------
658  nCutIndex = 44;
659  cutVariableValue = ptD0;
660  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
661  if(!bPassedCut && !fGetCutInfo) return 0;
662  //---------------------------------------------------------------------
663 
664  // "Pt Pion [GeV/c]" --------------------------------------------------
665  nCutIndex = 45;
666  cutVariableValue = ptPion;
667  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
668  if(!bPassedCut && !fGetCutInfo) return 0;
669  //---------------------------------------------------------------------
670 
671  // "d0 BPlus [cm]" -------------------------------------------------------
672  nCutIndex = 46;
673  cutVariableValue = d0Mother;
674  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
675  if(!bPassedCut && !fGetCutInfo) return 0;
676  //---------------------------------------------------------------------
677 
678  // "d0 D0 [cm]"-----------------------------------------------------
679  nCutIndex = 47;
680  cutVariableValue = d0firstTrack;
681  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
682  if(!bPassedCut && !fGetCutInfo) return 0;
683  //---------------------------------------------------------------------
684 
685  // "d0 Pion [cm]" -----------------------------------------------------
686  nCutIndex = 48;
687  cutVariableValue = d0secondTrack;
688  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
689  if(!bPassedCut && !fGetCutInfo) return 0;
690  //---------------------------------------------------------------------
691 
692  // "d0d0 [cm^2]" ------------------------------------------------------
693  nCutIndex = 49;
694  cutVariableValue = impactProduct;
695  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
696  if(!bPassedCut && !fGetCutInfo) return 0;
697  //---------------------------------------------------------------------
698 
699  // "d0d0 XY [cm^2]" ---------------------------------------------------
700  nCutIndex = 50;
701  cutVariableValue = impactProductXY;
702  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
703  if(!bPassedCut && !fGetCutInfo) return 0;
704  //---------------------------------------------------------------------
705 
706  // "angle between both daughters" -------------------------------------
707  nCutIndex = 51;
708  cutVariableValue = angleBetweenBothDaughters;
709  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
710  if(!bPassedCut && !fGetCutInfo) return 0;
711  //---------------------------------------------------------------------
712 
713  // "angle mother with first daughter" ---------------------------------
714  nCutIndex = 52;
715  cutVariableValue = angleMotherFirstDaughter;
716  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
717  if(!bPassedCut && !fGetCutInfo) return 0;
718  //---------------------------------------------------------------------
719 
720  // "angle mother with second daughter" --------------------------------
721  nCutIndex = 53;
722  cutVariableValue = angleMotherSecondDaughter;
723  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
724  if(!bPassedCut && !fGetCutInfo) return 0;
725  //---------------------------------------------------------------------
726 
727  // "cosThetaStar" -----------------------------------------------------
728  nCutIndex = 54;
729  cutVariableValue = cosThetaStar;
730  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
731  if(!bPassedCut && !fGetCutInfo) return 0;
732  //---------------------------------------------------------------------
733 
734  // "vertexDistance" ---------------------------------------------------
735  nCutIndex = 55;
736  cutVariableValue = vertexDistance;
737  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
738  if(!bPassedCut && !fGetCutInfo) return 0;
739  //---------------------------------------------------------------------
740 
741  // "pseudoProperDecayTime" --------------------------------------------
742  nCutIndex = 56;
743  cutVariableValue = pseudoProperDecayTime;
744  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
745  if(!bPassedCut && !fGetCutInfo) return 0;
746  //---------------------------------------------------------------------
747 
748  // "DecayTime" --------------------------------------------------------
749  nCutIndex = 57;
750  cutVariableValue = decayTime;
751  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
752  if(!bPassedCut && !fGetCutInfo) return 0;
753  //---------------------------------------------------------------------
754 
755  // "normalizedDecayTime" ----------------------------------------------------
756  nCutIndex = 58;
757  cutVariableValue = normalizedDecayTime;
758  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
759  if(!bPassedCut && !fGetCutInfo) return 0;
760  //---------------------------------------------------------------------
761 
762  // "normDecayLength" --------------------------------------------------
763  nCutIndex = 59;
764  cutVariableValue = normDecayLength;
765  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
766  if(!bPassedCut && !fGetCutInfo) return 0;
767  //---------------------------------------------------------------------
768 
769  // "topomatic first daughter" -----------------------------------------
770  nCutIndex = 60;
771  cutVariableValue = dd0pr1;
772  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
773  if(!bPassedCut && !fGetCutInfo) return 0;
774  //---------------------------------------------------------------------
775 
776  // "topomatic second daughter" ----------------------------------------
777  nCutIndex = 61;
778  cutVariableValue = dd0pr2;
779  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
780  if(!bPassedCut && !fGetCutInfo) return 0;
781  //---------------------------------------------------------------------
782 
783  // "topomatic max" ----------------------------------------------------
784  nCutIndex = 62;
785  cutVariableValue = dd0max;
786  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
787  if(!bPassedCut && !fGetCutInfo) return 0;
788  //---------------------------------------------------------------------
789 
790  // "topomatic min" ----------------------------------------------------
791  nCutIndex = 63;
792  cutVariableValue = dd0min;
793  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
794  if(!bPassedCut && !fGetCutInfo) return 0;
795  //---------------------------------------------------------------------
796 
797  // "pointing angle XY" ----------------------------------------------------
798  nCutIndex = 64;
799  cutVariableValue = cosPointingAngleXY;
800  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
801  if(!bPassedCut && !fGetCutInfo) return 0;
802  //---------------------------------------------------------------------
803 
804  // "vertex distance XY" ----------------------------------------------------
805  nCutIndex = 65;
806  cutVariableValue = distanceXYToVertex;
807  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
808  if(!bPassedCut && !fGetCutInfo) return 0;
809  //---------------------------------------------------------------------
810 
811  // "normalized decay length XY" ----------------------------------------------------
812  nCutIndex = 66;
813  cutVariableValue = normalizedDecayLengthXY;
814  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
815  if(!bPassedCut && !fGetCutInfo) return 0;
816  //---------------------------------------------------------------------
817 
818  // "chi squared per NDF" ----------------------------------------------------
819  nCutIndex = 67;
820  cutVariableValue = chi2Vertex;
821  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
822  if(!bPassedCut && !fGetCutInfo) return 0;
823  //---------------------------------------------------------------------
824 
825  // select D0
826  bPassedCut = IsD0FromBPlusSelected(ptMother,candidateBPlus,selectionLevel,aod,bCutArray);
827  if(!bPassedCut && !fGetCutInfo) return 0;
828  }
829 
830  if(bPassedCut==kFALSE)
831  {
832  returnvalue = 0;
833  } else
834  {
835  for (Int_t i = 39; i < 68; ++i)
836  {
837  if(bCutArray[i]==kTRUE){
838  returnvalue = 0;
839  break;
840  }
841  }
842  }
843 
844  fIsSelectedCuts = returnvalue;
845 
846  return returnvalue;
847 }
848 //_________________________________________________________________________________________________
849 Int_t AliRDHFCutsBPlustoD0Pi::IsD0FromBPlusSelected(Double_t ptBPlus, TObject* obj,Int_t selectionLevel, AliAODEvent* aod, Bool_t bCutArray[68]) {
850  //
851  // 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.
852  //
853 
854  if(!fCutsRD){
855  cout<<"Cut matrice not inizialized. Exit..."<<endl;
856  return 0;
857  }
858 
859  AliAODRecoDecayHF2Prong* candidateBPlus = (AliAODRecoDecayHF2Prong*)obj;
860  if(!candidateBPlus){
861  cout<<"candidateBPlus null"<<endl;
862  return 0;
863  }
864 
865  AliAODRecoDecayHF2Prong* candidateD0 = (AliAODRecoDecayHF2Prong*)candidateBPlus->GetDaughter(1);
866  if(!candidateD0){
867  cout<<"candidateD0 null"<<endl;
868  return 0;
869  }
870 
871  AliAODTrack *candidatePion = (AliAODTrack*)candidateD0->GetDaughter(0);
872  if(!candidatePion){
873  cout<<"candidatePion null 2"<<endl;
874  return 0;
875  }
876 
877  AliAODTrack *candidateKaon = (AliAODTrack*)candidateD0->GetDaughter(1);
878  if(!candidateKaon){
879  cout<<"candidateKaon null"<<endl;
880  return 0;
881  }
882 
883  AliAODVertex * vertexBPlus = candidateBPlus->GetSecondaryVtx();
884  if(!vertexBPlus){
885  cout<<"vertexBPlus null"<<endl;
886  return 0;
887  }
888 
889  AliAODVertex * vertexD0 = candidateD0->GetSecondaryVtx();
890  if(!vertexD0){
891  cout<<"vertexD0 null"<<endl;
892  return 0;
893  }
894 
895  AliAODVertex * primaryVertex = aod->GetPrimaryVertex();
896  if(!primaryVertex){
897  cout<<"primaryVertex null"<<endl;
898  return 0;
899  }
900 
901  Int_t returnvalue=1;
902  Bool_t bPassedCut = kFALSE;
903 
904  //get the magnetic field
905  Double_t bz = (Double_t)aod->GetMagneticField();
906 
907 
908  // selection on candidate
909  if(selectionLevel==AliRDHFCuts::kAll ||
910  selectionLevel==AliRDHFCuts::kCandidate) {
911 
912  Int_t ptbin=PtBin(ptBPlus);
913 
914  // D0mass
915  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
916 
917  // D0 window - invariant mass
918  Int_t chargeBPlus = candidateBPlus->Charge();
919  UInt_t prongs[2];
920  if(chargeBPlus==1)
921  {
922  prongs[0] = 211;
923  prongs[1] = 321;
924  }
925  else if (chargeBPlus==-1)
926  {
927  prongs[1] = 211;
928  prongs[0] = 321;
929  }
930  else
931  {
932  cout << "Wrong charge BPlus." << endl;
933  return 0;
934  }
935  Double_t invMassD0 = candidateD0->InvMass(2,prongs);
936  Double_t invMassDifference = TMath::Abs(mD0PDG - invMassD0);
937 
938  Double_t pointingAngle = candidateD0->CosPointingAngle();
939  Double_t dcaMother = candidateD0->GetDCA();
940  Double_t ptMother = candidateD0->Pt();
941  Double_t momentumMother = candidateD0->P();
942  Double_t ptPion = candidatePion->Pt();
943  Double_t ptKaon = candidateKaon->Pt();
944 
945  AliExternalTrackParam motherTrack;
946  motherTrack.CopyFromVTrack(candidateD0);
947  Double_t d0z0[2],covd0z0[3],d0[2];
948  motherTrack.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
949  d0[0] = d0z0[0];
950  Double_t d0Mother = TMath::Abs(d0[0]);
951  Double_t d0firstTrack = TMath::Abs(candidateD0->Getd0Prong(0));
952  Double_t d0secondTrack = TMath::Abs(candidateD0->Getd0Prong(1));
953 
954  Double_t impactProduct = candidateD0->Prodd0d0();
955  Double_t impactProductXY = TMath::Abs(candidateD0->ImpParXY());
956 
957  Double_t angleBetweenBothDaughters = (candidateKaon->Px() * candidatePion->Px() + candidateKaon->Py() * candidatePion->Py() + candidateKaon->Pz() * candidatePion->Pz()) /(candidateKaon->P() * candidatePion->P());
958  Double_t angleMotherFirstDaughter = (candidateD0->Px() * candidatePion->Px() + candidateD0->Py() * candidatePion->Py() + candidateD0->Pz() * candidatePion->Pz()) /(candidateD0->P() * candidatePion->P());
959  Double_t angleMotherSecondDaughter = (candidateD0->Px() * candidateKaon->Px() + candidateD0->Py() * candidateKaon->Py() + candidateD0->Pz() * candidateKaon->Pz()) /(candidateD0->P() * candidateKaon->P());
960 
961  Double_t cosThetaStar = candidateD0->CosThetaStar(0,421,211,321);
962  Double_t vertexDistance = vertexD0->DistanceToVertex(primaryVertex);
963  Double_t normDecayLength = candidateD0->NormalizedDecayLength();
964  Double_t pdgMassMother = TDatabasePDG::Instance()->GetParticle(421)->Mass();
965  Double_t pseudoProperDecayLength = ((vertexD0->GetX() - primaryVertex->GetX()) * candidateD0->Px() / TMath::Abs(candidateD0->Pt())) + ((vertexD0->GetY() - primaryVertex->GetY()) * candidateD0->Py() / TMath::Abs(candidateD0->Pt()));
966  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother/ptMother;
967  Double_t decayTime = vertexDistance / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother/(momentumMother*momentumMother)) + 1)));
968 
969  Double_t phi = candidateD0->Phi();
970  Double_t theta = candidateD0->Theta();
971  Double_t covMatrix[21];
972  candidateD0->GetCovarianceXYZPxPyPz(covMatrix);
973 
974  Double_t cp = TMath::Cos(phi);
975  Double_t sp = TMath::Sin(phi);
976  Double_t ct = TMath::Cos(theta);
977  Double_t st = TMath::Sin(theta);
978 
979  Double_t errorMomentum = covMatrix[9]*cp*cp*ct*ct // GetCovPxPx
980  +covMatrix[13]*2.*cp*sp*ct*ct // GetCovPxPy
981  +covMatrix[18]*2.*cp*ct*st // GetCovPxPz
982  +covMatrix[14]*sp*sp*ct*ct // GetCovPyPy
983  +covMatrix[19]*2.*sp*ct*st // GetCovPyPz
984  +covMatrix[20]*st*st; // GetCovPzPz
985  Double_t normalizedDecayTime = candidateD0->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother*errorMomentum*errorMomentum/(momentumMother*momentumMother)) + 1)));
986 
987  Double_t cosPointingAngleXY = candidateD0->CosPointingAngleXY();
988  Double_t distanceXYToVertex = vertexD0->DistanceXYToVertex(primaryVertex);
989  Double_t normalizedDecayLengthXY = candidateD0->NormalizedDecayLengthXY();
990  Double_t chi2Vertex = vertexD0->GetChi2perNDF();
991 
992 
993  //Topomatic
994  Double_t dd0pr1=0.;
995  Double_t dd0pr2=0.;
996  Double_t dd0max=0.;
997  Double_t dd0min=0.;
998  for(Int_t ipr=0; ipr<2; ipr++)
999  {
1000  Double_t diffIP, errdiffIP;
1001  candidateD0->Getd0MeasMinusExpProng(ipr,bz,diffIP,errdiffIP);
1002  Double_t normdd0=0.;
1003  if(errdiffIP>0.) normdd0=diffIP/errdiffIP;
1004  if(ipr==0) dd0pr1=normdd0;
1005  if(ipr==1) dd0pr2=normdd0;
1006  }
1007  if(TMath::Abs(dd0pr1)>TMath::Abs(dd0pr2)) {dd0max=dd0pr1; dd0min=dd0pr2;}
1008  else {dd0max=dd0pr2; dd0min=dd0pr1;}
1009 
1010 
1011  // We apply the cuts
1012  Int_t nCutIndex = 0;
1013  Double_t cutVariableValue = 0.0;
1014 
1015  // "inv. mass width [GeV]" --------------------------------------------
1016  nCutIndex = 0;
1017  cutVariableValue = invMassDifference;
1018  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1019  if(!bPassedCut && !fGetCutInfo) return 0;
1020  //---------------------------------------------------------------------
1021 
1022  // "delta mass width [GeV]" -------------------------------------------
1023  // nCutIndex = 1; // not used for D0
1024  // cutVariableValue = invMassDelta;
1025  // bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1026  // if(!bPassedCut && !fGetCutInfo) return 0;
1027  //---------------------------------------------------------------------
1028 
1029  // "pointing angle [Cos(theta)]" --------------------------------------
1030  nCutIndex = 2;
1031  cutVariableValue = pointingAngle;
1032  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1033  if(!bPassedCut && !fGetCutInfo) return 0;
1034  //---------------------------------------------------------------------
1035 
1036  // "dca [cm]" ---------------------------------------------------------
1037  nCutIndex = 3;
1038  cutVariableValue = dcaMother;
1039  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1040  if(!bPassedCut && !fGetCutInfo) return 0;
1041  //---------------------------------------------------------------------
1042 
1043  // "Pt D0 [GeV/c]" ----------------------------------------------------
1044  nCutIndex = 4;
1045  cutVariableValue = ptMother;
1046  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1047  if(!bPassedCut && !fGetCutInfo) return 0;
1048  //---------------------------------------------------------------------
1049 
1050  // "Pt Kaon [GeV/c]" -------------------------------------------------
1051  nCutIndex = 5;
1052  cutVariableValue = ptKaon;
1053  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1054  if(!bPassedCut && !fGetCutInfo) return 0;
1055  //---------------------------------------------------------------------
1056 
1057  // "Pt Pion [GeV/c]" --------------------------------------------------
1058  nCutIndex = 6;
1059  cutVariableValue = ptPion;
1060  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1061  if(!bPassedCut && !fGetCutInfo) return 0;
1062  //---------------------------------------------------------------------
1063 
1064  // "d0 D0 [cm]" -------------------------------------------------------
1065  nCutIndex = 7;
1066  cutVariableValue = d0Mother;
1067  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1068  if(!bPassedCut && !fGetCutInfo) return 0;
1069  //---------------------------------------------------------------------
1070 
1071  // "d0 Kaon [cm]"-----------------------------------------------------
1072  nCutIndex = 8;
1073  cutVariableValue = d0firstTrack;
1074  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1075  if(!bPassedCut && !fGetCutInfo) return 0;
1076  //---------------------------------------------------------------------
1077 
1078  // "d0 Pion [cm]" -----------------------------------------------------
1079  nCutIndex = 9;
1080  cutVariableValue = d0secondTrack;
1081  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1082  if(!bPassedCut && !fGetCutInfo) return 0;
1083  //---------------------------------------------------------------------
1084 
1085  // "d0d0 [cm^2]" ------------------------------------------------------
1086  nCutIndex = 10;
1087  cutVariableValue = impactProduct;
1088  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1089  if(!bPassedCut && !fGetCutInfo) return 0;
1090  //---------------------------------------------------------------------
1091 
1092  // "d0d0 XY [cm^2]" ---------------------------------------------------
1093  nCutIndex = 11;
1094  cutVariableValue = impactProductXY;
1095  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1096  if(!bPassedCut && !fGetCutInfo) return 0;
1097  //---------------------------------------------------------------------
1098 
1099  // "angle between both daughters" -------------------------------------
1100  nCutIndex = 12;
1101  cutVariableValue = angleBetweenBothDaughters;
1102  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1103  if(!bPassedCut && !fGetCutInfo) return 0;
1104  //---------------------------------------------------------------------
1105 
1106  // "angle mother with first daughter" ---------------------------------
1107  nCutIndex = 13;
1108  cutVariableValue = angleMotherFirstDaughter;
1109  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1110  if(!bPassedCut && !fGetCutInfo) return 0;
1111  //---------------------------------------------------------------------
1112 
1113  // "angle mother with second daughter" --------------------------------
1114  nCutIndex = 14;
1115  cutVariableValue = angleMotherSecondDaughter;
1116  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1117  if(!bPassedCut && !fGetCutInfo) return 0;
1118  //---------------------------------------------------------------------
1119 
1120  // "cosThetaStar" -----------------------------------------------------
1121  nCutIndex = 15;
1122  cutVariableValue = cosThetaStar;
1123  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1124  if(!bPassedCut && !fGetCutInfo) return 0;
1125  //---------------------------------------------------------------------
1126 
1127  // "vertexDistance" ---------------------------------------------------
1128  nCutIndex = 16;
1129  cutVariableValue = vertexDistance;
1130  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1131  if(!bPassedCut && !fGetCutInfo) return 0;
1132  //---------------------------------------------------------------------
1133 
1134  // "pseudoProperDecayTime" --------------------------------------------
1135  nCutIndex = 17;
1136  cutVariableValue = pseudoProperDecayTime;
1137  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1138  if(!bPassedCut && !fGetCutInfo) return 0;
1139  //---------------------------------------------------------------------
1140 
1141  // "DecayTime" --------------------------------------------------------
1142  nCutIndex = 18;
1143  cutVariableValue = decayTime;
1144  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1145  if(!bPassedCut && !fGetCutInfo) return 0;
1146  //---------------------------------------------------------------------
1147 
1148  // "normalizedDecayTime" ----------------------------------------------------
1149  nCutIndex = 19;
1150  cutVariableValue = normalizedDecayTime;
1151  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1152  if(!bPassedCut && !fGetCutInfo) return 0;
1153  //---------------------------------------------------------------------
1154 
1155  // "normDecayLength" --------------------------------------------------
1156  nCutIndex = 20;
1157  cutVariableValue = normDecayLength;
1158  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1159  if(!bPassedCut && !fGetCutInfo) return 0;
1160  //---------------------------------------------------------------------
1161 
1162  // "topomatic first daughter" -----------------------------------------
1163  nCutIndex = 21;
1164  cutVariableValue = dd0pr1;
1165  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1166  if(!bPassedCut && !fGetCutInfo) return 0;
1167  //---------------------------------------------------------------------
1168 
1169  // "topomatic second daughter" ----------------------------------------
1170  nCutIndex = 22;
1171  cutVariableValue = dd0pr2;
1172  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1173  if(!bPassedCut && !fGetCutInfo) return 0;
1174  //---------------------------------------------------------------------
1175 
1176  // "topomatic max" ----------------------------------------------------
1177  nCutIndex = 23;
1178  cutVariableValue = dd0max;
1179  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1180  if(!bPassedCut && !fGetCutInfo) return 0;
1181  //---------------------------------------------------------------------
1182 
1183  // "topomatic min" ----------------------------------------------------
1184  nCutIndex = 24;
1185  cutVariableValue = dd0min;
1186  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1187  if(!bPassedCut && !fGetCutInfo) return 0;
1188  //---------------------------------------------------------------------
1189 
1190  // "pointing angle XY" ----------------------------------------------------
1191  nCutIndex = 25;
1192  cutVariableValue = cosPointingAngleXY;
1193  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1194  if(!bPassedCut && !fGetCutInfo) return 0;
1195  //---------------------------------------------------------------------
1196 
1197  // "vertex distance XY" ----------------------------------------------------
1198  nCutIndex = 26;
1199  cutVariableValue = distanceXYToVertex;
1200  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1201  if(!bPassedCut && !fGetCutInfo) return 0;
1202  //---------------------------------------------------------------------
1203 
1204  // "normalized decay length XY" ----------------------------------------------------
1205  nCutIndex = 27;
1206  cutVariableValue = normalizedDecayLengthXY;
1207  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1208  if(!bPassedCut && !fGetCutInfo) return 0;
1209  //---------------------------------------------------------------------
1210 
1211  // "chi squared per NDF" ----------------------------------------------------
1212  nCutIndex = 28;
1213  cutVariableValue = chi2Vertex;
1214  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1215  if(!bPassedCut && !fGetCutInfo) return 0;
1216  //---------------------------------------------------------------------
1217 
1218 
1219 
1220  AliAODRecoDecay* candidateD0toBPlus = (AliAODRecoDecay*)candidateD0;
1221  AliExternalTrackParam pionD0Track;
1222  AliExternalTrackParam kaonD0Track;
1223 
1224  Double_t d0z0DSVert[2],covd0z0DSVert[3],d0DSVert[2];
1225 
1226  pionD0Track.CopyFromVTrack(candidatePion);
1227  pionD0Track.PropagateToDCA(vertexBPlus,bz,100.,d0z0DSVert,covd0z0DSVert);
1228  d0DSVert[0] = d0z0DSVert[0];
1229 
1230  kaonD0Track.CopyFromVTrack(candidateKaon);
1231  kaonD0Track.PropagateToDCA(vertexBPlus,bz,100.,d0z0DSVert,covd0z0DSVert);
1232  d0DSVert[1] = d0z0DSVert[0];
1233 
1234  AliExternalTrackParam D0Track;
1235  D0Track.CopyFromVTrack(candidateD0);
1236  Double_t d0z0D0DSVert[2],covd0z0D0DSVert[3],d0D0DSVert;
1237  motherTrack.PropagateToDCA(vertexBPlus,bz,100.,d0z0D0DSVert,covd0z0D0DSVert);
1238  d0D0DSVert = TMath::Abs(d0z0D0DSVert[0]);
1239 
1240  Double_t impactProductToBPlus = d0DSVert[0]*d0DSVert[1];
1241  Double_t impactProductXYToBPlus = candidateD0toBPlus->ImpParXY(vertexBPlus);
1242 
1243  Double_t pointingAngleToBPlus = candidateD0toBPlus->CosPointingAngle(vertexBPlus);
1244  Double_t d0FirstDaughterToBPlus = TMath::Abs(d0DSVert[0]);
1245  Double_t d0SecondDaughterToBPlus = TMath::Abs(d0DSVert[1]);
1246  Double_t normDecayLengthToBPlus = candidateD0toBPlus->NormalizedDecayLength(vertexBPlus);
1247 
1248  Double_t pseudoProperDecayLengthDSVert = ((vertexD0->GetX() - vertexBPlus->GetX()) * candidateD0->Px() / TMath::Abs(candidateD0->Pt())) + ((vertexD0->GetY() - vertexBPlus->GetY()) * candidateD0->Py() / TMath::Abs(candidateD0->Pt()));
1249  Double_t pseudoProperDecayTimeToBPlus = pseudoProperDecayLengthDSVert * pdgMassMother/ptMother;
1250  Double_t DecayTimeToBPlus = vertexDistance / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother/(momentumMother*momentumMother)) + 1)));
1251 
1252  Double_t phiDSVert = candidateD0->Phi();
1253  Double_t thetaDSVert = candidateD0->Theta();
1254  Double_t covMatrixDSVert[21];
1255  candidateD0->GetCovarianceXYZPxPyPz(covMatrixDSVert);
1256 
1257  cp = TMath::Cos(phiDSVert);
1258  sp = TMath::Sin(phiDSVert);
1259  ct = TMath::Cos(thetaDSVert);
1260  st = TMath::Sin(thetaDSVert);
1261 
1262  errorMomentum = covMatrix[9]*cp*cp*ct*ct // GetCovPxPx
1263  +covMatrix[13]*2.*cp*sp*ct*ct // GetCovPxPy
1264  +covMatrix[18]*2.*cp*ct*st // GetCovPxPz
1265  +covMatrix[14]*sp*sp*ct*ct // GetCovPyPy
1266  +covMatrix[19]*2.*sp*ct*st // GetCovPyPz
1267  +covMatrix[20]*st*st; // GetCovPzPz
1268  Double_t normalizedDecayTimeToBPlus = candidateD0toBPlus->NormalizedDecayLength(vertexBPlus) / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother*errorMomentum*errorMomentum/(momentumMother*momentumMother)) + 1)));
1269 
1270  // "pointingAngleToBPlus" ---------------------------------------------
1271  nCutIndex = 29;
1272  cutVariableValue = pointingAngleToBPlus;
1273  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1274  if(!bPassedCut && !fGetCutInfo) return 0;
1275  //---------------------------------------------------------------------
1276 
1277  // "d0MotherToBPlus" --------------------------------------------------
1278  nCutIndex = 30;
1279  cutVariableValue = d0D0DSVert;
1280  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1281  if(!bPassedCut && !fGetCutInfo) return 0;
1282  //---------------------------------------------------------------------
1283 
1284  // "d0FirstDaughterToBPlus" -------------------------------------------
1285  nCutIndex = 31;
1286  cutVariableValue = d0FirstDaughterToBPlus;
1287  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1288  if(!bPassedCut && !fGetCutInfo) return 0;
1289  //---------------------------------------------------------------------
1290 
1291  // "d0SecondDaughterToBPlus" ------------------------------------------
1292  nCutIndex = 32;
1293  cutVariableValue = d0SecondDaughterToBPlus;
1294  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1295  if(!bPassedCut && !fGetCutInfo) return 0;
1296  //---------------------------------------------------------------------
1297 
1298  // "impactProductToBPlus" ---------------------------------------------
1299  nCutIndex = 33;
1300  cutVariableValue = impactProductToBPlus;
1301  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1302  if(!bPassedCut && !fGetCutInfo) return 0;
1303  //---------------------------------------------------------------------
1304 
1305  // "impactProductXYToBPlus" -------------------------------------------
1306  nCutIndex = 34;
1307  cutVariableValue = impactProductXYToBPlus;
1308  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1309  if(!bPassedCut && !fGetCutInfo) return 0;
1310  //---------------------------------------------------------------------
1311 
1312  // "normDecayLengthToBPlus" -------------------------------------------
1313  nCutIndex = 35;
1314  cutVariableValue = normDecayLengthToBPlus;
1315  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1316  if(!bPassedCut && !fGetCutInfo) return 0;
1317  //---------------------------------------------------------------------
1318 
1319  // "pseudoProperDecayTimeToBPlus" -------------------------------------
1320  nCutIndex = 36;
1321  cutVariableValue = pseudoProperDecayTimeToBPlus;
1322  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1323  if(!bPassedCut && !fGetCutInfo) return 0;
1324  //---------------------------------------------------------------------
1325 
1326  // "DecayTimeToBPlus" -------------------------------------------------
1327  nCutIndex = 37;
1328  cutVariableValue = DecayTimeToBPlus;
1329  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1330  if(!bPassedCut && !fGetCutInfo) return 0;
1331  //---------------------------------------------------------------------
1332 
1333  // "normalizedDecayTimeToBPlus" ---------------------------------------------
1334  nCutIndex = 38;
1335  cutVariableValue = normalizedDecayTimeToBPlus;
1336  bPassedCut = ApplyCutOnVariable(nCutIndex,ptbin,cutVariableValue,bCutArray);
1337  if(!bPassedCut && !fGetCutInfo) return 0;
1338  //---------------------------------------------------------------------
1339  }
1340 
1341  for (Int_t i = 0; i < 39; ++i)
1342  {
1343  if(bCutArray[i]==kTRUE){
1344  returnvalue = 0;
1345  break;
1346  }
1347  }
1348 
1349  return returnvalue;
1350 }
1351 //----------------------------------------------------------------------------------
1353  //
1354  // Apply selection on D0 candidate.
1355  //
1356 
1357  if(!fCutsRDD0forD0ptbin){
1358  cout<<"Cut matrice not inizialized. Exit..."<<endl;
1359  return 0;
1360  }
1361 
1362  AliAODRecoDecayHF2Prong* candidateD0 = (AliAODRecoDecayHF2Prong*)obj;
1363  if(!candidateD0){
1364  cout<<"candidateD0 null"<<endl;
1365  return 0;
1366  }
1367 
1368  AliAODTrack *candidatePion = (AliAODTrack*)candidateD0->GetDaughter(0);
1369  if(!candidatePion){
1370  cout<<"candidatePion null 3"<<endl;
1371  return 0;
1372  }
1373 
1374  AliAODTrack *candidateKaon = (AliAODTrack*)candidateD0->GetDaughter(1);
1375  if(!candidateKaon){
1376  cout<<"candidateKaon null"<<endl;
1377  return 0;
1378  }
1379 
1380  AliAODVertex * vertexD0 = candidateD0->GetSecondaryVtx();
1381  if(!vertexD0){
1382  cout<<"vertexD0 null"<<endl;
1383  return 0;
1384  }
1385 
1386  AliAODVertex * primaryVertex = aod->GetPrimaryVertex();
1387  if(!primaryVertex){
1388  cout<<"primaryVertex null"<<endl;
1389  return 0;
1390  }
1391 
1392  Int_t returnvalue=1;
1393  Bool_t bPassedCut = kFALSE;
1394 
1395  //get the magnetic field
1396  Double_t bz = (Double_t)aod->GetMagneticField();
1397 
1398 
1399  // selection on candidate
1400  if(selectionLevel==AliRDHFCuts::kAll ||
1401  selectionLevel==AliRDHFCuts::kCandidate) {
1402 
1403  Int_t ptbin=PtBinD0forD0ptbin(candidateD0->Pt());
1404  if(ptbin==-1) return -1;
1405 
1406  // D0mass
1407  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1408 
1409  // D0 window - invariant mass
1410  UInt_t prongs[2];
1411  prongs[0] = 211; prongs[1] = 321;
1412  Double_t invMassD0 = candidateD0->InvMass(2,prongs);
1413  Double_t invMassDifference = TMath::Abs(mD0PDG - invMassD0);
1414 
1415  UInt_t prongs2[2];
1416  prongs2[1] = 211; prongs2[0] = 321;
1417  Double_t invMassD02 = candidateD0->InvMass(2,prongs2);
1418  Double_t invMassDifference2 = TMath::Abs(mD0PDG - invMassD02);
1419 
1420  if(invMassDifference > invMassDifference2) invMassDifference = invMassDifference2;
1421 
1422  Double_t pointingAngle = candidateD0->CosPointingAngle();
1423  Double_t dcaMother = candidateD0->GetDCA();
1424  Double_t ptMother = candidateD0->Pt();
1425  Double_t momentumMother = candidateD0->P();
1426  Double_t ptPion = candidatePion->Pt();
1427  Double_t ptKaon = candidateKaon->Pt();
1428 
1429  AliExternalTrackParam motherTrack;
1430  motherTrack.CopyFromVTrack(candidateD0);
1431  Double_t d0z0[2],covd0z0[3],d0[2];
1432  motherTrack.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
1433  d0[0] = d0z0[0];
1434  Double_t d0Mother = TMath::Abs(d0[0]);
1435  Double_t d0firstTrack = TMath::Abs(candidateD0->Getd0Prong(0));
1436  Double_t d0secondTrack = TMath::Abs(candidateD0->Getd0Prong(1));
1437 
1438  Double_t impactProduct = candidateD0->Prodd0d0();
1439  Double_t impactProductXY = TMath::Abs(candidateD0->ImpParXY());
1440 
1441  Double_t angleBetweenBothDaughters = (candidateKaon->Px() * candidatePion->Px() + candidateKaon->Py() * candidatePion->Py() + candidateKaon->Pz() * candidatePion->Pz()) /(candidateKaon->P() * candidatePion->P());
1442  Double_t angleMotherFirstDaughter = (candidateD0->Px() * candidatePion->Px() + candidateD0->Py() * candidatePion->Py() + candidateD0->Pz() * candidatePion->Pz()) /(candidateD0->P() * candidatePion->P());
1443  Double_t angleMotherSecondDaughter = (candidateD0->Px() * candidateKaon->Px() + candidateD0->Py() * candidateKaon->Py() + candidateD0->Pz() * candidateKaon->Pz()) /(candidateD0->P() * candidateKaon->P());
1444 
1445  Double_t cosThetaStar = candidateD0->CosThetaStar(0,421,211,321);
1446  Double_t vertexDistance = vertexD0->DistanceToVertex(primaryVertex);
1447  Double_t normDecayLength = candidateD0->NormalizedDecayLength();
1448  Double_t pdgMassMother = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1449  Double_t pseudoProperDecayLength = ((vertexD0->GetX() - primaryVertex->GetX()) * candidateD0->Px() / TMath::Abs(candidateD0->Pt())) + ((vertexD0->GetY() - primaryVertex->GetY()) * candidateD0->Py() / TMath::Abs(candidateD0->Pt()));
1450  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother/ptMother;
1451  Double_t decayTime = vertexDistance / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother/(momentumMother*momentumMother)) + 1)));
1452 
1453  Double_t phi = candidateD0->Phi();
1454  Double_t theta = candidateD0->Theta();
1455  Double_t covMatrix[21];
1456  candidateD0->GetCovarianceXYZPxPyPz(covMatrix);
1457 
1458  Double_t cp = TMath::Cos(phi);
1459  Double_t sp = TMath::Sin(phi);
1460  Double_t ct = TMath::Cos(theta);
1461  Double_t st = TMath::Sin(theta);
1462 
1463  Double_t errorMomentum = covMatrix[9]*cp*cp*ct*ct // GetCovPxPx
1464  +covMatrix[13]*2.*cp*sp*ct*ct // GetCovPxPy
1465  +covMatrix[18]*2.*cp*ct*st // GetCovPxPz
1466  +covMatrix[14]*sp*sp*ct*ct // GetCovPyPy
1467  +covMatrix[19]*2.*sp*ct*st // GetCovPyPz
1468  +covMatrix[20]*st*st; // GetCovPzPz
1469  Double_t normalizedDecayTime = candidateD0->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother*errorMomentum*errorMomentum/(momentumMother*momentumMother)) + 1)));
1470 
1471  Double_t cosPointingAngleXY = candidateD0->CosPointingAngleXY();
1472  Double_t distanceXYToVertex = vertexD0->DistanceXYToVertex(primaryVertex);
1473  Double_t normalizedDecayLengthXY = candidateD0->NormalizedDecayLengthXY();
1474  Double_t chi2Vertex = vertexD0->GetChi2perNDF();
1475 
1476  //Topomatic
1477  Double_t dd0pr1=0.;
1478  Double_t dd0pr2=0.;
1479  Double_t dd0max=0.;
1480  Double_t dd0min=0.;
1481  for(Int_t ipr=0; ipr<2; ipr++)
1482  {
1483  Double_t diffIP, errdiffIP;
1484  candidateD0->Getd0MeasMinusExpProng(ipr,bz,diffIP,errdiffIP);
1485  Double_t normdd0=0.;
1486  if(errdiffIP>0.) normdd0=diffIP/errdiffIP;
1487  if(ipr==0) dd0pr1=normdd0;
1488  if(ipr==1) dd0pr2=normdd0;
1489  }
1490  if(TMath::Abs(dd0pr1)>TMath::Abs(dd0pr2)) {dd0max=dd0pr1; dd0min=dd0pr2;}
1491  else {dd0max=dd0pr2; dd0min=dd0pr1;}
1492 
1493 
1494  // We apply the cuts
1495  Int_t nCutIndex = 0;
1496  Double_t cutVariableValue = 0.0;
1497 
1498  // "inv. mass width [GeV]" --------------------------------------------
1499  nCutIndex = 0;
1500  cutVariableValue = invMassDifference;
1501  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1502  if(!bPassedCut && !fGetCutInfo) return 0;
1503  //---------------------------------------------------------------------
1504 
1505  // "delta mass width [GeV]" -------------------------------------------
1506  // nCutIndex = 1; // not used for D0
1507  // cutVariableValue = invMassDelta;
1508  // bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1509  // if(!bPassedCut && !fGetCutInfo) return 0;
1510  //---------------------------------------------------------------------
1511 
1512  // "pointing angle [Cos(theta)]" --------------------------------------
1513  nCutIndex = 2;
1514  cutVariableValue = pointingAngle;
1515  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1516  if(!bPassedCut && !fGetCutInfo) return 0;
1517  //---------------------------------------------------------------------
1518 
1519  // "dca [cm]" ---------------------------------------------------------
1520  nCutIndex = 3;
1521  cutVariableValue = dcaMother;
1522  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1523  if(!bPassedCut && !fGetCutInfo) return 0;
1524  //---------------------------------------------------------------------
1525 
1526  // "Pt D0 [GeV/c]" ----------------------------------------------------
1527  nCutIndex = 4;
1528  cutVariableValue = ptMother;
1529  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1530  if(!bPassedCut && !fGetCutInfo) return 0;
1531  //---------------------------------------------------------------------
1532 
1533  // "Pt Kaon [GeV/c]" -------------------------------------------------
1534  nCutIndex = 5;
1535  cutVariableValue = ptKaon;
1536  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1537  if(!bPassedCut && !fGetCutInfo) return 0;
1538  //---------------------------------------------------------------------
1539 
1540  // "Pt Pion [GeV/c]" --------------------------------------------------
1541  nCutIndex = 6;
1542  cutVariableValue = ptPion;
1543  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1544  if(!bPassedCut && !fGetCutInfo) return 0;
1545  //---------------------------------------------------------------------
1546 
1547  // "d0 D0 [cm]" -------------------------------------------------------
1548  nCutIndex = 7;
1549  cutVariableValue = d0Mother;
1550  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1551  if(!bPassedCut && !fGetCutInfo) return 0;
1552  //---------------------------------------------------------------------
1553 
1554  // "d0 Kaon [cm]"-----------------------------------------------------
1555  nCutIndex = 8;
1556  cutVariableValue = d0firstTrack;
1557  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1558  if(!bPassedCut && !fGetCutInfo) return 0;
1559  //---------------------------------------------------------------------
1560 
1561  // "d0 Pion [cm]" -----------------------------------------------------
1562  nCutIndex = 9;
1563  cutVariableValue = d0secondTrack;
1564  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1565  if(!bPassedCut && !fGetCutInfo) return 0;
1566  //---------------------------------------------------------------------
1567 
1568  // "d0d0 [cm^2]" ------------------------------------------------------
1569  nCutIndex = 10;
1570  cutVariableValue = impactProduct;
1571  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1572  if(!bPassedCut && !fGetCutInfo) return 0;
1573  //---------------------------------------------------------------------
1574 
1575  // "d0d0 XY [cm^2]" ---------------------------------------------------
1576  nCutIndex = 11;
1577  cutVariableValue = impactProductXY;
1578  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1579  if(!bPassedCut && !fGetCutInfo) return 0;
1580  //---------------------------------------------------------------------
1581 
1582  // "angle between both daughters" -------------------------------------
1583  nCutIndex = 12;
1584  cutVariableValue = angleBetweenBothDaughters;
1585  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1586  if(!bPassedCut && !fGetCutInfo) return 0;
1587  //---------------------------------------------------------------------
1588 
1589  // "angle mother with first daughter" ---------------------------------
1590  nCutIndex = 13;
1591  cutVariableValue = angleMotherFirstDaughter;
1592  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1593  if(!bPassedCut && !fGetCutInfo) return 0;
1594  //---------------------------------------------------------------------
1595 
1596  // "angle mother with second daughter" --------------------------------
1597  nCutIndex = 14;
1598  cutVariableValue = angleMotherSecondDaughter;
1599  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1600  if(!bPassedCut && !fGetCutInfo) return 0;
1601  //---------------------------------------------------------------------
1602 
1603  // "cosThetaStar" -----------------------------------------------------
1604  nCutIndex = 15;
1605  cutVariableValue = cosThetaStar;
1606  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1607  if(!bPassedCut && !fGetCutInfo) return 0;
1608  //---------------------------------------------------------------------
1609 
1610  // "vertexDistance" ---------------------------------------------------
1611  nCutIndex = 16;
1612  cutVariableValue = vertexDistance;
1613  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1614  if(!bPassedCut && !fGetCutInfo) return 0;
1615  //---------------------------------------------------------------------
1616 
1617  // "pseudoProperDecayTime" --------------------------------------------
1618  nCutIndex = 17;
1619  cutVariableValue = pseudoProperDecayTime;
1620  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1621  if(!bPassedCut && !fGetCutInfo) return 0;
1622  //---------------------------------------------------------------------
1623 
1624  // "DecayTime" --------------------------------------------------------
1625  nCutIndex = 18;
1626  cutVariableValue = decayTime;
1627  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1628  if(!bPassedCut && !fGetCutInfo) return 0;
1629  //---------------------------------------------------------------------
1630 
1631  // "normalizedDecayTime" ----------------------------------------------------
1632  nCutIndex = 19;
1633  cutVariableValue = normalizedDecayTime;
1634  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1635  if(!bPassedCut && !fGetCutInfo) return 0;
1636  //---------------------------------------------------------------------
1637 
1638  // "normDecayLength" --------------------------------------------------
1639  nCutIndex = 20;
1640  cutVariableValue = normDecayLength;
1641  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1642  if(!bPassedCut && !fGetCutInfo) return 0;
1643  //---------------------------------------------------------------------
1644 
1645  // "topomatic first daughter" -----------------------------------------
1646  nCutIndex = 21;
1647  cutVariableValue = dd0pr1;
1648  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1649  if(!bPassedCut && !fGetCutInfo) return 0;
1650  //---------------------------------------------------------------------
1651 
1652  // "topomatic second daughter" ----------------------------------------
1653  nCutIndex = 22;
1654  cutVariableValue = dd0pr2;
1655  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1656  if(!bPassedCut && !fGetCutInfo) return 0;
1657  //---------------------------------------------------------------------
1658 
1659  // "topomatic max" ----------------------------------------------------
1660  nCutIndex = 23;
1661  cutVariableValue = dd0max;
1662  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1663  if(!bPassedCut && !fGetCutInfo) return 0;
1664  //---------------------------------------------------------------------
1665 
1666  // "topomatic min" ----------------------------------------------------
1667  nCutIndex = 24;
1668  cutVariableValue = dd0min;
1669  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1670  if(!bPassedCut && !fGetCutInfo) return 0;
1671  //---------------------------------------------------------------------
1672 
1673  // "pointing angle XY" ----------------------------------------------------
1674  nCutIndex = 25;
1675  cutVariableValue = cosPointingAngleXY;
1676  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1677  if(!bPassedCut && !fGetCutInfo) return 0;
1678  //---------------------------------------------------------------------
1679 
1680  // "vertex distance XY" ----------------------------------------------------
1681  nCutIndex = 26;
1682  cutVariableValue = distanceXYToVertex;
1683  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1684  if(!bPassedCut && !fGetCutInfo) return 0;
1685  //---------------------------------------------------------------------
1686 
1687  // "normalized decay length XY" ----------------------------------------------------
1688  nCutIndex = 27;
1689  cutVariableValue = normalizedDecayLengthXY;
1690  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1691  if(!bPassedCut && !fGetCutInfo) return 0;
1692  //---------------------------------------------------------------------
1693 
1694  // "chi squared per NDF" ----------------------------------------------------
1695  nCutIndex = 28;
1696  cutVariableValue = chi2Vertex;
1697  bPassedCut = ApplyCutOnVariableD0forD0ptbin(nCutIndex,ptbin,cutVariableValue,bCutArray);
1698  if(!bPassedCut && !fGetCutInfo) return 0;
1699  //---------------------------------------------------------------------
1700 
1701  }
1702 
1703  for (Int_t i = 0; i < 29; ++i)
1704  {
1705  if(bCutArray[i]==kTRUE){
1706  returnvalue = 0;
1707  break;
1708  }
1709  }
1710 
1711  return returnvalue;
1712 }
1713 //----------------------------------------------------------------------------------
1715 {
1716  //
1717  // D* fiducial acceptance region // not (yet) used for the BPlus
1718  //
1719 
1720  // if(fMaxRapidityCand>-998.){
1721  // if(TMath::Abs(y) > fMaxRapidityCand) return kFALSE;
1722  // else return kTRUE;
1723  // }
1724 
1725  // if(pt > 5.) {
1726  // // applying cut for pt > 5 GeV
1727  // AliDebug(4,Form("pt of D* = %f (> 5), cutting at |y| < 0.8\n",pt));
1728  // if (TMath::Abs(y) > 0.8){
1729  // return kFALSE;
1730  // }
1731  // } else {
1732  // // appliying smooth cut for pt < 5 GeV
1733  // Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
1734  // Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
1735  // AliDebug(2,Form("pt of D* = %f (< 5), cutting according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
1736  // if (y < minFiducialY || y > maxFiducialY){
1737  // return kFALSE;
1738  // }
1739  // }
1740 
1741  return kTRUE;
1742 }
1743 //_______________________________________________________________________________-
1745 {
1746  //
1747  // PID method, n sigma approach default // not used for BPlus, done seperately for each daughter
1748  //
1749 
1750  // AliAODRecoDecayHF2Prong* dstar = (AliAODRecoDecayHF2Prong*)obj;
1751  // if(!dstar){
1752  // cout<<"AliAODRecoDecayHF2Prong null"<<endl;
1753  // return 0;
1754  // }
1755 
1756  // if(!fUsePID || dstar->Pt() > fMaxPtPid) return 3;
1757 
1758  // AliAODRecoDecayHF2Prong* d0 = (AliAODRecoDecayHF2Prong*)dstar->Get2Prong();
1759  // if(!d0){
1760  // cout<<"AliAODRecoDecayHF2Prong null"<<endl;
1761  // return 0;
1762  // }
1763 
1764  // // here the PID
1765  // AliAODTrack *pos = (AliAODTrack*)dstar->Get2Prong()->GetDaughter(0);
1766  // AliAODTrack *neg = (AliAODTrack*)dstar->Get2Prong()->GetDaughter(1);
1767 
1768  // if (dstar->Charge()>0){
1769  // if(!SelectPID(pos,2)) return 0;//pion+
1770  // if(!SelectPID(neg,3)) return 0;//kaon-
1771  // }else{
1772  // if(!SelectPID(pos,3)) return 0;//kaon+
1773  // if(!SelectPID(neg,2)) return 0;//pion-
1774  // }
1775 
1776  // if ((fPidHF->GetMatch() == 10 || fPidHF->GetMatch() == 11) && fPidHF->GetITS()) { //ITS n sigma band
1777  // AliAODTrack *softPion = (AliAODTrack*) dstar->GetBachelor();
1778 
1779  // if (fPidHF->CheckBands(AliPID::kPion, AliPIDResponse::kITS, softPion) == -1) {
1780  // return 0;
1781  // }
1782  // }
1783 
1784  return 3;
1785 }
1786 //_______________________________________________________________________________-
1788 {
1789  //
1790  // here the PID
1791 
1792  Bool_t isParticle=kTRUE;
1793  if(!fUsePID) return isParticle;
1794  Int_t match = fPidHF->GetMatch();
1795 
1796  if(match == 1){//n-sigma
1797  Bool_t TPCon=TMath::Abs(2)>1e-4?kTRUE:kFALSE;
1798  Bool_t TOFon=TMath::Abs(3)>1e-4?kTRUE:kFALSE;
1799 
1800  Bool_t isTPC=kTRUE;
1801  Bool_t isTOF=kTRUE;
1802 
1803  if (TPCon){//TPC
1804  if(fPidHF->CheckStatus(track,"TPC")){
1805  if(type==2) isTPC=fPidHF->IsPionRaw(track,"TPC");
1806  if(type==3) isTPC=fPidHF->IsKaonRaw(track,"TPC");
1807  }
1808  }
1809  if (TOFon){//TOF
1810  if(fPidHF->CheckStatus(track,"TOF")){
1811  if(type==2) isTOF=fPidHF->IsPionRaw(track,"TOF");
1812  if(type==3) isTOF=fPidHF->IsKaonRaw(track,"TOF");
1813  }
1814  }
1815 
1816  //--------------------------------
1817  // cut on high momentum in the TPC
1818  //--------------------------------
1819  Double_t pPIDcut = track->P();
1820  if(pPIDcut>fTPCflag) isTPC=1;
1821 
1822  isParticle = isTPC&&isTOF;
1823  }
1824 
1825  if(match == 2){//bayesian
1826  //Double_t priors[5]={0.01,0.001,0.3,0.3,0.3};
1827  Double_t prob[5]={1.,1.,1.,1.,1.};
1828 
1829  //fPidHF->SetPriors(priors,5);
1830  // fPidHF->BayesianProbability(track,prob);
1831 
1832  Double_t max=0.;
1833  Int_t k=-1;
1834  for (Int_t i=0; i<5; i++) {
1835  if (prob[i]>max) {k=i; max=prob[i];}
1836  }
1837  isParticle = Bool_t(k==type);
1838  }
1839 
1840  if (match == 10 || match == 11) { //Assymetric PID using histograms
1841  Int_t checkTPC = fPidHF->CheckBands((AliPID::EParticleType) type, AliPIDResponse::kTPC, track);
1842  Int_t checkTOF = fPidHF->CheckBands((AliPID::EParticleType) type, AliPIDResponse::kTOF, track);
1843 
1844  isParticle = checkTPC >= 0 && checkTOF >= 0 ? kTRUE : kFALSE; //Standard: both at least compatible
1845  if (match == 11) { //Extra requirement: at least one identified
1846  isParticle = isParticle && checkTPC+checkTOF >= 1 ? kTRUE : kFALSE;
1847  }
1848  }
1849 
1850  if (match == 12) { //Circular cut
1851  Double_t nSigmaTOF = 0;
1852  Double_t nSigmaTPC = 0;
1853 
1854  Double_t radius = fCircRadius;
1855 
1856  isParticle = kTRUE;
1857  if (radius > 0) {
1858  Int_t TPCok = fPidHF->GetnSigmaTPC(track, type, nSigmaTPC);
1859  Int_t TOFok = fPidHF->GetnSigmaTOF(track, type, nSigmaTOF);
1860  if (TPCok != -1 && TOFok != -1) {
1861  //Both detectors gave PID information
1862  isParticle = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF) <= radius ? kTRUE : kFALSE;
1863  }
1864  else {
1865  //Maximum one detector gave PID information
1866  if (TPCok != -1) {
1867  isParticle = nSigmaTPC <= radius ? kTRUE : kFALSE;
1868  }
1869  if (TOFok != -1) {
1870  isParticle = nSigmaTOF <= radius ? kTRUE : kFALSE;
1871  }
1872  }
1873  }
1874  }
1875 
1876  return isParticle;
1877 
1878 }
1879 //-------------------------------------------------------------------------------------
1881 {
1885 
1886  AliAODRecoDecayHF2Prong * D0 = (AliAODRecoDecayHF2Prong*)BPlus->GetDaughter(1);
1887 
1888  Double_t e[3] = {0};
1889  if (BPlus->Charge() == -1)
1890  {
1891  e[0] = D0->EProng(0, 211);
1892  e[1] = D0->EProng(1, 321);
1893  } else if (BPlus->Charge() == 1) {
1894  e[0] = D0->EProng(0, 321);
1895  e[1] = D0->EProng(1, 211);
1896  }
1897  e[2] = BPlus->EProng(0, 211);
1898 
1899  Double_t esum = e[0] + e[1] + e[2];
1900  Double_t invMassBPlus = TMath::Sqrt(esum * esum - BPlus->P2());
1901 
1902  Double_t invMassD0 = -1;
1903 
1904  if (BPlus->Charge() == -1) {invMassD0 = D0->InvMassD0();}
1905  else {invMassD0 = D0->InvMassD0bar();}
1906  if (invMassD0 == -1) {cout << "wrong invmass delta D0 BPlus" << endl;}
1907 
1908  return invMassBPlus - invMassD0;
1909 }
1910 
1911 
1912 //---------------------------------------------------------------------------
1913 //
1914 // DO for D0 pt bin functions
1915 //
1916 //---------------------------------------------------------------------------
1917 
1919  //
1920  // store the cuts
1921  //
1922 
1923  if(nVars!=fnVarsD0forD0ptbin) {
1924  printf("Wrong number of variables: it has to be %d\n",fnVarsD0forD0ptbin);
1925  AliFatal("exiting");
1926  }
1927  if(nPtBins!=fnPtBinsD0forD0ptbin) {
1928  printf("Wrong number of pt bins: it has to be %d\n",fnPtBinsD0forD0ptbin);
1929  AliFatal("exiting");
1930  }
1931 
1933 
1934 
1935  for(Int_t iv=0; iv<fnVarsD0forD0ptbin; iv++)
1936  {
1937  for(Int_t ib=0; ib<fnPtBinsD0forD0ptbin; ib++)
1938  {
1939  //check
1940 
1942  {
1943  cout<<"Overflow, exit..."<<endl;
1944  return;
1945  }
1946 
1947  fCutsRDD0forD0ptbin[GetGlobalIndexD0forD0ptbin(iv,ib)] = cutsRDD0forD0ptbin[iv][ib];
1948 
1949  }
1950  }
1951 
1952  return;
1953 }
1954 //---------------------------------------------------------------------------
1956  //
1957  // store the cuts
1958  //
1959  if(glIndex != fGlobalIndexD0forD0ptbin){
1960  cout<<"Wrong array size: it has to be "<<fGlobalIndexD0forD0ptbin<<endl;
1961  AliFatal("exiting");
1962  }
1964 
1965  for(Int_t iGl=0;iGl<fGlobalIndexD0forD0ptbin;iGl++){
1966  fCutsRDD0forD0ptbin[iGl] = cutsRDD0forD0ptbin[iGl];
1967  }
1968  return;
1969 }
1970 //---------------------------------------------------------------------------
1972  //
1973  //give the pt bin where the pt lies.
1974  //
1975  Int_t ptbin=-1;
1976  if(pt<fPtBinLimitsD0forD0ptbin[0])return ptbin;
1977  for (Int_t i=0;i<fnPtBinsD0forD0ptbin;i++){
1978  if(pt<fPtBinLimitsD0forD0ptbin[i+1]) {
1979  ptbin=i;
1980  break;
1981  }
1982  }
1983  return ptbin;
1984 }
1985 //---------------------------------------------------------------------------
1987  // Set the pt bins
1988 
1990  delete [] fPtBinLimitsD0forD0ptbin;
1991  fPtBinLimitsD0forD0ptbin = NULL;
1992  printf("Changing the pt bins\n");
1993  }
1994 
1995  if(nPtBinLimits != fnPtBinsD0forD0ptbin+1){
1996  cout<<"Warning: ptBinLimits dimension "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBinsD0forD0ptbin+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
1997  SetNPtBinsD0forD0ptbin(nPtBinLimits-1);
1998  }
1999 
2000  fnPtBinLimitsD0forD0ptbin = nPtBinLimits;
2002  //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
2004  for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimitsD0forD0ptbin[ib]=ptBinLimits[ib];
2005  for(Int_t ib=0; ib<nPtBinLimits; ib++) std::cout << "limit " << ib << " = " << fPtBinLimitsD0forD0ptbin[ib] << std::endl;
2006 
2007  return;
2008 }
2009 //---------------------------------------------------------------------------
2010 Int_t AliRDHFCutsBPlustoD0Pi::ApplyCutOnVariable(Int_t nCutIndex, Int_t ptbin, Float_t cutVariableValue, Bool_t bCutArray[68]){
2011 
2012  if(GetIsCutUsed(nCutIndex,ptbin)==kTRUE)
2013  {
2014  Bool_t bCut = kFALSE;
2015  if(GetIsUpperCut(nCutIndex)==kTRUE)
2016  {
2017  if(cutVariableValue > fCutsRD[GetGlobalIndex(nCutIndex,ptbin)]) bCut = kTRUE;
2018  } else
2019  {
2020  if(cutVariableValue < fCutsRD[GetGlobalIndex(nCutIndex,ptbin)]) bCut = kTRUE;
2021  }
2022  if(bCut == kTRUE) {bCutArray[nCutIndex] = 1; return 0;}
2023  }
2024  return 1;
2025 }
2026 //---------------------------------------------------------------------------
2027 Int_t AliRDHFCutsBPlustoD0Pi::ApplyCutOnVariableD0forD0ptbin(Int_t nCutIndex, Int_t ptbin, Float_t cutVariableValue, Bool_t bCutArray[29]){
2028 
2029  // std::cout << "index: " << nCutIndex << ", ptbin: " << ptbin << ", used: " << GetIsCutUsedD0forD0ptbin(nCutIndex,ptbin) << ", upper: " << GetIsUpperCutD0forD0ptbin(nCutIndex) << ", value: " << cutVariableValue << ", cutvalue: " << fCutsRDD0forD0ptbin[GetGlobalIndexD0forD0ptbin(nCutIndex,ptbin)] << std::endl;
2030  if(GetIsCutUsedD0forD0ptbin(nCutIndex,ptbin)==kTRUE)
2031  {
2032  Bool_t bCut = kFALSE;
2033  if(GetIsUpperCutD0forD0ptbin(nCutIndex)==kTRUE)
2034  {
2035  if(cutVariableValue > fCutsRDD0forD0ptbin[GetGlobalIndexD0forD0ptbin(nCutIndex,ptbin)]) bCut = kTRUE;
2036  } else
2037  {
2038  if(cutVariableValue < fCutsRDD0forD0ptbin[GetGlobalIndexD0forD0ptbin(nCutIndex,ptbin)]) bCut = kTRUE;
2039  }
2040  if(bCut == kTRUE) {bCutArray[nCutIndex] = 1; return 0;}
2041  }
2042  return 1;
2043 }
2044 //---------------------------------------------------------------------------
2046  // Set the variable names
2047 
2048  if(fVarNamesD0forD0ptbin) {
2049  delete [] fVarNamesD0forD0ptbin;
2050  fVarNamesD0forD0ptbin = NULL;
2051  //printf("Changing the variable names\n");
2052  }
2053  if(nVars!=fnVarsD0forD0ptbin){
2054  printf("Wrong number of variables: it has to be %d\n",fnVarsD0forD0ptbin);
2055  return;
2056  }
2057  //fnVars=nVars;
2058  fVarNamesD0forD0ptbin = new TString[nVars];
2059  fIsUpperCutD0forD0ptbin = new Bool_t[nVars];
2060  for(Int_t iv=0; iv<nVars; iv++) {
2061  fVarNamesD0forD0ptbin[iv] = varNames[iv];
2062  fIsUpperCutD0forD0ptbin[iv] = isUpperCut[iv];
2063  }
2064 
2065  return;
2066 }
2067 //---------------------------------------------------------------------------
2069  // 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.
2070 
2071  if(fIsCutUsed) {
2072  delete [] fIsCutUsed;
2073  fIsCutUsed = NULL;
2074  }
2075 
2076  fIsCutUsed = new Bool_t[(GetNPtBins())*(GetNVars())];
2077  for(Int_t iv=0; iv<(GetNPtBins())*(GetNVars()); iv++) {
2078  fIsCutUsed[iv] = kFALSE;
2079  }
2080 
2081  if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
2082 
2083 
2084  for(Int_t iv=0; iv<fnVars; iv++) {
2085 
2086  for(Int_t ib=0; ib<fnPtBins; ib++) {
2087 
2088  //check
2089  if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
2090  cout<<"Overflow, exit..."<<endl;
2091  return;
2092  }
2093 
2094  fCutsRD[GetGlobalIndex(iv,ib)] = 0;
2095 
2096  }
2097  }
2098 
2099  // D0 for D0 pt bin
2100 
2102  delete [] fIsCutUsedD0forD0ptbin;
2103  fIsCutUsedD0forD0ptbin = NULL;
2104  }
2105 
2107  for(Int_t iv=0; iv<(GetNPtBinsD0forD0ptbin())*(GetNVarsD0forD0ptbin()); iv++) {
2108  fIsCutUsedD0forD0ptbin[iv] = kFALSE;
2109 
2110  }
2111 
2113 
2114  for(Int_t iv=0; iv<fnVarsD0forD0ptbin; iv++)
2115  {
2116  for(Int_t ib=0; ib<fnPtBinsD0forD0ptbin; ib++)
2117  {
2118  //check
2119 
2121  {
2122  cout<<"Overflow, exit..."<<endl;
2123  return;
2124  }
2125 
2127  }
2128  }
2129 
2130  return;
2131 }
2132 //---------------------------------------------------------------------------
2134  // Set the cut value and direction
2135 
2136  fIsCutUsed[GetGlobalIndex(nCutIndex,ptBin)] = kTRUE;
2137  fIsUpperCut[nCutIndex] = cutDirection;
2138  fCutsRD[GetGlobalIndex(nCutIndex,ptBin)] = cutValue;
2139 
2140  return;
2141 }
2142 //---------------------------------------------------------------------------
2144  // Set the cut value and direction
2145 
2146  fIsCutUsedD0forD0ptbin[GetGlobalIndexD0forD0ptbin(nCutIndex,ptBin)] = kTRUE;
2147  fIsUpperCutD0forD0ptbin[nCutIndex] = cutDirection;
2148  fCutsRDD0forD0ptbin[GetGlobalIndexD0forD0ptbin(nCutIndex,ptBin)] = cutValue;
2149 
2150  return;
2151 }
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:455
Int_t IsD0FromBPlusSelected(Double_t ptBPlus, TObject *obj, Int_t selectionLevel, AliAODEvent *aod, Bool_t bCutArray[68])
Double_t NormalizedDecayLength() const
Int_t GetNVars() const
Definition: AliRDHFCuts.h:250
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod, Bool_t *bCutArray)
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:433
Bool_t IsKaonRaw(AliAODTrack *track, TString detector) const
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
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)
Bool_t fHardSelectionArrayITSD0SecondDaughter[7]
void SetNVars(Int_t nVars)
Definition: AliRDHFCuts.h:404
Double_t CosPointingAngleXY() const
void SetIsCutUsedD0forD0ptbin(Int_t nCutIndex, Int_t ptbin, Bool_t isCutUsed)
Int_t IsD0forD0ptbinSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod, Bool_t *bCutArray)
AliRDHFCuts & operator=(const AliRDHFCuts &source)
Bool_t fUsePID
Definition: AliRDHFCuts.h:434
virtual void GetCutVarsForOpt(AliAODRecoDecayHF *d, Float_t *vars, Int_t nvars, Int_t *pdgdaughters)
const Int_t nPtBins
Int_t ApplyCutOnVariableD0forD0ptbin(Int_t nCutIndex, Int_t ptbin, Float_t cutVariableValue, Bool_t bCutArray[29])
Int_t CheckBands(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, AliAODTrack *track)
void SetNPtBinsD0forD0ptbin(Int_t nptBins)
Int_t PtBinD0forD0ptbin(Double_t pt) const
int Int_t
Definition: External.C:63
Bool_t * GetIsUpperCut() const
Definition: AliRDHFCuts.h:261
unsigned int UInt_t
Definition: External.C:33
virtual Int_t IsSelectedPID(AliAODRecoDecayHF *rd)
void SetSoftSelectionArrayITSD0SecondDaughter(const Bool_t array[7]=0)
float Float_t
Definition: External.C:68
void SetCut(Int_t nCutIndex, Int_t ptBin, AliRDHFCutsBPlustoD0Pi::EUpperCut cutDirection, Float_t cutValue)
Bool_t CheckStatus(AliAODTrack *track, TString detectors) const
void SetIsCutUsed(Int_t nCutIndex, Int_t ptbin, Bool_t isCutUsed)
AliRDHFCutsBPlustoD0Pi(const char *name="BPlustoD0PiCuts")
virtual Int_t SelectPID(AliAODTrack *track, Int_t type)
void SetSoftSelectionArrayITSBPlusPion(const Bool_t array[7]=0)
Int_t fIsSelectedPID
outcome of cuts selection
Definition: AliRDHFCuts.h:456
Float_t * fCutsRD
fnVars*fnPtBins
Definition: AliRDHFCuts.h:432
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)
Bool_t IsPionRaw(AliAODTrack *track, TString detector) const
void SetPtBins(Int_t nPtBinLimits, Float_t *ptBinLimits)
void SetHardSelectionArrayITSBPlusPion(const Bool_t array[7]=0)
void SetSoftSelectionArrayITSD0FirstDaughter(const Bool_t array[7]=0)
Int_t GetMatch() const
Definition: AliAODPidHF.h:145
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:249
void SetNVarsD0forD0ptbin(Int_t nVars)
void SetPtBinsD0forD0ptbin(Int_t nPtBinLimits, Float_t *ptBinLimits)
AliRDHFCutsBPlustoD0Pi & operator=(const AliRDHFCutsBPlustoD0Pi &source)
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
Int_t fnPtBins
cuts on the candidate
Definition: AliRDHFCuts.h:424
AliAODPidHF * fPidHF
enable AOD049 centrality cleanup
Definition: AliRDHFCuts.h:436
Int_t GetNPtBinsD0forD0ptbin() const
Int_t fGlobalIndex
Definition: AliRDHFCuts.h:431
Int_t PtBin(Double_t pt) const
Bool_t fSoftSelectionArrayITSD0SecondDaughter[7]
void SetCutD0forD0ptbin(Int_t nCutIndex, Int_t ptBin, AliRDHFCutsBPlustoD0Pi::EUpperCut cutDirection, Float_t cutValue)
Int_t GetGlobalIndex(Int_t iVar, Int_t iPtBin) const
Int_t ApplyCutOnVariable(Int_t nCutIndex, Int_t ptbin, Float_t cutVariableValue, Bool_t bCutArray[68])
void SetCutsD0forD0ptbin(Int_t nVars, Int_t nPtBins, Float_t **cutsRDD0forD0ptbin)
Bool_t fGetCutInfo
Radius for circular PID nsigma cut.
Int_t fnVars
Definition: AliRDHFCuts.h:427