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