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