AliPhysics  cdeda5a (cdeda5a)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliRDHFCutsD0toKpi.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$ */
17 
19 //
20 // Class for cuts on AOD reconstructed D0->Kpi
21 //
22 // Author: A.Dainese, andrea.dainese@pd.infn.it
24 
25 #include <TDatabasePDG.h>
26 #include <Riostream.h>
27 
28 #include "AliRDHFCutsD0toKpi.h"
30 #include "AliAODTrack.h"
31 #include "AliESDtrack.h"
32 #include "AliAODPid.h"
33 #include "AliTPCPIDResponse.h"
34 #include "AliAODVertex.h"
35 #include "AliKFVertex.h"
36 #include "AliKFParticle.h"
37 
38 using std::cout;
39 using std::endl;
40 
41 
45 
46 
47 //--------------------------------------------------------------------------
49  AliRDHFCuts(name),
50  fUseImpParDCut(0),
51  fMaxImpParD(0x0),
52  fUsed0MeasMinusExpCut(0),
53  fMaxd0MeasMinusExp(0x0),
54  fUseSpecialCuts(kFALSE),
55  fLowPt(kTRUE),
56  fDefaultPID(kFALSE),
57  fUseKF(kFALSE),
58  fPtLowPID(2.),
59  fPtMaxSpecialCuts(9999.),
60  fmaxPtrackForPID(9999),
61  fCombPID(kFALSE),
62  fnSpecies(AliPID::kSPECIES),
63  fWeightsPositive(0),
64  fWeightsNegative(0),
65  fProbThreshold(0.5),
66  fBayesianStrategy(0),
67  fBayesianCondition(0)
68 {
69  //
70  // Default Constructor
71  //
72  Int_t nvars=11;
73  SetNVars(nvars);
74  TString varNames[11]={"inv. mass [GeV]",
75  "dca [cm]",
76  "cosThetaStar",
77  "pTK [GeV/c]",
78  "pTPi [GeV/c]",
79  "d0K [cm]",
80  "d0Pi [cm]",
81  "d0d0 [cm^2]",
82  "cosThetaPoint",
83  "|cosThetaPointXY|",
84  "NormDecayLenghtXY"};
85  Bool_t isUpperCut[11]={kTRUE,
86  kTRUE,
87  kTRUE,
88  kFALSE,
89  kFALSE,
90  kTRUE,
91  kTRUE,
92  kTRUE,
93  kFALSE,
94  kFALSE,
95  kFALSE};
96  SetVarNames(nvars,varNames,isUpperCut);
97  Bool_t forOpt[11]={kFALSE,
98  kTRUE,
99  kTRUE,
100  kFALSE,
101  kFALSE,
102  kFALSE,
103  kFALSE,
104  kTRUE,
105  kTRUE,
106  kFALSE,
107  kFALSE};
108  SetVarsForOpt(4,forOpt);
109  Float_t limits[2]={0,999999999.};
110  SetPtBins(2,limits);
111 
112  fWeightsNegative = new Double_t[AliPID::kSPECIES];
113  fWeightsPositive = new Double_t[AliPID::kSPECIES];
114 
115  for (Int_t i = 0; i<AliPID::kSPECIES; i++) {
116  fWeightsPositive[i] = 0.;
117  fWeightsNegative[i] = 0.;
118  }
119 
120 }
121 //--------------------------------------------------------------------------
123  AliRDHFCuts(source),
124  fUseImpParDCut(source.fUseImpParDCut),
125  fMaxImpParD(0x0),
126  fUsed0MeasMinusExpCut(source.fUsed0MeasMinusExpCut),
127  fMaxd0MeasMinusExp(0x0),
128  fUseSpecialCuts(source.fUseSpecialCuts),
129  fLowPt(source.fLowPt),
130  fDefaultPID(source.fDefaultPID),
131  fUseKF(source.fUseKF),
132  fPtLowPID(source.fPtLowPID),
133  fPtMaxSpecialCuts(source.fPtMaxSpecialCuts),
134  fmaxPtrackForPID(source.fmaxPtrackForPID),
135  fCombPID(source.fCombPID),
136  fnSpecies(source.fnSpecies),
137  fWeightsPositive(source.fWeightsPositive),
138  fWeightsNegative(source.fWeightsNegative),
139  fProbThreshold(source.fProbThreshold),
140  fBayesianStrategy(source.fBayesianStrategy),
141  fBayesianCondition(source.fBayesianCondition)
142 {
143  //
144  // Copy constructor
145  //
147  if(source.fUseImpParDCut&&source.fMaxImpParD) SetImpParDCut(source.fUseImpParDCut,source.fMaxImpParD);
148 }
149 //--------------------------------------------------------------------------
151 {
152  //
153  // assignment operator
154  //
155  if(&source == this) return *this;
156 
157  AliRDHFCuts::operator=(source);
159  if(source.fUseImpParDCut&&source.fMaxImpParD) SetImpParDCut(source.fUseImpParDCut,source.fMaxImpParD);
163  fLowPt=source.fLowPt;
164  fDefaultPID=source.fDefaultPID;
165  fUseKF=source.fUseKF;
166  fPtLowPID=source.fPtLowPID;
169  fCombPID = source.fCombPID;
170  fnSpecies = source.fnSpecies;
176  return *this;
177 }
178 
179 
180 
181 //---------------------------------------------------------------------------
183 {
184  //
185  // Destructor
186  //
187  if (fWeightsNegative) {
188  delete[] fWeightsNegative;
189  fWeightsNegative = 0;
190  }
191  if (fWeightsPositive) {
192  delete[] fWeightsNegative;
193  fWeightsNegative = 0;
194  }
195  if (fMaxd0MeasMinusExp){
196  delete [] fMaxd0MeasMinusExp;
198  }
199  if (fMaxImpParD){
200  delete [] fMaxImpParD;
201  fMaxImpParD =0;
202  }
203 
204 }
205 
206 //---------------------------------------------------------------------------
208  // set the same selection on normalized d0Meas-d0Exp for all pt bins
209  Float_t cuts[fnPtBins];
210  for(Int_t ip=0;ip<fnPtBins;ip++){
211  cuts[ip]=value;
212  }
213  Setd0MeasMinusExpCut(fnPtBins,cuts);
214 }
215 
216 //---------------------------------------------------------------------------
218  //
219  // store the cuts
220  //
221  if(nPtBins!=fnPtBins) {
222  AliFatal(Form("Wrong number of pt bins: it has to be %d, exiting",fnPtBins));
223  }
225  else {
226  delete [] fMaxd0MeasMinusExp;
228  }
229  for(Int_t ib=0; ib<fnPtBins; ib++) fMaxd0MeasMinusExp[ib] = cutval[ib];
231  return;
232 }
233 //---------------------------------------------------------------------------
235  // set the same selection on normalized d0Meas-d0Exp for all pt bins
236  Float_t cuts[fnPtBins];
237  for(Int_t ip=0;ip<fnPtBins;ip++){
238  cuts[ip]=value;
239  }
240  SetImpParDCut(fnPtBins,cuts);
241 }
242 
243 //---------------------------------------------------------------------------
245  //
246  // store the cuts
247  //
248  if(nPtBins!=fnPtBins) {
249  AliFatal(Form("Wrong number of pt bins: it has to be %d\n",fnPtBins));
250  }
252  else {
253  delete [] fMaxImpParD;
255  }
256  for(Int_t ib=0; ib<fnPtBins; ib++) fMaxImpParD[ib] = cutval[ib];
258  return;
259 }
260 //_____________________________________________________________________________
262  //
263  // Fills in vars the values of the variables
264  //
265 
266  if(nvars!=fnVarsForOpt) {
267  printf("AliRDHFCutsD0toKpi::GetCutsVarsForOpt: wrong number of variables\n");
268  return;
269  }
270 
272 
273  //recalculate vertex w/o daughters
274  Bool_t cleanvtx=kFALSE;
275  AliAODVertex *origownvtx=0x0;
277  if(dd->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*dd->GetOwnPrimaryVtx());
278  cleanvtx=kTRUE;
279  if(!RecalcOwnPrimaryVtx(dd,aod)) {
280  CleanOwnPrimaryVtx(dd,aod,origownvtx);
281  cleanvtx=kFALSE;
282  }
283  }
284 
285  Int_t iter=-1;
286  if(fVarsForOpt[0]){
287  iter++;
288  if(TMath::Abs(pdgdaughters[0])==211) {
289  vars[iter]=dd->InvMassD0();
290  } else {
291  vars[iter]=dd->InvMassD0bar();
292  }
293  }
294  if(fVarsForOpt[1]){
295  iter++;
296  vars[iter]=dd->GetDCA();
297  }
298  if(fVarsForOpt[2]){
299  iter++;
300  if(TMath::Abs(pdgdaughters[0])==211) {
301  vars[iter] = dd->CosThetaStarD0();
302  } else {
303  vars[iter] = dd->CosThetaStarD0bar();
304  }
305  }
306  if(fVarsForOpt[3]){
307  iter++;
308  if(TMath::Abs(pdgdaughters[0])==321) {
309  vars[iter]=dd->PtProng(0);
310  }
311  else{
312  vars[iter]=dd->PtProng(1);
313  }
314  }
315  if(fVarsForOpt[4]){
316  iter++;
317  if(TMath::Abs(pdgdaughters[0])==211) {
318  vars[iter]=dd->PtProng(0);
319  }
320  else{
321  vars[iter]=dd->PtProng(1);
322  }
323  }
324  if(fVarsForOpt[5]){
325  iter++;
326  if(TMath::Abs(pdgdaughters[0])==321) {
327  vars[iter]=dd->Getd0Prong(0);
328  }
329  else{
330  vars[iter]=dd->Getd0Prong(1);
331  }
332  }
333  if(fVarsForOpt[6]){
334  iter++;
335  if(TMath::Abs(pdgdaughters[0])==211) {
336  vars[iter]=dd->Getd0Prong(0);
337  }
338  else{
339  vars[iter]=dd->Getd0Prong(1);
340  }
341  }
342  if(fVarsForOpt[7]){
343  iter++;
344  vars[iter]= dd->Prodd0d0();
345  }
346  if(fVarsForOpt[8]){
347  iter++;
348  vars[iter]=dd->CosPointingAngle();
349  }
350 
351  if(fVarsForOpt[9]){
352  iter++;
353  vars[iter]=TMath::Abs(dd->CosPointingAngleXY());
354  }
355 
356  if(fVarsForOpt[10]){
357  iter++;
358  vars[iter]=dd->NormalizedDecayLengthXY();
359  }
360 
361  if(cleanvtx)CleanOwnPrimaryVtx(dd,aod,origownvtx);
362 
363  return;
364 }
365 //---------------------------------------------------------------------------
367  //
368  // Apply selection
369  //
370 
371  fIsSelectedCuts=0;
372  fIsSelectedPID=0;
373 
375  AliFatal(Form("Wrong number of pt bins for D imp par cut: it has to be %d\n",fnPtBins));
376  }
378  AliFatal(Form("Wrong number of pt bins for D imp par cut: it has to be %d\n",fnPtBins));
379  }
380 
381  if(!fCutsRD){
382  cout<<"Cut matrice not inizialized. Exit..."<<endl;
383  return 0;
384  }
385  //PrintAll();
387  if(!d){
388  cout<<"AliAODRecoDecayHF2Prong null"<<endl;
389  return 0;
390  }
391 
392  if(fKeepSignalMC) if(IsSignalMC(d,aod,421)) return 3;
393 
394  Double_t ptD=d->Pt();
395  if(ptD<fMinPtCand) return 0;
396  if(ptD>fMaxPtCand) return 0;
397 
399 
400  // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
401  Int_t returnvaluePID=3;
402  Int_t returnvalueCuts=3;
403 
404  // selection on candidate
405  if(selectionLevel==AliRDHFCuts::kAll ||
406  selectionLevel==AliRDHFCuts::kCandidate) {
407 
408  if(!fUseKF) {
409 
410  //recalculate vertex w/o daughters
411  AliAODVertex *origownvtx=0x0;
413  if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
414  if(!RecalcOwnPrimaryVtx(d,aod)) {
415  CleanOwnPrimaryVtx(d,aod,origownvtx);
416  return 0;
417  }
418  }
419 
420  if(fUseMCVertex) {
421  if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
422  if(!SetMCPrimaryVtx(d,aod)) {
423  CleanOwnPrimaryVtx(d,aod,origownvtx);
424  return 0;
425  }
426  }
427 
428  Double_t pt=d->Pt();
429 
430  Int_t okD0=0,okD0bar=0;
431 
432  Int_t ptbin=PtBin(pt);
433  if (ptbin==-1) {
434  CleanOwnPrimaryVtx(d,aod,origownvtx);
435  return 0;
436  }
437 
438  Double_t mD0,mD0bar,ctsD0,ctsD0bar;
439  okD0=1; okD0bar=1;
440 
441  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
442 
443  d->InvMassD0(mD0,mD0bar);
444  if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
445  if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
446  if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
447 
448  if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
449 
450  if(d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(3,ptbin)]*fCutsRD[GetGlobalIndex(3,ptbin)] || d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(4,ptbin)]*fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
451  if(d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(3,ptbin)]*fCutsRD[GetGlobalIndex(3,ptbin)] || d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(4,ptbin)]*fCutsRD[GetGlobalIndex(4,ptbin)]) okD0bar = 0;
452  if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
453 
454  if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
455  TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
456  if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
457  TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
458  if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
459 
460  if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
461 
462  d->CosThetaStarD0(ctsD0,ctsD0bar);
463  if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
464  if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
465  if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
466 
467  if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
468 
469  if(TMath::Abs(d->CosPointingAngleXY()) < fCutsRD[GetGlobalIndex(9,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
470 
471  Double_t normalDecayLengXY=d->NormalizedDecayLengthXY();
472  if (normalDecayLengXY < fCutsRD[GetGlobalIndex(10, ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
473 
474  if (returnvalueCuts!=0) {
475  if (okD0) returnvalueCuts=1; //cuts passed as D0
476  if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar
477  if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
478  }
479  //Topomatic
481  Double_t dd0max=0.;
482  for(Int_t ipr=0; ipr<2; ipr++) {
483  Double_t diffIP, errdiffIP;
484  d->Getd0MeasMinusExpProng(ipr,aod->GetMagneticField(),diffIP,errdiffIP);
485  Double_t normdd0=0.;
486  if(errdiffIP>0.) normdd0=diffIP/errdiffIP;
487  if(ipr==0) dd0max=normdd0;
488  else if(TMath::Abs(normdd0)>TMath::Abs(dd0max)) dd0max=normdd0;
489  }
490  if(fMaxd0MeasMinusExp[ptbin]>=0.){
491  if(TMath::Abs(dd0max)>fMaxd0MeasMinusExp[ptbin]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
492  }else{
493  if(TMath::Abs(dd0max)<TMath::Abs(fMaxd0MeasMinusExp[ptbin])) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
494  }
495  }
496  //ImpactParameterCand
497  if(fUseImpParDCut){
498  Double_t d0=d->ImpParXY()*10000.;//cut value in cm.
499  if(fMaxImpParD[ptbin]>=0.){//keep the region within the cut
500  if(TMath::Abs(d0)>fMaxImpParD[ptbin]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
501  }else{//keep the region outside the cut
502  if(TMath::Abs(d0)<TMath::Abs(fMaxImpParD[ptbin])) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
503  }
504  }//impact parameter cut
505 
506  // call special cuts
507  Int_t special=1;
509  if(!special) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
510 
511  // unset recalculated primary vertex when not needed any more
512  CleanOwnPrimaryVtx(d,aod,origownvtx);
513 
514  } else {
515  // go to selection with Kalman vertexing, if requested
516  returnvalueCuts = IsSelectedKF(d,aod);
517  }
518 
519  fIsSelectedCuts=returnvalueCuts;
520  if(!returnvalueCuts) return 0;
521  }
522 
523  // selection on PID
524  if(selectionLevel==AliRDHFCuts::kAll ||
525  selectionLevel==AliRDHFCuts::kCandidate ||
526  selectionLevel==AliRDHFCuts::kPID) {
527  if (!fCombPID) {
528  returnvaluePID = IsSelectedPID(d);
529  fIsSelectedPID=returnvaluePID;
530  if(!returnvaluePID) return 0;
531  } else {
532  returnvaluePID = IsSelectedCombPID(d);
533  if(!returnvaluePID) return 0;
534  }
535  }
536 
537  Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
538 
539  if(!returnvalueComb) return 0;
540 
541  // selection on daughter tracks
542  if(selectionLevel==AliRDHFCuts::kAll ||
543  selectionLevel==AliRDHFCuts::kTracks) {
544  if(!AreDaughtersSelected(d,aod)) return 0;
545  }
546 
547  // cout<<"Pid = "<<returnvaluePID<<endl;
548  return returnvalueComb;
549 }
550 
551 //------------------------------------------------------------------------------------------
553  //
554  // Apply selection using KF-vertexing
555  //
556 
557  AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
558  AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
559 
560  if(!track0 || !track1) {
561  cout<<"one or two D0 daughters missing!"<<endl;
562  return 0;
563  }
564 
565  // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
566  Int_t returnvalueCuts=3;
567 
568  // candidate selection with AliKF
569  AliKFParticle::SetField(aod->GetMagneticField()); // set the magnetic field
570 
571  Int_t okD0=0,okD0bar=0;
572  okD0=1; okD0bar=1;
573 
574  // convert tracks into AliKFParticles
575 
576  AliKFParticle negPiKF(*track1,-211); // neg pion kandidate
577  AliKFParticle negKKF(*track1,-321); // neg kaon kandidate
578  AliKFParticle posPiKF(*track0,211); // pos pion kandidate
579  AliKFParticle posKKF(*track0,321); // pos kaon kandidate
580 
581  // build D0 candidates
582 
583  AliKFParticle d0c(negKKF,posPiKF); // D0 candidate
584  AliKFParticle ad0c(posKKF,negPiKF); // D0bar candidate
585 
586  // create kf primary vertices
587 
588  AliAODVertex *vtx1 = aod->GetPrimaryVertex();
589  AliKFVertex primVtx(*vtx1);
590  AliKFVertex aprimVtx(*vtx1);
591 
592  if(primVtx.GetNContributors()<=0) okD0 = 0;
593  if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
594  if(!okD0 && !okD0bar) returnvalueCuts=0;
595 
596  // calculate mass
597 
598  Double_t d0mass = d0c.GetMass();
599  Double_t ad0mass = ad0c.GetMass();
600 
601  // calculate P of D0 and D0bar
602  Double_t d0P = d0c.GetP();
603  Double_t d0Px = d0c.GetPx();
604  Double_t d0Py = d0c.GetPy();
605  Double_t d0Pz = d0c.GetPz();
606  Double_t ad0P = ad0c.GetP();
607  Double_t ad0Px = ad0c.GetPx();
608  Double_t ad0Py = ad0c.GetPy();
609  Double_t ad0Pz = ad0c.GetPz();
610 
611  //calculate Pt of D0 and D0bar
612 
613  Double_t pt=d0c.GetPt();
614  Double_t apt=ad0c.GetPt();
615 
616  // remove D0 daughters from primary vertices (if used in vertex fit) and add D0-candidates
617 
618  if(track0->GetUsedForPrimVtxFit()) {
619  primVtx -= posPiKF;
620  aprimVtx -= posKKF;
621  }
622 
623  if(track1->GetUsedForPrimVtxFit()) {
624  primVtx -= negKKF;
625  aprimVtx -= negPiKF;
626  }
627 
628  primVtx += d0c;
629  aprimVtx += ad0c;
630 
631  if(primVtx.GetNContributors()<=0) okD0 = 0;
632  if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
633  if(!okD0 && !okD0bar) returnvalueCuts=0;
634 
635  //calculate cut variables
636 
637  // calculate impact params of daughters w.r.t recalculated vertices
638 
639  Double_t impactPi = posPiKF.GetDistanceFromVertexXY(primVtx);
640  Double_t aimpactPi = negPiKF.GetDistanceFromVertexXY(aprimVtx);
641  Double_t impactKa = negKKF.GetDistanceFromVertexXY(primVtx);
642  Double_t aimpactKa = posKKF.GetDistanceFromVertexXY(aprimVtx);
643 
644  // calculate Product of Impact Params
645 
646  Double_t prodParam = impactPi*impactKa;
647  Double_t aprodParam = aimpactPi*aimpactKa;
648 
649  // calculate cosine of pointing angles
650 
651  TVector3 mom(d0c.GetPx(),d0c.GetPy(),d0c.GetPz());
652  TVector3 fline(d0c.GetX()-primVtx.GetX(),
653  d0c.GetY()-primVtx.GetY(),
654  d0c.GetZ()-primVtx.GetZ());
655  Double_t pta = mom.Angle(fline);
656  Double_t cosP = TMath::Cos(pta); // cosine of pta for D0 candidate
657 
658  TVector3 amom(ad0c.GetPx(),ad0c.GetPy(),ad0c.GetPz());
659  TVector3 afline(ad0c.GetX()-aprimVtx.GetX(),
660  ad0c.GetY()-aprimVtx.GetY(),
661  ad0c.GetZ()-aprimVtx.GetZ());
662  Double_t apta = amom.Angle(afline);
663  Double_t acosP = TMath::Cos(apta); // cosine of pta for D0bar candidate
664 
665  // calculate P of Pions at Decay Position of D0 and D0bar candidates
666  negKKF.TransportToParticle(d0c);
667  posPiKF.TransportToParticle(d0c);
668  posKKF.TransportToParticle(ad0c);
669  negPiKF.TransportToParticle(ad0c);
670 
671  Double_t pxPi = posPiKF.GetPx();
672  Double_t pyPi = posPiKF.GetPy();
673  Double_t pzPi = posPiKF.GetPz();
674  Double_t ptPi = posPiKF.GetPt();
675 
676  Double_t apxPi = negPiKF.GetPx();
677  Double_t apyPi = negPiKF.GetPy();
678  Double_t apzPi = negPiKF.GetPz();
679  Double_t aptPi = negPiKF.GetPt();
680 
681  // calculate Pt of Kaons at Decay Position of D0 and D0bar candidates
682 
683  Double_t ptK = negKKF.GetPt();
684  Double_t aptK = posKKF.GetPt();
685 
686  //calculate cos(thetastar)
687  Double_t massvtx = TDatabasePDG::Instance()->GetParticle(421)->Mass();
688  Double_t massp[2];
689  massp[0] = TDatabasePDG::Instance()->GetParticle(321)->Mass();
690  massp[1] = TDatabasePDG::Instance()->GetParticle(211)->Mass();
691  Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)
692  -4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
693 
694  // cos(thetastar) for D0 and Pion
695 
696  Double_t d0E = TMath::Sqrt(massvtx*massvtx + d0P*d0P);
697  Double_t beta = d0P/d0E;
698  Double_t gamma = d0E/massvtx;
699  TVector3 momPi(pxPi,pyPi,pzPi);
700  TVector3 momTot(d0Px,d0Py,d0Pz);
701  Double_t q1 = momPi.Dot(momTot)/momTot.Mag();
702  Double_t cts = (q1/gamma-beta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
703 
704  // cos(thetastar) for D0bar and Pion
705 
706  Double_t ad0E = TMath::Sqrt(massvtx*massvtx + ad0P*ad0P);
707  Double_t abeta = ad0P/ad0E;
708  Double_t agamma = ad0E/massvtx;
709  TVector3 amomPi(apxPi,apyPi,apzPi);
710  TVector3 amomTot(ad0Px,ad0Py,ad0Pz);
711  Double_t aq1 = amomPi.Dot(amomTot)/amomTot.Mag();
712  Double_t acts = (aq1/agamma-abeta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
713 
714  // calculate reduced Chi2 for the full D0 fit
715  d0c.SetProductionVertex(primVtx);
716  ad0c.SetProductionVertex(aprimVtx);
717  negKKF.SetProductionVertex(d0c);
718  posPiKF.SetProductionVertex(d0c);
719  posKKF.SetProductionVertex(ad0c);
720  negPiKF.SetProductionVertex(ad0c);
721  d0c.TransportToProductionVertex();
722  ad0c.TransportToProductionVertex();
723 
724  // calculate the decay length
725  Double_t decayLengthD0 = d0c.GetDecayLength();
726  Double_t adecayLengthD0 = ad0c.GetDecayLength();
727 
728  Double_t chi2D0 = 50.;
729  if(d0c.GetNDF() > 0 && d0c.GetChi2() >= 0) {
730  chi2D0 = d0c.GetChi2()/d0c.GetNDF();
731  }
732 
733  Double_t achi2D0 = 50.;
734  if(ad0c.GetNDF() > 0 && ad0c.GetChi2() >= 0) {
735  achi2D0 = ad0c.GetChi2()/ad0c.GetNDF();
736  }
737 
738  // Get the Pt-bins
739  Int_t ptbin=PtBin(pt);
740  Int_t aptbin=PtBin(apt);
741 
742  if(ptbin < 0) okD0 = 0;
743  if(aptbin < 0) okD0bar = 0;
744  if(!okD0 && !okD0bar) returnvalueCuts=0;
745 
746  if(ptK < fCutsRD[GetGlobalIndex(3,ptbin)] || ptPi < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
747  if(aptK < fCutsRD[GetGlobalIndex(3,aptbin)] || aptPi < fCutsRD[GetGlobalIndex(4,aptbin)]) okD0bar = 0;
748  if(!okD0 && !okD0bar) returnvalueCuts=0;
749 
750 
751  if(TMath::Abs(impactKa) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
752  TMath::Abs(impactPi) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
753 
754  if(TMath::Abs(aimpactKa) > fCutsRD[GetGlobalIndex(5,aptbin)] ||
755  TMath::Abs(aimpactPi) > fCutsRD[GetGlobalIndex(6,aptbin)]) okD0bar = 0;
756 
757  if(!okD0 && !okD0bar) returnvalueCuts=0;
758 
759  // for the moment via the standard method due to bug in AliKF
760  if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) okD0 = 0;
761  if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,aptbin)]) okD0bar = 0;
762  if(!okD0 && !okD0bar) returnvalueCuts=0;
763 
764  if(TMath::Abs(d0mass-massvtx) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
765  if(TMath::Abs(ad0mass-massvtx) > fCutsRD[GetGlobalIndex(0,aptbin)]) okD0bar = 0;
766  if(!okD0 && !okD0bar) returnvalueCuts=0;
767 
768  if(TMath::Abs(cts) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
769  if(TMath::Abs(acts) > fCutsRD[GetGlobalIndex(2,aptbin)]) okD0bar = 0;
770  if(!okD0 && !okD0bar) returnvalueCuts=0;
771 
772  if(prodParam > fCutsRD[GetGlobalIndex(7,ptbin)]) okD0 = 0;
773  if(aprodParam > fCutsRD[GetGlobalIndex(7,aptbin)]) okD0bar = 0;
774  if(!okD0 && !okD0bar) returnvalueCuts=0;
775 
776  if(cosP < fCutsRD[GetGlobalIndex(8,ptbin)]) okD0 = 0;
777  if(acosP < fCutsRD[GetGlobalIndex(8,aptbin)]) okD0bar = 0;
778  if(!okD0 && !okD0bar) returnvalueCuts=0;
779 
780  if(chi2D0 > fCutsRD[GetGlobalIndex(10,ptbin)]) okD0 = 0;
781  if(achi2D0 > fCutsRD[GetGlobalIndex(10,aptbin)]) okD0bar = 0;
782  if(!okD0 && !okD0bar) returnvalueCuts=0;
783 
784  if(decayLengthD0 < fCutsRD[GetGlobalIndex(9,ptbin)]) okD0 = 0;
785  if(adecayLengthD0 < fCutsRD[GetGlobalIndex(9,aptbin)]) okD0bar = 0;
786  if(!okD0 && !okD0bar) returnvalueCuts=0;
787 
788  if(returnvalueCuts!=0) {
789  if(okD0) returnvalueCuts=1; //cuts passed as D0
790  if(okD0bar) returnvalueCuts=2; //cuts passed as D0bar
791  if(okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
792  }
793 
794  return returnvalueCuts;
795 }
796 
797 //---------------------------------------------------------------------------
799 {
800  //
801  // Checking if D0 is in fiducial acceptance region
802  //
803 
804  if(fMaxRapidityCand>-998.){
805  if(TMath::Abs(y) > fMaxRapidityCand) return kFALSE;
806  else return kTRUE;
807  }
808 
809  if(pt > 5.) {
810  // applying cut for pt > 5 GeV
811  AliDebug(2,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt));
812  if (TMath::Abs(y) > 0.8){
813  return kFALSE;
814  }
815  } else {
816  // appliying smooth cut for pt < 5 GeV
817  Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
818  Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
819  AliDebug(2,Form("pt of D0 = %f (< 5), cutting according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
820  if (y < minFiducialY || y > maxFiducialY){
821  return kFALSE;
822  }
823  }
824 
825  return kTRUE;
826 }
827 //---------------------------------------------------------------------------
829 {
830  // ############################################################
831  //
832  // Apply PID selection
833  //
834  //
835  // ############################################################
836 
837  if(!fUsePID) return 3;
838  if(fDefaultPID) return IsSelectedPIDdefault(d);
839  if (fCombPID) return IsSelectedCombPID(d); //to use Bayesian method if applicable
840  fWhyRejection=0;
841  Int_t isD0D0barPID[2]={1,2};
842  Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
843  // same for prong 2
844  // values convention -1 = discarded
845  // 0 = not identified (but compatible) || No PID (->hasPID flag)
846  // 1 = identified
847  // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
848  // Initial hypothesis: unknwon (but compatible)
849  combinedPID[0][0]=0; // prima figlia, Kaon
850  combinedPID[0][1]=0; // prima figlia, pione
851  combinedPID[1][0]=0; // seconda figlia, Kaon
852  combinedPID[1][1]=0; // seconda figlia, pion
853 
854  Bool_t checkPIDInfo[2]={kTRUE,kTRUE};
855  Double_t sigma_tmp[3]={fPidHF->GetSigma(0),fPidHF->GetSigma(1),fPidHF->GetSigma(2)};
856  Bool_t isTOFused=fPidHF->GetTOF(),isCompat=fPidHF->GetCompat();
857  AliAODTrack *aodtrack1=(AliAODTrack*)d->GetDaughter(0);
858  AliAODTrack *aodtrack2=(AliAODTrack*)d->GetDaughter(1);
859  Short_t relativeSign = aodtrack1->Charge() * aodtrack2->Charge();
860  for(Int_t daught=0;daught<2;daught++){
861  //Loop con prongs
862  AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
863  if(fPidHF->IsTOFPiKexcluded(aodtrack,5.)) return 0;
864 
865  if(!(fPidHF->CheckStatus(aodtrack,"TPC")) && !(fPidHF->CheckStatus(aodtrack,"TOF"))) {
866  checkPIDInfo[daught]=kFALSE;
867  continue;
868  }
869  Double_t pProng=aodtrack->P();
870 
871  // identify kaon
872  if(pProng<fmaxPtrackForPID){
873  combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
874  }
875  // identify pion
876  if(pProng<fmaxPtrackForPID){
877  if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
878  combinedPID[daught][1]=0;
879  }else{
880  fPidHF->SetTOF(kFALSE);
881  combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
882  if(isTOFused)fPidHF->SetTOF(kTRUE);
883  if(isCompat)fPidHF->SetCompat(kTRUE);
884  }
885  }
886 
887  if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded
888  isD0D0barPID[0]=0;
889  isD0D0barPID[1]=0;
890  }
891  else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){
892  if((relativeSign == -1 && aodtrack->Charge() == -1) || (relativeSign == 1 && daught == 1)) isD0D0barPID[1]=0;//if K- D0bar excluded
893  else isD0D0barPID[0]=0;// if K+ D0 excluded
894  }
895  /* else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){
896  isD0D0barPID[0]=0;
897  isD0D0barPID[1]=0;
898  }
899  */
900  else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
901  if((relativeSign == -1 && aodtrack->Charge() == -1) || (relativeSign == 1 && daught == 1)) isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
902  else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
903  }
904  else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
905  if((relativeSign == -1 && aodtrack->Charge() == -1) || (relativeSign == 1 && daught == 1)) isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
906  else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
907  }
908 
909  if(fLowPt && d->Pt()<fPtLowPID){
910  Double_t sigmaTPC[3]={3.,2.,0.};
911  fPidHF->SetSigmaForTPC(sigmaTPC);
912  // identify kaon
913  combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
914 
915  if(pProng<0.6){
916  fPidHF->SetCompat(kFALSE);
917  combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
918  fPidHF->SetCompat(kTRUE);
919  }
920 
921  if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
922  combinedPID[daught][1]=0;
923  }else{
924  fPidHF->SetTOF(kFALSE);
925  Double_t sigmaTPCpi[3]={3.,3.,0.};
926  fPidHF->SetSigmaForTPC(sigmaTPCpi);
927  combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
928  fPidHF->SetTOF(kTRUE);
929  if(pProng<0.8){
930  Bool_t isTPCpion=fPidHF->IsPionRaw(aodtrack,"TPC");
931  if(isTPCpion){
932  combinedPID[daught][1]=1;
933  }else{
934  combinedPID[daught][1]=-1;
935  }
936  }
937  }
938 
939  }
940  fPidHF->SetSigmaForTPC(sigma_tmp);
941  }// END OF LOOP ON DAUGHTERS
942 
943  if(!checkPIDInfo[0] && !checkPIDInfo[1]) {
944  if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
945  return 0;
946  }
947 
948 
949  // FURTHER PID REQUEST (both daughter info is needed)
950  if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
951  fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
952  if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
953  return 0;
954  }
955 
956  if(fLowPt && d->Pt()<fPtLowPID){
957  if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
958  fWhyRejection=32;// reject cases where the Kaon is not identified
959  fPidHF->SetSigmaForTPC(sigma_tmp);
960  return 0;
961  }
962  }
963  if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
964 
965  // cout<<"Why? "<<fWhyRejection<<endl;
966  return isD0D0barPID[0]+isD0D0barPID[1];
967 }
968 
969 //---------------------------------------------------------------------------
971 {
972  // ############################################################
973  //
974  // Apply PID selection
975  //
976  //
977  // temporary selection: PID AS USED FOR D0 by Andrea Rossi (up to 28/06/2010)
978  //
979  // d must be a AliAODRecoDecayHF2Prong object
980  // returns 0 if both D0 and D0bar are rejectecd
981  // 1 if D0 is accepted while D0bar is rejected
982  // 2 if D0bar is accepted while D0 is rejected
983  // 3 if both are accepted
984  // fWhyRejection variable (not returned for the moment, print it if needed)
985  // keeps some information on why a candidate has been
986  // rejected according to the following (unfriendly?) scheme
987  // if more rejection cases are considered interesting, just add numbers
988  //
989  // TO BE CONSIDERED WITH A GRAIN OF SALT (the order in which cut are applied is relevant)
990  // from 20 to 30: "detector" selection (PID acceptance)
991  // 26: TPC refit
992  // 27: ITS refit
993  // 28: no (TOF||TPC) pid information (no kTOFpid,kTOFout,kTIME,kTPCpid,...)
994  //
995  // from 30 to 40: PID selection
996  // 31: no Kaon compatible tracks found between daughters
997  // 32: no Kaon identified tracks found (strong sel. at low momenta)
998  // 33: both mass hypotheses are rejected
999  //
1000  // ############################################################
1001 
1002  if(!fUsePID) return 3;
1003  fWhyRejection=0;
1004  Int_t isD0D0barPID[2]={1,2};
1005  Double_t nsigmaTPCpi=-1., nsigmaTPCK=-1.; //used for TPC pid
1006  Double_t tofSig,times[5];// used fot TOF pid
1007  Int_t hasPID[2]={2,2};// flag to count how many detectors give PID info for the daughters
1008  Int_t isKaonPionTOF[2][2],isKaonPionTPC[2][2];
1009  Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
1010  // same for prong 2
1011  // values convention -1 = discarded
1012  // 0 = not identified (but compatible) || No PID (->hasPID flag)
1013  // 1 = identified
1014  // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
1015  // Initial hypothesis: unknwon (but compatible)
1016  isKaonPionTOF[0][0]=0;
1017  isKaonPionTOF[0][1]=0;
1018  isKaonPionTOF[1][0]=0;
1019  isKaonPionTOF[1][1]=0;
1020 
1021  isKaonPionTPC[0][0]=0;
1022  isKaonPionTPC[0][1]=0;
1023  isKaonPionTPC[1][0]=0;
1024  isKaonPionTPC[1][1]=0;
1025 
1026  combinedPID[0][0]=0;
1027  combinedPID[0][1]=0;
1028  combinedPID[1][0]=0;
1029  combinedPID[1][1]=0;
1030 
1031  for(Int_t daught=0;daught<2;daught++){
1032  //Loop con prongs
1033 
1034  // ########### Step 0- CHECKING minimal PID "ACCEPTANCE" ####################
1035 
1036  AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
1037 
1038  if(!(aodtrack->GetStatus()&AliESDtrack::kTPCrefit)){
1039  fWhyRejection=26;
1040  return 0;
1041  }
1042  if(!(aodtrack->GetStatus()&AliESDtrack::kITSrefit)){
1043  fWhyRejection=27;
1044  return 0;
1045  }
1046 
1047  AliAODPid *pid=aodtrack->GetDetPid();
1048  if(!pid) {
1049  //delete esdtrack;
1050  hasPID[daught]--;
1051  continue;
1052  }
1053 
1054  // ########### Step 1- Check of TPC and TOF response ####################
1055 
1056  Double_t ptrack=aodtrack->P();
1057  //#################### TPC PID #######################
1058  if (!(aodtrack->GetStatus()&AliESDtrack::kTPCpid )){
1059  // NO TPC PID INFO FOR THIS TRACK
1060  hasPID[daught]--;
1061  }
1062  else {
1063  static AliTPCPIDResponse theTPCpid;
1064  AliAODPid *pidObj = aodtrack->GetDetPid();
1065  Double_t ptProng=pidObj->GetTPCmomentum();
1066  nsigmaTPCpi = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kPion);
1067  nsigmaTPCK = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kKaon);
1068  //if(ptrack<0.6){
1069  if(ptProng<0.6){
1070  if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
1071  else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1072  if(TMath::Abs(nsigmaTPCpi)<2.)isKaonPionTPC[daught][1]=1;
1073  else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1074  }
1075  //else if(ptrack<.8){
1076  else if(ptProng<.8){
1077  if(TMath::Abs(nsigmaTPCK)<1.)isKaonPionTPC[daught][0]=1;
1078  else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1079  if(TMath::Abs(nsigmaTPCpi)<1.)isKaonPionTPC[daught][1]=1;
1080  else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1081  }
1082  else {
1083  // if(nsigmaTPCK>-2.&&nsigmaTPCK<1.)isKaonPionTPC[daught][0]=1;
1084  if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1085  //if(nsigmaTPCpi>-1.&&nsigmaTPCpi<2.)isKaonPionTPC[daught][1]=1;
1086  if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1087  }
1088  }
1089 
1090  // ##### TOF PID: do not ask nothing for pion/protons ############
1091  if(!((aodtrack->GetStatus()&AliESDtrack::kTOFpid)&&(aodtrack->GetStatus()&AliESDtrack::kTOFout)&&(aodtrack->GetStatus()&AliESDtrack::kTIME))){
1092  // NO TOF PID INFO FOR THIS TRACK
1093  hasPID[daught]--;
1094  }
1095  else{
1096  tofSig=pid->GetTOFsignal();
1097  pid->GetIntegratedTimes(times);
1098  if((tofSig-times[3])>5.*160.)return 0;// PROTON REJECTION
1099  if(TMath::Abs(tofSig-times[3])>3.*160.){
1100  isKaonPionTOF[daught][0]=-1;
1101  }
1102  else {
1103  if(ptrack<1.5){
1104  isKaonPionTOF[daught][0]=1;
1105  }
1106  }
1107  }
1108 
1109  //######### Step 2: COMBINE TOF and TPC PID ###############
1110  // we apply the following convention: if TPC and TOF disagree (discarded Vs identified) -> unknown
1111  combinedPID[daught][0]=isKaonPionTOF[daught][0]+isKaonPionTPC[daught][0];
1112  combinedPID[daught][1]=isKaonPionTOF[daught][1]+isKaonPionTPC[daught][1];
1113 
1114  //######### Step 3: USE PID INFO
1115 
1116  if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){// if not a K- and not a pi- both D0 and D0bar excluded
1117  isD0D0barPID[0]=0;
1118  isD0D0barPID[1]=0;
1119  }
1120  else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){// if in conflict (both pi- and K-), if k for both TPC and TOF -> is K
1121  if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
1122  else isD0D0barPID[0]=0;// if K+ D0 excluded
1123  }
1124  else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){// if in conflict (both pi- and K-) and k- only for TPC or TOF -> reject
1125  isD0D0barPID[0]=0;
1126  isD0D0barPID[1]=0;
1127  }
1128  else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
1129  if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
1130  else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
1131  }
1132  else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
1133  if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
1134  else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
1135  }
1136 
1137  // ########## ALSO DIFFERENT TPC PID REQUEST FOR LOW pt D0: request of K identification ###############################
1138  // ########## more tolerant criteria for single particle ID-> more selective criteria for D0 ##############################
1139  // ############### NOT OPTIMIZED YET ###################################
1140  if(d->Pt()<2.){
1141  isKaonPionTPC[daught][0]=0;
1142  isKaonPionTPC[daught][1]=0;
1143  AliAODPid *pidObj = aodtrack->GetDetPid();
1144  Double_t ptProng=pidObj->GetTPCmomentum();
1145  //if(ptrack<0.6){
1146  if(ptProng<0.6){
1147  if(TMath::Abs(nsigmaTPCK)<3.)isKaonPionTPC[daught][0]=1;
1148  else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1149  if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
1150  else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1151  }
1152  //else if(ptrack<.8){
1153  else if(ptProng<.8){
1154  if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
1155  else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1156  if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
1157  else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1158  }
1159  else {
1160  if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1161  if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1162  }
1163  }
1164 
1165  }// END OF LOOP ON DAUGHTERS
1166 
1167  // FURTHER PID REQUEST (both daughter info is needed)
1168  if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
1169  fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
1170  return 0;
1171  }
1172  else if(hasPID[0]==0&&hasPID[1]==0){
1173  fWhyRejection=28;// reject cases in which no PID info is available
1174  return 0;
1175  }
1176  if(d->Pt()<2.){
1177  // request of K identification at low D0 pt
1178  combinedPID[0][0]=0;
1179  combinedPID[0][1]=0;
1180  combinedPID[1][0]=0;
1181  combinedPID[1][1]=0;
1182 
1183  combinedPID[0][0]=isKaonPionTOF[0][0]+isKaonPionTPC[0][0];
1184  combinedPID[0][1]=isKaonPionTOF[0][1]+isKaonPionTPC[0][1];
1185  combinedPID[1][0]=isKaonPionTOF[1][0]+isKaonPionTPC[1][0];
1186  combinedPID[1][1]=isKaonPionTOF[1][1]+isKaonPionTPC[1][1];
1187 
1188  if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
1189  fWhyRejection=32;// reject cases where the Kaon is not identified
1190  return 0;
1191  }
1192  }
1193 
1194  // cout<<"Why? "<<fWhyRejection<<endl;
1195  return isD0D0barPID[0]+isD0D0barPID[1];
1196 }
1197 
1198 //---------------------------------------------------------------------------
1200  Int_t selectionvalCand,
1201  Int_t selectionvalPID) const
1202 {
1203  //
1204  // This method combines the tracks, PID and cuts selection results
1205  //
1206  if(selectionvalTrack==0) return 0;
1207 
1208  Int_t returnvalue;
1209 
1210  switch(selectionvalPID) {
1211  case 0:
1212  returnvalue=0;
1213  break;
1214  case 1:
1215  returnvalue=((selectionvalCand==1 || selectionvalCand==3) ? 1 : 0);
1216  break;
1217  case 2:
1218  returnvalue=((selectionvalCand==2 || selectionvalCand==3) ? 2 : 0);
1219  break;
1220  case 3:
1221  returnvalue=selectionvalCand;
1222  break;
1223  default:
1224  returnvalue=0;
1225  break;
1226  }
1227 
1228  return returnvalue;
1229 }
1230 
1231 //----------------------------------------------------------------------------
1233 {
1234  //
1235  // Note: this method is temporary
1236  // Additional cuts on decay lenght and lower cut for d0 norm are applied using vertex without candidate's daughters
1237  //
1238 
1239  //apply cuts
1240 
1241  Float_t normDecLengthCut=1.,decLengthCut=TMath::Min(d->P()*0.0066+0.01,0.06/*cm*/), normd0Cut=0.5;
1242  // "decay length" expo law with tau' = beta*gamma*ctau= p/m*ctau =p*0.0123/1.864~p*0.0066
1243  // decay lenght > ctau' implies to retain (1-1/e) (for signal without considering detector resolution),
1244 
1245  Int_t returnvalue=3; //cut passed
1246  for(Int_t i=0;i<2/*prongs*/;i++){
1247  if(TMath::Abs(d->Normalizedd0Prong(i))<normd0Cut) return 0; //normd0Cut not passed
1248  }
1249  if(d->DecayLength2()<decLengthCut*decLengthCut) return 0; //decLengthCut not passed
1250  if(d->NormalizedDecayLength2()<normDecLengthCut*normDecLengthCut) return 0; //decLengthCut not passed
1251 
1252  return returnvalue;
1253 }
1254 
1255 //----------------------------------------------
1257 {
1258  //switch on candidate selection via AliKFparticle
1259  if(!useKF) return;
1260  if(useKF){
1261  fUseKF=useKF;
1262  Int_t nvarsKF=11;
1263  SetNVars(nvarsKF);
1264  TString varNamesKF[11]={"inv. mass [GeV]",
1265  "dca [cm]",
1266  "cosThetaStar",
1267  "pTK [GeV/c]",
1268  "pTPi [GeV/c]",
1269  "d0K [cm]",
1270  "d0Pi [cm]",
1271  "d0d0 [cm^2]",
1272  "cosThetaPoint"
1273  "DecayLength[cm]",
1274  "RedChi2"};
1275  Bool_t isUpperCutKF[11]={kTRUE,
1276  kTRUE,
1277  kTRUE,
1278  kFALSE,
1279  kFALSE,
1280  kTRUE,
1281  kTRUE,
1282  kTRUE,
1283  kFALSE,
1284  kFALSE,
1285  kTRUE};
1286  SetVarNames(nvarsKF,varNamesKF,isUpperCutKF);
1287  SetGlobalIndex();
1288  Bool_t forOpt[11]={kFALSE,
1289  kTRUE,
1290  kTRUE,
1291  kFALSE,
1292  kFALSE,
1293  kFALSE,
1294  kFALSE,
1295  kTRUE,
1296  kTRUE,
1297  kFALSE,
1298  kFALSE};
1299  SetVarsForOpt(4,forOpt);
1300  }
1301  return;
1302 }
1303 
1304 //---------------------------------------------------------------------------
1306  //
1307  //STANDARD CUTS USED FOR 2010 pp analysis
1308  //dca cut will be enlarged soon to 400 micron
1309  //
1310 
1311  SetName("D0toKpiCutsStandard");
1312  SetTitle("Standard Cuts for D0 analysis");
1313 
1314  // PILE UP REJECTION
1316 
1317  // EVENT CUTS
1318  SetMinVtxContr(1);
1319  // MAX Z-VERTEX CUT
1320  SetMaxVtxZ(10.);
1321 
1322  // TRACKS ON SINGLE TRACKS
1323  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1324  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1325  esdTrackCuts->SetRequireTPCRefit(kTRUE);
1326  esdTrackCuts->SetRequireITSRefit(kTRUE);
1327  // esdTrackCuts->SetMinNClustersITS(4);
1328  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1329  esdTrackCuts->SetMinDCAToVertexXY(0.);
1330  esdTrackCuts->SetEtaRange(-0.8,0.8);
1331  esdTrackCuts->SetPtRange(0.3,1.e10);
1332 
1333  AddTrackCuts(esdTrackCuts);
1334  delete esdTrackCuts;
1335  esdTrackCuts=NULL;
1336 
1337  const Int_t nptbins =14;
1338  const Double_t ptmax = 9999.;
1339  const Int_t nvars=11;
1340  Float_t ptbins[nptbins+1];
1341  ptbins[0]=0.;
1342  ptbins[1]=0.5;
1343  ptbins[2]=1.;
1344  ptbins[3]=2.;
1345  ptbins[4]=3.;
1346  ptbins[5]=4.;
1347  ptbins[6]=5.;
1348  ptbins[7]=6.;
1349  ptbins[8]=7.;
1350  ptbins[9]=8.;
1351  ptbins[10]=12.;
1352  ptbins[11]=16.;
1353  ptbins[12]=20.;
1354  ptbins[13]=24.;
1355  ptbins[14]=ptmax;
1356 
1357  SetGlobalIndex(nvars,nptbins);
1358  SetPtBins(nptbins+1,ptbins);
1359 
1360  Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,350.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,0.},/* pt<0.5*/
1361  {0.400,350.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,0.},/* 0.5<pt<1*/
1362  {0.400,300.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.80,0.,0.},/* 1<pt<2 */
1363  {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-15000.*1E-8,0.85,0.,0.},/* 2<pt<3 */
1364  {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 3<pt<4 */
1365  {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 4<pt<5 */
1366  {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 5<pt<6 */
1367  {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 6<pt<7 */
1368  {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-7000.*1E-8,0.85,0.,0.},/* 7<pt<8 */
1369  {0.400,300.*1E-4,0.9,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.,0.},/* 8<pt<12 */
1370  {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,10000.*1E-8,0.85,0.,0.},/* 12<pt<16 */
1371  {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,999999.*1E-8,0.85,0.,0.},/* 16<pt<20 */
1372  {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,999999.*1E-8,0.85,0.,0.},/* 20<pt<24 */
1373  {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,999999.*1E-8,0.85,0.,0.}};/* pt>24 */
1374 
1375  //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1376  Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1377  for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1378 
1379  for (Int_t ibin=0;ibin<nptbins;ibin++){
1380  for (Int_t ivar = 0; ivar<nvars; ivar++){
1381  cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1382  }
1383  }
1384 
1385  SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1386  SetUseSpecialCuts(kTRUE);
1388 
1389  for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1390  delete [] cutsMatrixTransposeStand;
1391  cutsMatrixTransposeStand=NULL;
1392 
1393  // PID SETTINGS
1394  AliAODPidHF* pidObj=new AliAODPidHF();
1395  //pidObj->SetName("pid4D0");
1396  Int_t mode=1;
1397  const Int_t nlims=2;
1398  Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1399  Bool_t compat=kTRUE; //effective only for this mode
1400  Bool_t asym=kTRUE;
1401  Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1402  pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1403  pidObj->SetMatch(mode);
1404  pidObj->SetPLimit(plims,nlims);
1405  pidObj->SetSigma(sigmas);
1406  pidObj->SetCompat(compat);
1407  pidObj->SetTPC(kTRUE);
1408  pidObj->SetTOF(kTRUE);
1409  pidObj->SetPCompatTOF(1.5);
1410  pidObj->SetSigmaForTPCCompat(3.);
1411  pidObj->SetSigmaForTOFCompat(3.);
1412  pidObj->SetOldPid(kTRUE);
1413 
1414  SetPidHF(pidObj);
1415  SetUsePID(kTRUE);
1416  SetUseDefaultPID(kFALSE);
1417 
1418  SetLowPt(kFALSE);
1419 
1420  PrintAll();
1421 
1422  delete pidObj;
1423  pidObj=NULL;
1424 
1425  return;
1426 }
1427 
1428 //___________________________________________________________________________
1430  //
1431  // Cuts for 2010 pp 7 TeV data analysis in Ntracklets bins (vs multiplicity)
1432  //
1433  SetName("D0toKpiCuts");
1434  SetTitle("Cuts for D0 analysis in 2010-data pp 7 TeV vs multiplicity");
1435 
1436  //
1437  // Track cuts
1438  //
1439  AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
1440  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1441  //default
1442  esdTrackCuts->SetRequireTPCRefit(kTRUE);
1443  esdTrackCuts->SetRequireITSRefit(kTRUE);
1444  esdTrackCuts->SetEtaRange(-0.8,0.8);
1445  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1446  AliESDtrackCuts::kAny);
1447  // default is kBoth, otherwise kAny
1448  esdTrackCuts->SetMinDCAToVertexXY(0.);
1449  esdTrackCuts->SetPtRange(0.3,1.e10);
1450 
1451  AddTrackCuts(esdTrackCuts);
1452  delete esdTrackCuts;
1453  esdTrackCuts=NULL;
1454 
1455  //
1456  // Cut values per pt bin
1457  //
1458  const Int_t nvars=11;
1459  const Int_t nptbins=14;
1460  Float_t* ptbins;
1461  ptbins=new Float_t[nptbins+1];
1462  ptbins[0]=0.;
1463  ptbins[1]=0.5;
1464  ptbins[2]=1.;
1465  ptbins[3]=2.;
1466  ptbins[4]=3.;
1467  ptbins[5]=4.;
1468  ptbins[6]=5.;
1469  ptbins[7]=6.;
1470  ptbins[8]=7.;
1471  ptbins[9]=8.;
1472  ptbins[10]=12.;
1473  ptbins[11]=16.;
1474  ptbins[12]=20.;
1475  ptbins[13]=24.;
1476  ptbins[14]=9999.;
1477 
1478  SetPtBins(nptbins+1,ptbins);
1479 
1480  //setting cut values
1481  Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,350.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,3.},/* pt<0.5*/
1482  {0.400,350.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,3.},/* 0.5<pt<1*/
1483  {0.400,300.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-30000.*1E-8,0.80,0.,4.},/* 1<pt<2 */
1484  {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-20000.*1E-8,0.85,0.,4.},/* 2<pt<3 */
1485  {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-15000.*1E-8,0.85,0.,4.},/* 3<pt<4 */
1486  {0.400,300.*1E-4,0.75,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.875,0.,0.},/* 4<pt<5 */
1487  {0.400,300.*1E-4,0.75,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.875,0.,0.},/* 5<pt<6 */
1488  {0.400,400.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 6<pt<7 */
1489  {0.400,400.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 7<pt<8 */
1490  {0.400,0.06,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.00001,0.85,0.,0.},/* 8<pt<12 */
1491  {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 12<pt<16 */
1492  {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 16<pt<20 */
1493  {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 20<pt<24 */
1494  {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.}};/* pt>24 */
1495 
1496  //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1497  Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1498  for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1499 
1500  for (Int_t ibin=0;ibin<nptbins;ibin++){
1501  for (Int_t ivar = 0; ivar<nvars; ivar++){
1502  cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1503  }
1504  }
1505 
1506  SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1507  for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1508  delete [] cutsMatrixTransposeStand;
1509 
1510  SetUseSpecialCuts(kTRUE);
1512 
1513  //Do recalculate the vertex
1515 
1516  //
1517  // Pid settings
1518  //
1519  Bool_t pidflag=kTRUE;
1520  SetUsePID(pidflag);
1521  if(pidflag) cout<<"PID is used"<<endl;
1522  else cout<<"PID is not used"<<endl;
1523  //
1524  AliAODPidHF* pidObj=new AliAODPidHF();
1525  Int_t mode=1;
1526  const Int_t nlims=2;
1527  Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1528  Bool_t compat=kTRUE; //effective only for this mode
1529  Bool_t asym=kTRUE;
1530  Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1531  pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1532  pidObj->SetMatch(mode);
1533  pidObj->SetPLimit(plims,nlims);
1534  pidObj->SetSigma(sigmas);
1535  pidObj->SetCompat(compat);
1536  pidObj->SetTPC(kTRUE);
1537  pidObj->SetTOF(kTRUE);
1538  pidObj->SetPCompatTOF(1.5);
1539  pidObj->SetSigmaForTPCCompat(3.);
1540  pidObj->SetSigmaForTOFCompat(3.);
1541  pidObj->SetOldPid(kTRUE);
1542  // pidObj->SetOldPid(kFALSE);
1543  SetPidHF(pidObj);
1544 
1545  SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
1546 
1547  SetLowPt(kFALSE);
1548 
1549  //activate pileup rejection (for pp)
1551 
1552  PrintAll();
1553 
1554  delete pidObj;
1555  pidObj=NULL;
1556  delete [] ptbins;
1557  ptbins=NULL;
1558 
1559  return;
1560 }
1561 
1562 //---------------------------------------------------------------------------
1564  //
1565  // STANDARD CUTS USED FOR 2011 pp analysis at 2.76TeV
1566  //
1567 
1568  SetName("D0toKpiCutsStandard");
1569  SetTitle("Standard Cuts for D0 analysis in pp2011 at 2.76TeV run");
1570 
1571  //
1572  // Track cuts
1573  //
1574  AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
1575  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1576  //default
1577  esdTrackCuts->SetRequireTPCRefit(kTRUE);
1578  esdTrackCuts->SetRequireITSRefit(kTRUE);
1579  esdTrackCuts->SetEtaRange(-0.8,0.8);
1580  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1581  AliESDtrackCuts::kAny);
1582  // default is kBoth, otherwise kAny
1583  esdTrackCuts->SetMinDCAToVertexXY(0.);
1584  esdTrackCuts->SetPtRange(0.3,1.e10);
1585 
1586  esdTrackCuts->SetMaxDCAToVertexXY(1.);
1587  esdTrackCuts->SetMaxDCAToVertexZ(1.);
1588  esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1589 
1590  AddTrackCuts(esdTrackCuts);
1591  delete esdTrackCuts;
1592  esdTrackCuts=NULL;
1593 
1594  const Int_t nvars=11;
1595  const Int_t nptbins=13;
1596  Float_t ptbins[nptbins+1];
1597  ptbins[0]=0.;
1598  ptbins[1]=0.5;
1599  ptbins[2]=1.;
1600  ptbins[3]=2.;
1601  ptbins[4]=3.;
1602  ptbins[5]=4.;
1603  ptbins[6]=5.;
1604  ptbins[7]=6.;
1605  ptbins[8]=8.;
1606  ptbins[9]=12.;
1607  ptbins[10]=16.;
1608  ptbins[11]=20.;
1609  ptbins[12]=24.;
1610  ptbins[13]=9999.;
1611 
1612  SetPtBins(nptbins+1,ptbins);
1613 
1614  Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,0.04,0.75,0.3,0.3,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* pt<0.5*/
1615  {0.400,0.04,0.75,0.3,0.3,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* 0.5<pt<1*/
1616  {0.400,0.03,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.8,0.,0.},/* 1<pt<2 */
1617  {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.0003,0.9,0.,0.},/* 2<pt<3 */
1618  {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.0002,0.9,0.,0.},/* 3<pt<4 */
1619  {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.00015,0.9,0.,0.},/* 4<pt<5 */
1620  {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.0001,0.9,0.,0.},/* 5<pt<6 */
1621  {0.400,0.09,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* 6<pt<8 */
1622  {0.400,0.06,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.00001,0.85,0.,0.},/* 8<pt<12 */
1623  {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 12<pt<16 */
1624  {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 16<pt<20 */
1625  {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 20<pt<24 */
1626  {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.}};/* pt>24 */
1627 
1628  //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1629  Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1630  for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1631  for (Int_t ibin=0;ibin<nptbins;ibin++){
1632  for (Int_t ivar = 0; ivar<nvars; ivar++){
1633  cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1634  }
1635  }
1636  SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1637  for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1638  delete [] cutsMatrixTransposeStand;
1639 
1640  //pid settings
1641  AliAODPidHF* pidObj=new AliAODPidHF();
1642  Int_t mode=1;
1643  const Int_t nlims=2;
1644  Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1645  Bool_t compat=kTRUE; //effective only for this mode
1646  Bool_t asym=kTRUE;
1647  Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1648  pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1649  pidObj->SetMatch(mode);
1650  pidObj->SetPLimit(plims,nlims);
1651  pidObj->SetSigma(sigmas);
1652  pidObj->SetCompat(compat);
1653  pidObj->SetTPC(kTRUE);
1654  pidObj->SetTOF(kTRUE);
1655  pidObj->SetOldPid(kTRUE);
1656  SetPidHF(pidObj);
1657 
1658  SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
1659 
1660  SetUsePID(kTRUE);
1661 
1662  //activate pileup rejection (for pp)
1664 
1665  //Do recalculate the vertex
1667 
1668  // Cut on the zvtx
1669  SetMaxVtxZ(10.);
1670 
1671  // Use the kFirst selection for tracks with candidate D with pt<2
1672  SetSelectCandTrackSPDFirst(kTRUE,2.);
1673 
1674  // Use special cuts for D candidates with pt<2
1675  SetUseSpecialCuts(kTRUE);
1677 
1678  PrintAll();
1679 
1680  delete pidObj;
1681  pidObj=NULL;
1682 
1683  return;
1684 }
1685 
1686 //---------------------------------------------------------------------------
1688  //
1689  // Final CUTS USED FOR 2010 PbPb analysis of central events
1690  // REMEMBER TO EVENTUALLY SWITCH ON
1691  // SetUseAOD049(kTRUE);
1692  //
1693 
1694  SetName("D0toKpiCutsStandard");
1695  SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
1696 
1697  // CENTRALITY SELECTION
1698  SetMinCentrality(0.);
1699  SetMaxCentrality(80.);
1701 
1702  // EVENT CUTS
1703  SetMinVtxContr(1);
1704  // MAX Z-VERTEX CUT
1705  SetMaxVtxZ(10.);
1706  // PILE UP
1708  // PILE UP REJECTION
1709  //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1710 
1711  // TRACKS ON SINGLE TRACKS
1712  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1713  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1714  esdTrackCuts->SetRequireTPCRefit(kTRUE);
1715  esdTrackCuts->SetRequireITSRefit(kTRUE);
1716  // esdTrackCuts->SetMinNClustersITS(4);
1717  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1718  esdTrackCuts->SetMinDCAToVertexXY(0.);
1719  esdTrackCuts->SetEtaRange(-0.8,0.8);
1720  esdTrackCuts->SetPtRange(0.7,1.e10);
1721 
1722  esdTrackCuts->SetMaxDCAToVertexXY(1.);
1723  esdTrackCuts->SetMaxDCAToVertexZ(1.);
1724  esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1725 
1726  AddTrackCuts(esdTrackCuts);
1727  delete esdTrackCuts;
1728  esdTrackCuts=NULL;
1729  // SPD k FIRST for Pt<3 GeV/c
1730  SetSelectCandTrackSPDFirst(kTRUE, 3);
1731 
1732  // CANDIDATE CUTS
1733  const Int_t nptbins =13;
1734  const Double_t ptmax = 9999.;
1735  const Int_t nvars=11;
1736  Float_t ptbins[nptbins+1];
1737  ptbins[0]=0.;
1738  ptbins[1]=0.5;
1739  ptbins[2]=1.;
1740  ptbins[3]=2.;
1741  ptbins[4]=3.;
1742  ptbins[5]=4.;
1743  ptbins[6]=5.;
1744  ptbins[7]=6.;
1745  ptbins[8]=8.;
1746  ptbins[9]=12.;
1747  ptbins[10]=16.;
1748  ptbins[11]=20.;
1749  ptbins[12]=24.;
1750  ptbins[13]=ptmax;
1751 
1752  SetGlobalIndex(nvars,nptbins);
1753  SetPtBins(nptbins+1,ptbins);
1754  SetMinPtCandidate(2.);
1755 
1756  Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.85,0.99,2.},/* pt<0.5*/
1757  {0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.9,0.99,2.},/* 0.5<pt<1*/
1758  {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-43000.*1E-8,0.85,0.995,8.},/* 1<pt<2 */
1759  {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-45000.*1E-8,0.95,0.998,7.},/* 2<pt<3 */
1760  {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-36000.*1E-8,0.95,0.998,5.},/* 3<pt<4 */
1761  {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-27000.*1E-8,0.95,0.998,5.},/* 4<pt<5 */
1762  {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-21000.*1E-8,0.92,0.998,5.},/* 5<pt<6 */
1763  {0.400,270.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-14000.*1E-8,0.88,0.998,5.},/* 6<pt<8 */
1764  {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.998,5.},/* 8<pt<12 */
1765  {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.83,0.998,8.},/* 12<pt<16 */
1766  {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.82,0.995,6.},/* 16<pt<20 */
1767  {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.81,0.995,6.},/* 20<pt<24 */
1768  {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.8,0.99,2.}};/* pt>24 */
1769 
1770  //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1771  Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1772  for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1773 
1774  for (Int_t ibin=0;ibin<nptbins;ibin++){
1775  for (Int_t ivar = 0; ivar<nvars; ivar++){
1776  cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1777  }
1778  }
1779 
1780  SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1781  SetUseSpecialCuts(kTRUE);
1782  SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1783  for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1784  delete [] cutsMatrixTransposeStand;
1785  cutsMatrixTransposeStand=NULL;
1786 
1787  // PID SETTINGS
1788  AliAODPidHF* pidObj=new AliAODPidHF();
1789  //pidObj->SetName("pid4D0");
1790  Int_t mode=1;
1791  const Int_t nlims=2;
1792  Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1793  Bool_t compat=kTRUE; //effective only for this mode
1794  Bool_t asym=kTRUE;
1795  Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1796  pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1797  pidObj->SetMatch(mode);
1798  pidObj->SetPLimit(plims,nlims);
1799  pidObj->SetSigma(sigmas);
1800  pidObj->SetCompat(compat);
1801  pidObj->SetTPC(kTRUE);
1802  pidObj->SetTOF(kTRUE);
1803  pidObj->SetPCompatTOF(2.);
1804  pidObj->SetSigmaForTPCCompat(3.);
1805  pidObj->SetSigmaForTOFCompat(3.);
1806  pidObj->SetOldPid(kTRUE);
1807 
1808  SetPidHF(pidObj);
1809  SetUsePID(kTRUE);
1810  SetUseDefaultPID(kFALSE);
1811 
1812  PrintAll();
1813 
1814  delete pidObj;
1815  pidObj=NULL;
1816 
1817  return;
1818 }
1819 
1820 //---------------------------------------------------------------------------
1822  // CUTS USED FOR D0 RAA for the CENTRALITY RANGE 20-80%
1823 
1824  SetName("D0toKpiCutsStandard");
1825  SetTitle("Standard Cuts for D0 analysis in PbPb2010 run, for peripheral events");
1826 
1827  // CENTRALITY SELECTION
1828  SetMinCentrality(40.);
1829  SetMaxCentrality(80.);
1831 
1832  // EVENT CUTS
1833  SetMinVtxContr(1);
1834  // MAX Z-VERTEX CUT
1835  SetMaxVtxZ(10.);
1836  // PILE UP
1838  // PILE UP REJECTION
1839  //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1840 
1841  // TRACKS ON SINGLE TRACKS
1842  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1843  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1844  esdTrackCuts->SetRequireTPCRefit(kTRUE);
1845  esdTrackCuts->SetRequireITSRefit(kTRUE);
1846  // esdTrackCuts->SetMinNClustersITS(4);
1847  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1848  esdTrackCuts->SetMinDCAToVertexXY(0.);
1849  esdTrackCuts->SetEtaRange(-0.8,0.8);
1850  esdTrackCuts->SetPtRange(0.5,1.e10);
1851 
1852  esdTrackCuts->SetMaxDCAToVertexXY(1.);
1853  esdTrackCuts->SetMaxDCAToVertexZ(1.);
1854  esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0025*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1855 
1856  AddTrackCuts(esdTrackCuts);
1857  delete esdTrackCuts;
1858  esdTrackCuts=NULL;
1859  // SPD k FIRST for Pt<3 GeV/c
1860  SetSelectCandTrackSPDFirst(kTRUE, 3);
1861 
1862  // CANDIDATE CUTS
1863  const Int_t nptbins =13;
1864  const Double_t ptmax = 9999.;
1865  const Int_t nvars=11;
1866  Float_t ptbins[nptbins+1];
1867  ptbins[0]=0.;
1868  ptbins[1]=0.5;
1869  ptbins[2]=1.;
1870  ptbins[3]=2.;
1871  ptbins[4]=3.;
1872  ptbins[5]=4.;
1873  ptbins[6]=5.;
1874  ptbins[7]=6.;
1875  ptbins[8]=8.;
1876  ptbins[9]=12.;
1877  ptbins[10]=16.;
1878  ptbins[11]=20.;
1879  ptbins[12]=24.;
1880  ptbins[13]=ptmax;
1881 
1882  SetGlobalIndex(nvars,nptbins);
1883  SetPtBins(nptbins+1,ptbins);
1884  SetMinPtCandidate(0.);
1885 
1886  Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-20000.*1E-8,0.7,0.993,2.},/* pt<0.5*/
1887  {0.400,400.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.85,0.993,2.},/* 0.5<pt<1*/
1888  {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.85,0.995,6.},/* 1<pt<2 */
1889  {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.95,0.991,5.},/* 2<pt<3 */
1890  {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-36000.*1E-8,0.95,0.993,5.},/* 3<pt<4 */
1891  {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-27000.*1E-8,0.95,0.995,5.},/* 4<pt<5 */
1892  {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-21000.*1E-8,0.92,0.998,5.},/* 5<pt<6 */
1893  {0.400,270.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-14000.*1E-8,0.88,0.998,5.},/* 6<pt<8 */
1894  {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.995,5.},/* 8<pt<12 */
1895  {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.83,0.995,4.},/* 12<pt<16 */
1896  {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.82,0.995,4.},/* 16<pt<20 */
1897  {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.81,0.995,4.},/* 20<pt<24 */
1898  {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.8,0.995,4.}};/* pt>24 */
1899 
1900  //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1901  Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1902  for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1903 
1904  for (Int_t ibin=0;ibin<nptbins;ibin++){
1905  for (Int_t ivar = 0; ivar<nvars; ivar++){
1906  cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1907  }
1908  }
1909 
1910  SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1911  SetUseSpecialCuts(kTRUE);
1912  SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1913  for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1914  delete [] cutsMatrixTransposeStand;
1915  cutsMatrixTransposeStand=NULL;
1916 
1917  // PID SETTINGS
1918  AliAODPidHF* pidObj=new AliAODPidHF();
1919  //pidObj->SetName("pid4D0");
1920  Int_t mode=1;
1921  const Int_t nlims=2;
1922  Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1923  Bool_t compat=kTRUE; //effective only for this mode
1924  Bool_t asym=kTRUE;
1925  Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1926  pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1927  pidObj->SetMatch(mode);
1928  pidObj->SetPLimit(plims,nlims);
1929  pidObj->SetSigma(sigmas);
1930  pidObj->SetCompat(compat);
1931  pidObj->SetTPC(kTRUE);
1932  pidObj->SetTOF(kTRUE);
1933  pidObj->SetPCompatTOF(2.);
1934  pidObj->SetSigmaForTPCCompat(3.);
1935  pidObj->SetSigmaForTOFCompat(3.);
1936  pidObj->SetOldPid(kTRUE);
1937 
1938  SetPidHF(pidObj);
1939  SetUsePID(kTRUE);
1940  SetUseDefaultPID(kFALSE);
1941 
1942  SetLowPt(kTRUE,2.);
1943 
1944  PrintAll();
1945 
1946  delete pidObj;
1947  pidObj=NULL;
1948 
1949  return;
1950 }
1951 
1952 
1953 //---------------------------------------------------------------------------
1955 {
1956  // Default 2010 PbPb cut object
1958  AliAODPidHF *pidobj=GetPidHF();
1959 
1960  pidobj->SetOldPid(kFALSE);
1961 
1962  //
1963  // Enable all 2011 PbPb run triggers
1964  //
1965  SetTriggerClass("");
1969 
1970 }
1971 
1972 //---------------------------------------------------------------------------
1974 {
1975  // ############################################################
1976  //
1977  // Apply Bayesian PID selection
1978  // Implementation of Bayesian PID: Jeremy Wilkinson
1979  //
1980  // ############################################################
1981 
1982  if (!fUsePID || !d) return 3;
1983 
1984 
1986  //WeightNoFilter: Accept all particles (no PID cut) but fill mass histos with weights in task
1988  return 3;
1989  }
1990 
1991  Int_t isD0 = 0;
1992  Int_t isD0bar = 0;
1993  Int_t returnvalue = 0;
1994 
1995  Int_t isPosKaon = 0, isNegKaon = 0, isPosPion = 0, isNegPion = 0;
1996 
1997  //Bayesian methods used here check for ID of kaon, and whether it is positive or negative.
1998 
1999  Bool_t checkPIDInfo[2] = {kTRUE, kTRUE};
2000  AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
2001 
2002  if ((aodtrack[0]->Charge()*aodtrack[1]->Charge()) != -1) {
2003  return 0; //reject if charges not opposite
2004  }
2005  Double_t momentumpositive=0., momentumnegative=0.; //Used in "prob > prior" method
2006  for (Int_t daught = 0; daught < 2; daught++) {
2007  //Loop over prongs
2008 
2009  if (aodtrack[daught]->Charge() == -1) {
2010  momentumnegative = aodtrack[daught]->P();
2011  }
2012  if (aodtrack[daught]->Charge() == +1) {
2013  momentumpositive = aodtrack[daught]->P();
2014  }
2015  if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
2016  checkPIDInfo[daught] = kFALSE;
2017  continue;
2018  }
2019 
2020  }
2021 
2022  //Loop over daughters ends here
2023 
2024  if (!checkPIDInfo[0] && !checkPIDInfo[1]) {
2025  return 0; //Reject if both daughters lack both TPC and TOF info
2026  }
2027 
2028  CalculateBayesianWeights(d); //Calculates all Bayesian probabilities for both positive and negative tracks
2029  // Double_t prob0[AliPID::kSPECIES];
2030  // fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), prob0);
2032  switch (fBayesianCondition) {
2034  case kMaxProb:
2035  if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kKaon]) { //If highest probability lies with kaon
2036  isPosKaon = 1; //flag [daught] as a kaon
2037  }
2038 
2039  if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kPion]) { //If highest probability lies with pion
2040  isPosPion = 1; //flag [daught] as a pion
2041  }
2042 
2043  if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kKaon]) { //If highest probability lies with kaon
2044  isNegKaon = 1; //flag [daught] as a kaon
2045  }
2046 
2047  if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kPion]) { //If highest probability lies with kaon
2048  isNegPion = 1; //flag [daught] as a pion
2049  }
2050 
2051  break;
2053  case kAbovePrior:
2054 
2055  if (fWeightsNegative[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
2056  GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumnegative)))) { //Retrieves relevant prior, gets value at momentum
2057  isNegKaon = 1;
2058  }
2059  if (fWeightsPositive[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
2060  GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumpositive)))) { //Retrieves relevant prior, gets value at momentum
2061  isPosKaon = 1;
2062  }
2063 
2064  break;
2065 
2067  case kThreshold:
2068  if (fWeightsNegative[AliPID::kKaon] > fProbThreshold) {
2069  isNegKaon = 1;
2070  }
2071  if (fWeightsNegative[AliPID::kPion] > fProbThreshold) {
2072  isNegPion = 1;
2073  }
2074 
2075  if (fWeightsPositive[AliPID::kKaon] > fProbThreshold) {
2076  isPosKaon = 1;
2077  }
2078 
2079  if (fWeightsPositive[AliPID::kPion] > fProbThreshold) {
2080  isPosPion = 1;
2081  }
2082 
2083  break;
2084  }
2085 
2086  //Momentum-based selection (also applied to filtered weighted method)
2087 
2089  if (isNegKaon && isPosKaon) { // If both are kaons, reject
2090  isD0 = 1;
2091  isD0bar = 1;
2092  } else if (isNegKaon && isPosPion) { //If negative kaon present, D0
2093  isD0 = 1;
2094  } else if (isPosKaon && isNegPion) { //If positive kaon present, D0bar
2095  isD0bar = 1;
2096  } else { //If neither ID'd as kaon, subject to extra tests
2097  isD0 = 1;
2098  isD0bar = 1;
2099  if (aodtrack[0]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[0], "TOF"))) { //If it's low momentum and not ID'd as a kaon, we assume it must be a pion
2100  if (aodtrack[0]->Charge() == -1) {
2101  isD0 = 0; //Check charge and reject based on pion hypothesis
2102  }
2103  if (aodtrack[0]->Charge() == 1) {
2104  isD0bar = 0;
2105  }
2106  }
2107  if (aodtrack[1]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[1], "TOF"))) {
2108  if (aodtrack[1]->Charge() == -1) {
2109  isD0 = 0;
2110  }
2111  if (aodtrack[1]->Charge() == 1) {
2112  isD0bar = 0;
2113  }
2114  }
2115  }
2116 
2117  if (isD0 && isD0bar) {
2118  returnvalue = 3;
2119  }
2120  if (isD0 && !isD0bar) {
2121  returnvalue = 1;
2122  }
2123  if (!isD0 && isD0bar) {
2124  returnvalue = 2;
2125  }
2126  if (!isD0 && !isD0bar) {
2127  returnvalue = 0;
2128  }
2129  }
2130 
2131  //Simple Bayesian
2133 
2134  if (isPosKaon && isNegKaon) { //If both are ID'd as kaons, accept as possible
2135  returnvalue = 3;
2136  } else if (isNegKaon && isPosPion) { //If negative kaon, D0
2137  returnvalue = 1;
2138  } else if (isPosKaon && isNegPion) { //If positive kaon, D0-bar
2139  returnvalue = 2;
2140  } else if (isPosPion && isNegPion) {
2141  returnvalue = 0; //If neither kaon, reject
2142  } else {returnvalue = 0;} //default
2143 
2144  }
2145 
2146  return returnvalue;
2147 }
2148 
2149 //---------------------------------------------------------------------------
2151 {
2152  //Function to compute weights for Bayesian method
2153 
2154  AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
2155  if ((aodtrack[0]->Charge() * aodtrack[1]->Charge()) != -1) {
2156  return; //Reject if charges do not oppose
2157  }
2158  for (Int_t daught = 0; daught < 2; daught++) {
2159  //Loop over prongs
2160 
2161  if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
2162  continue;
2163  }
2164 
2165  // identify kaon, define weights
2166  if (aodtrack[daught]->Charge() == +1) {
2167  fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsPositive);
2168  }
2169 
2170  if (aodtrack[daught]->Charge() == -1) {
2171  fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsNegative);
2172  }
2173  }
2174 }
2175 //___________________________________________________________
2177  //
2178  // print all cuts values
2179  //
2182  printf("Cuts on d0meas-d0exp:\n");
2183  for(Int_t ib=0;ib<fUsed0MeasMinusExpCut;ib++){
2184  printf("%f ",fMaxd0MeasMinusExp[ib]);
2185  }
2186  printf("\n");
2187  }
2188  if(fUseImpParDCut){
2189  printf("Cuts on d0:\n");
2190  for(Int_t ib=0;ib<fUseImpParDCut;ib++){
2191  printf("%f ",fMaxImpParD[ib]);
2192  }
2193  printf("\n");
2194  }
2195 }
Double_t NormalizedDecayLengthXY() const
Int_t fIsSelectedCuts
fix the daughter track references
Definition: AliRDHFCuts.h:422
void SetAsym(Bool_t asym)
Definition: AliAODPidHF.h:88
Int_t fnSpecies
switch for Bayesian
Int_t fWhyRejection
PID for heavy flavours manager.
Definition: AliRDHFCuts.h:405
void SetUseDefaultPID(Bool_t defPID)
#define P(T, U, S)
UInt_t fUseImpParDCut
transverse momentum below which the strong PID is applied
double Double_t
Definition: External.C:58
Bool_t IsSignalMC(AliAODRecoDecay *d, AliAODEvent *aod, Int_t pdg) const
Bool_t IsTOFPiKexcluded(AliAODTrack *track, Double_t nsigmaK)
general method to perform PID using raw signals
Bool_t fUseKF
flag to switch on/off the default pid
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
Bool_t fUseMCVertex
flag to switch on the removal of duaghters from the primary vertex computation
Definition: AliRDHFCuts.h:408
void SetSigmaForTPC(Double_t *sigma)
Definition: AliAODPidHF.h:43
Double_t GetSigma(Int_t idet) const
Definition: AliAODPidHF.h:134
void SetMaxVtxZ(Float_t z=1e6)
Definition: AliRDHFCuts.h:61
Bool_t SetMCPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod) const
Bool_t fCombPID
max momentum for applying PID
Bool_t fRemoveDaughtersFromPrimary
Definition: AliRDHFCuts.h:407
void SetUseCentrality(Int_t flag=1)
virtual void PrintAll() const
virtual Int_t IsSelectedCombPID(AliAODRecoDecayHF *d)
Double_t NormalizedDecayLength2() const
Double_t ImpParXY() const
Double_t fmaxPtrackForPID
transverse momentum below which the special cuts are applied
void SetPCompatTOF(Double_t pTOF)
Definition: AliAODPidHF.h:106
virtual void SetStandardCutsPP2011_276TeV()
void SetNVars(Int_t nVars)
Definition: AliRDHFCuts.h:375
virtual void SetStandardCutsPP2010()
Double_t CosPointingAngleXY() const
Double_t fMaxRapidityCand
minimum pt of the candidate
Definition: AliRDHFCuts.h:426
void EnableSemiCentralTrigger()
Definition: AliRDHFCuts.h:94
AliRDHFCutsD0toKpi(const char *name="CutsD0toKpi")
void SetGlobalIndex()
Definition: AliRDHFCuts.h:197
AliRDHFCuts & operator=(const AliRDHFCuts &source)
Bool_t fLowPt
flag to switch on/off special cuts
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
Bool_t fUsePID
Definition: AliRDHFCuts.h:402
void SetSelectCandTrackSPDFirst(Bool_t flag, Double_t ptmax)
Flag and pt-maximum to check if the candidate daughters fulfill the kFirst criteria.
Definition: AliRDHFCuts.h:358
void SetPidHF(AliAODPidHF *pidObj)
see enum below
Definition: AliRDHFCuts.h:211
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:232
void SetTOF(Bool_t tof)
Definition: AliAODPidHF.h:95
Bool_t HasBadDaughters() const
Bool_t fDefaultPID
flag to switch on/off different pid for low pt D0
const Int_t nPtBins
void SetMinVtxContr(Int_t contr=1)
Definition: AliRDHFCuts.h:59
void Setd0MeasMinusExpCut(UInt_t nPtBins, Float_t *cutval)
int Int_t
Definition: External.C:63
Bool_t fUseTrackSelectionWithFilterBits
flag to reject kink daughters
Definition: AliRDHFCuts.h:436
void SetCuts(Int_t nVars, Int_t nPtBins, Float_t **cutsRD)
Int_t fnVarsForOpt
Definition: AliRDHFCuts.h:397
void SetUseKF(Bool_t useKF)
void ResetMaskAndEnableMBTrigger()
Definition: AliRDHFCuts.h:74
void SetMaximumPtSpecialCuts(Double_t pt)
void SetMinCentrality(Float_t minCentrality=0.)
Definition: AliRDHFCuts.h:51
unsigned int UInt_t
Definition: External.C:33
AliPIDCombined * GetPidCombined() const
Definition: AliAODPidHF.h:161
Double_t CosThetaStarD0bar() const
angle of K
float Float_t
Definition: External.C:68
const Double_t ptmax
Double_t fMaxPtCand
minimum pt of the candidate
Definition: AliRDHFCuts.h:425
Bool_t CheckStatus(AliAODTrack *track, TString detectors) const
Double_t sigmas[nPtBins]
AliAODVertex * GetOwnPrimaryVtx() const
Bool_t fKeepSignalMC
max rapidity of candidate (if !=-999 overrides IsInFiducialAcceptance)
Definition: AliRDHFCuts.h:427
Float_t * fMaxImpParD
switch for cut on D0 ImpParXY; =0 –> not used, >0 value represents array size (it has to coincide wit...
virtual void CalculateBayesianWeights(AliAODRecoDecayHF *d)
Int_t fIsSelectedPID
outcome of cuts selection
Definition: AliRDHFCuts.h:423
void SetImpParDCut(UInt_t nPtBins, Float_t *cutval)
Bool_t GetTOF() const
Definition: AliAODPidHF.h:142
Int_t MakeRawPid(AliAODTrack *track, Int_t specie)
Int_t mode
Definition: anaM.C:40
AliPIDResponse * GetPidResponse() const
Definition: AliAODPidHF.h:160
void SetSigma(Double_t *sigma)
Definition: AliAODPidHF.h:39
void SetSigmaForTOFCompat(Double_t sigma)
Definition: AliAODPidHF.h:45
Float_t * fCutsRD
fnVars*fnPtBins
Definition: AliRDHFCuts.h:400
Int_t IsSelectedKF(AliAODRecoDecayHF2Prong *d, AliAODEvent *aod) const
Double_t DecayLength2() const
kinematics & topology
void SetSigmaForTPCCompat(Double_t sigma)
Definition: AliAODPidHF.h:44
void SetMaxCentrality(Float_t maxCentrality=100.)
Definition: AliRDHFCuts.h:52
AliRDHFCutsD0toKpi & operator=(const AliRDHFCutsD0toKpi &source)
short Short_t
Definition: External.C:23
Int_t IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
Int_t CombineSelectionLevels(Int_t selectionvalTrack, Int_t selectionvalCand, Int_t selectionvalPID) const
void SetVarsForOpt(Int_t nVars, Bool_t *forOpt)
void SetVarNames(Int_t nVars, TString *varNames, Bool_t *isUpperCut)
Float_t * fMaxd0MeasMinusExp
switch for cut on d0meas-d0exp; =0 –> not used, >0 value represents array size (it has to coincide wi...
Bool_t * fVarsForOpt
number of cut vars to be optimized for candidates
Definition: AliRDHFCuts.h:398
Bool_t AreDaughtersSelected(AliAODRecoDecayHF *rd, const AliAODEvent *aod=0x0) const
void EnableCentralTrigger()
Definition: AliRDHFCuts.h:82
Bool_t IsPionRaw(AliAODTrack *track, TString detector) const
void SetUsePID(Bool_t flag=kTRUE)
Definition: AliRDHFCuts.h:205
Double_t Normalizedd0Prong(Int_t ip) const
virtual void PrintAll() const
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
virtual Int_t IsSelectedPID(AliAODRecoDecayHF *rd)
void SetFlatd0MeasMinusExpCut(Float_t value)
void CleanOwnPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod, AliAODVertex *origownvtx) const
void SetRemoveDaughtersFromPrim(Bool_t removeDaughtersPrim)
Definition: AliRDHFCuts.h:215
Bool_t GetCompat() const
Definition: AliAODPidHF.h:147
void SetUseSpecialCuts(Bool_t useSpecialCuts)
void SetOldPid(Bool_t oldPid)
Definition: AliAODPidHF.h:108
void SetPtBins(Int_t nPtBinLimits, Float_t *ptBinLimits)
virtual void SetStandardCutsPP2010vsMult()
void SetLowPt(Bool_t lowpt, Double_t ptlow=2.)
void SetFlatImpParDCut(Float_t value)
Double_t * fWeightsPositive
number of species (used only for array declaration)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Int_t fBayesianCondition
Switch for which Bayesian PID strategy to use.
void SetMatch(Int_t match)
Definition: AliAODPidHF.h:98
void AddTrackCuts(const AliESDtrackCuts *cuts)
Definition: AliRDHFCuts.h:203
void SetPLimit(Double_t *plim, Int_t npLim)
void SetTPC(Bool_t tpc)
Definition: AliAODPidHF.h:94
void SetMinPtCandidate(Double_t ptCand=-1.)
Definition: AliRDHFCuts.h:216
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
Int_t fnPtBins
cuts on the candidate
Definition: AliRDHFCuts.h:392
void SetCompat(Bool_t comp)
Definition: AliAODPidHF.h:100
void SetTriggerClass(TString trclass0, TString trclass1="")
Definition: AliRDHFCuts.h:193
AliAODPidHF * fPidHF
enable AOD049 centrality cleanup
Definition: AliRDHFCuts.h:404
Bool_t RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod) const
virtual void GetCutVarsForOpt(AliAODRecoDecayHF *d, Float_t *vars, Int_t nvars, Int_t *pdgdaughters)
Int_t PtBin(Double_t pt) const
void SetOptPileup(Int_t opt=0)
Definition: AliRDHFCuts.h:219
virtual void SetStandardCutsPbPb2010()
Int_t GetGlobalIndex(Int_t iVar, Int_t iPtBin) const
Double_t fPtLowPID
flag to switch on/off D0 selection via KF
Int_t nptbins
Double_t fMinPtCand
outcome of PID selection
Definition: AliRDHFCuts.h:424
Int_t IsSelectedPIDdefault(AliAODRecoDecayHF *rd)
virtual void SetStandardCutsPbPb2011()