AliPhysics  b76e98e (b76e98e)
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 
43 ClassImp(AliRDHFCutsD0toKpi);
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),
125  fMaxImpParD(0x0),
127  fMaxd0MeasMinusExp(0x0),
129  fLowPt(source.fLowPt),
130  fDefaultPID(source.fDefaultPID),
131  fUseKF(source.fUseKF),
132  fPtLowPID(source.fPtLowPID),
135  fCombPID(source.fCombPID),
136  fnSpecies(source.fnSpecies),
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:455
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:437
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:440
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:439
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:404
virtual void SetStandardCutsPP2010()
Double_t CosPointingAngleXY() const
Double_t fMaxRapidityCand
minimum pt of the candidate
Definition: AliRDHFCuts.h:459
void EnableSemiCentralTrigger()
Definition: AliRDHFCuts.h:97
AliRDHFCutsD0toKpi(const char *name="CutsD0toKpi")
void SetGlobalIndex()
Definition: AliRDHFCuts.h:201
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:434
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:386
void SetPidHF(AliAODPidHF *pidObj)
see enum below
Definition: AliRDHFCuts.h:215
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:247
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:469
void SetCuts(Int_t nVars, Int_t nPtBins, Float_t **cutsRD)
Int_t fnVarsForOpt
Definition: AliRDHFCuts.h:429
void SetUseKF(Bool_t useKF)
void ResetMaskAndEnableMBTrigger()
Definition: AliRDHFCuts.h:77
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:458
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:460
Float_t * fMaxImpParD
switch for cut on D0 ImpParXY; =0 –> not used, >0 value represents array size (it has to coincide wi...
virtual void CalculateBayesianWeights(AliAODRecoDecayHF *d)
Int_t fIsSelectedPID
outcome of cuts selection
Definition: AliRDHFCuts.h:456
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:41
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:432
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 w...
Bool_t * fVarsForOpt
number of cut vars to be optimized for candidates
Definition: AliRDHFCuts.h:430
Bool_t AreDaughtersSelected(AliAODRecoDecayHF *rd, const AliAODEvent *aod=0x0) const
void EnableCentralTrigger()
Definition: AliRDHFCuts.h:85
Bool_t IsPionRaw(AliAODTrack *track, TString detector) const
void SetUsePID(Bool_t flag=kTRUE)
Definition: AliRDHFCuts.h:209
Double_t Normalizedd0Prong(Int_t ip) const
virtual void PrintAll() const
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:219
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:207
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:220
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
Int_t fnPtBins
cuts on the candidate
Definition: AliRDHFCuts.h:424
void SetCompat(Bool_t comp)
Definition: AliAODPidHF.h:100
void SetTriggerClass(TString trclass0, TString trclass1="")
Definition: AliRDHFCuts.h:197
AliAODPidHF * fPidHF
enable AOD049 centrality cleanup
Definition: AliRDHFCuts.h:436
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:223
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:457
Int_t IsSelectedPIDdefault(AliAODRecoDecayHF *rd)
virtual void SetStandardCutsPbPb2011()