1 /**************************************************************************
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /////////////////////////////////////////////////////////////
20 // Class for cuts on AOD reconstructed D0->Kpi
22 // Author: A.Dainese, andrea.dainese@pd.infn.it
23 /////////////////////////////////////////////////////////////
25 #include <TDatabasePDG.h>
26 #include <Riostream.h>
28 #include "AliRDHFCutsD0toKpi.h"
29 #include "AliAODRecoDecayHF2Prong.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"
41 ClassImp(AliRDHFCutsD0toKpi)
43 //--------------------------------------------------------------------------
44 AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const char* name) :
46 fUseSpecialCuts(kFALSE),
51 fPtMaxSpecialCuts(9999.),
52 fmaxPtrackForPID(9999),
54 fnSpecies(AliPID::kSPECIES),
62 // Default Constructor
66 TString varNames[11]={"inv. mass [GeV]",
77 Bool_t isUpperCut[11]={kTRUE,
88 SetVarNames(nvars,varNames,isUpperCut);
89 Bool_t forOpt[11]={kFALSE,
100 SetVarsForOpt(4,forOpt);
101 Float_t limits[2]={0,999999999.};
104 fWeightsNegative = new Double_t[AliPID::kSPECIES];
105 fWeightsPositive = new Double_t[AliPID::kSPECIES];
107 for (Int_t i = 0; i<AliPID::kSPECIES; i++) {
108 fWeightsPositive[i] = 0.;
109 fWeightsNegative[i] = 0.;
113 //--------------------------------------------------------------------------
114 AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :
116 fUseSpecialCuts(source.fUseSpecialCuts),
117 fLowPt(source.fLowPt),
118 fDefaultPID(source.fDefaultPID),
119 fUseKF(source.fUseKF),
120 fPtLowPID(source.fPtLowPID),
121 fPtMaxSpecialCuts(source.fPtMaxSpecialCuts),
122 fmaxPtrackForPID(source.fmaxPtrackForPID),
123 fCombPID(source.fCombPID),
124 fnSpecies(source.fnSpecies),
125 fWeightsPositive(source.fWeightsPositive),
126 fWeightsNegative(source.fWeightsNegative),
127 fProbThreshold(source.fProbThreshold),
128 fBayesianStrategy(source.fBayesianStrategy),
129 fBayesianCondition(source.fBayesianCondition)
136 //--------------------------------------------------------------------------
137 AliRDHFCutsD0toKpi &AliRDHFCutsD0toKpi::operator=(const AliRDHFCutsD0toKpi &source)
140 // assignment operator
142 if(&source == this) return *this;
144 AliRDHFCuts::operator=(source);
145 fUseSpecialCuts=source.fUseSpecialCuts;
146 fLowPt=source.fLowPt;
147 fDefaultPID=source.fDefaultPID;
148 fUseKF=source.fUseKF;
149 fPtLowPID=source.fPtLowPID;
150 fPtMaxSpecialCuts=source.fPtMaxSpecialCuts;
151 fmaxPtrackForPID=source.fmaxPtrackForPID;
152 fCombPID = source.fCombPID;
153 fnSpecies = source.fnSpecies;
154 fWeightsPositive = source.fWeightsPositive;
155 fWeightsNegative = source.fWeightsNegative;
156 fProbThreshold = source.fProbThreshold;
157 fBayesianStrategy = source.fBayesianStrategy;
158 fBayesianCondition = source.fBayesianCondition;
164 //---------------------------------------------------------------------------
165 AliRDHFCutsD0toKpi::~AliRDHFCutsD0toKpi()
170 if (fWeightsNegative) {
171 delete[] fWeightsNegative;
172 fWeightsNegative = 0;
174 if (fWeightsPositive) {
175 delete[] fWeightsNegative;
176 fWeightsNegative = 0;
180 //---------------------------------------------------------------------------
181 void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent *aod) {
183 // Fills in vars the values of the variables
186 if(nvars!=fnVarsForOpt) {
187 printf("AliRDHFCutsD0toKpi::GetCutsVarsForOpt: wrong number of variables\n");
191 AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)d;
193 //recalculate vertex w/o daughters
194 Bool_t cleanvtx=kFALSE;
195 AliAODVertex *origownvtx=0x0;
196 if(fRemoveDaughtersFromPrimary) {
197 if(dd->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*dd->GetOwnPrimaryVtx());
199 if(!RecalcOwnPrimaryVtx(dd,aod)) {
200 CleanOwnPrimaryVtx(dd,aod,origownvtx);
208 if(TMath::Abs(pdgdaughters[0])==211) {
209 vars[iter]=dd->InvMassD0();
211 vars[iter]=dd->InvMassD0bar();
216 vars[iter]=dd->GetDCA();
220 if(TMath::Abs(pdgdaughters[0])==211) {
221 vars[iter] = dd->CosThetaStarD0();
223 vars[iter] = dd->CosThetaStarD0bar();
228 if(TMath::Abs(pdgdaughters[0])==321) {
229 vars[iter]=dd->PtProng(0);
232 vars[iter]=dd->PtProng(1);
237 if(TMath::Abs(pdgdaughters[0])==211) {
238 vars[iter]=dd->PtProng(0);
241 vars[iter]=dd->PtProng(1);
246 if(TMath::Abs(pdgdaughters[0])==321) {
247 vars[iter]=dd->Getd0Prong(0);
250 vars[iter]=dd->Getd0Prong(1);
255 if(TMath::Abs(pdgdaughters[0])==211) {
256 vars[iter]=dd->Getd0Prong(0);
259 vars[iter]=dd->Getd0Prong(1);
264 vars[iter]= dd->Prodd0d0();
268 vars[iter]=dd->CosPointingAngle();
273 vars[iter]=TMath::Abs(dd->CosPointingAngleXY());
278 vars[iter]=dd->NormalizedDecayLengthXY();
281 if(cleanvtx)CleanOwnPrimaryVtx(dd,aod,origownvtx);
285 //---------------------------------------------------------------------------
286 Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod) {
296 cout<<"Cut matrice not inizialized. Exit..."<<endl;
300 AliAODRecoDecayHF2Prong* d=(AliAODRecoDecayHF2Prong*)obj;
302 cout<<"AliAODRecoDecayHF2Prong null"<<endl;
306 if(fKeepSignalMC) if(IsSignalMC(d,aod,421)) return 3;
308 Double_t ptD=d->Pt();
309 if(ptD<fMinPtCand) return 0;
310 if(ptD>fMaxPtCand) return 0;
312 if(fUseTrackSelectionWithFilterBits && d->HasBadDaughters()) return 0;
314 // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
315 Int_t returnvaluePID=3;
316 Int_t returnvalueCuts=3;
318 // selection on candidate
319 if(selectionLevel==AliRDHFCuts::kAll ||
320 selectionLevel==AliRDHFCuts::kCandidate) {
324 //recalculate vertex w/o daughters
325 AliAODVertex *origownvtx=0x0;
326 if(fRemoveDaughtersFromPrimary && !fUseMCVertex) {
327 if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
328 if(!RecalcOwnPrimaryVtx(d,aod)) {
329 CleanOwnPrimaryVtx(d,aod,origownvtx);
335 if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
336 if(!SetMCPrimaryVtx(d,aod)) {
337 CleanOwnPrimaryVtx(d,aod,origownvtx);
344 Int_t okD0=0,okD0bar=0;
346 Int_t ptbin=PtBin(pt);
348 CleanOwnPrimaryVtx(d,aod,origownvtx);
352 Double_t mD0,mD0bar,ctsD0,ctsD0bar;
355 Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
357 d->InvMassD0(mD0,mD0bar);
358 if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
359 if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
360 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
362 if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
365 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;
366 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;
367 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
370 if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
371 TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
372 if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
373 TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
374 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
376 if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
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;}
384 if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
387 if(TMath::Abs(d->CosPointingAngleXY()) < fCutsRD[GetGlobalIndex(9,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
389 Double_t normalDecayLengXY=d->NormalizedDecayLengthXY();
390 if (normalDecayLengXY < fCutsRD[GetGlobalIndex(10, ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
392 if (returnvalueCuts!=0) {
393 if (okD0) returnvalueCuts=1; //cuts passed as D0
394 if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar
395 if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
400 if(fUseSpecialCuts && (pt<fPtMaxSpecialCuts)) special=IsSelectedSpecialCuts(d);
401 if(!special) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
403 // unset recalculated primary vertex when not needed any more
404 CleanOwnPrimaryVtx(d,aod,origownvtx);
407 // go to selection with Kalman vertexing, if requested
408 returnvalueCuts = IsSelectedKF(d,aod);
411 fIsSelectedCuts=returnvalueCuts;
412 if(!returnvalueCuts) return 0;
417 if(selectionLevel==AliRDHFCuts::kAll ||
418 selectionLevel==AliRDHFCuts::kCandidate ||
419 selectionLevel==AliRDHFCuts::kPID) {
421 returnvaluePID = IsSelectedPID(d);
422 fIsSelectedPID=returnvaluePID;
423 if(!returnvaluePID) return 0;
425 returnvaluePID = IsSelectedCombPID(d);
426 if(!returnvaluePID) return 0;
432 Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
434 if(!returnvalueComb) return 0;
436 // selection on daughter tracks
437 if(selectionLevel==AliRDHFCuts::kAll ||
438 selectionLevel==AliRDHFCuts::kTracks) {
439 if(!AreDaughtersSelected(d)) return 0;
442 // cout<<"Pid = "<<returnvaluePID<<endl;
443 return returnvalueComb;
446 //------------------------------------------------------------------------------------------
447 Int_t AliRDHFCutsD0toKpi::IsSelectedKF(AliAODRecoDecayHF2Prong *d,
448 AliAODEvent* aod) const {
450 // Apply selection using KF-vertexing
453 AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
454 AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
456 if(!track0 || !track1) {
457 cout<<"one or two D0 daughters missing!"<<endl;
461 // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
462 Int_t returnvalueCuts=3;
464 // candidate selection with AliKF
465 AliKFParticle::SetField(aod->GetMagneticField()); // set the magnetic field
467 Int_t okD0=0,okD0bar=0;
470 // convert tracks into AliKFParticles
472 AliKFParticle negPiKF(*track1,-211); // neg pion kandidate
473 AliKFParticle negKKF(*track1,-321); // neg kaon kandidate
474 AliKFParticle posPiKF(*track0,211); // pos pion kandidate
475 AliKFParticle posKKF(*track0,321); // pos kaon kandidate
477 // build D0 candidates
479 AliKFParticle d0c(negKKF,posPiKF); // D0 candidate
480 AliKFParticle ad0c(posKKF,negPiKF); // D0bar candidate
482 // create kf primary vertices
484 AliAODVertex *vtx1 = aod->GetPrimaryVertex();
485 AliKFVertex primVtx(*vtx1);
486 AliKFVertex aprimVtx(*vtx1);
488 if(primVtx.GetNContributors()<=0) okD0 = 0;
489 if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
490 if(!okD0 && !okD0bar) returnvalueCuts=0;
494 Double_t d0mass = d0c.GetMass();
495 Double_t ad0mass = ad0c.GetMass();
497 // calculate P of D0 and D0bar
498 Double_t d0P = d0c.GetP();
499 Double_t d0Px = d0c.GetPx();
500 Double_t d0Py = d0c.GetPy();
501 Double_t d0Pz = d0c.GetPz();
502 Double_t ad0P = ad0c.GetP();
503 Double_t ad0Px = ad0c.GetPx();
504 Double_t ad0Py = ad0c.GetPy();
505 Double_t ad0Pz = ad0c.GetPz();
507 //calculate Pt of D0 and D0bar
509 Double_t pt=d0c.GetPt();
510 Double_t apt=ad0c.GetPt();
512 // remove D0 daughters from primary vertices (if used in vertex fit) and add D0-candidates
514 if(track0->GetUsedForPrimVtxFit()) {
519 if(track1->GetUsedForPrimVtxFit()) {
527 if(primVtx.GetNContributors()<=0) okD0 = 0;
528 if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
529 if(!okD0 && !okD0bar) returnvalueCuts=0;
531 //calculate cut variables
533 // calculate impact params of daughters w.r.t recalculated vertices
535 Double_t impactPi = posPiKF.GetDistanceFromVertexXY(primVtx);
536 Double_t aimpactPi = negPiKF.GetDistanceFromVertexXY(aprimVtx);
537 Double_t impactKa = negKKF.GetDistanceFromVertexXY(primVtx);
538 Double_t aimpactKa = posKKF.GetDistanceFromVertexXY(aprimVtx);
540 // calculate Product of Impact Params
542 Double_t prodParam = impactPi*impactKa;
543 Double_t aprodParam = aimpactPi*aimpactKa;
545 // calculate cosine of pointing angles
547 TVector3 mom(d0c.GetPx(),d0c.GetPy(),d0c.GetPz());
548 TVector3 fline(d0c.GetX()-primVtx.GetX(),
549 d0c.GetY()-primVtx.GetY(),
550 d0c.GetZ()-primVtx.GetZ());
551 Double_t pta = mom.Angle(fline);
552 Double_t cosP = TMath::Cos(pta); // cosine of pta for D0 candidate
554 TVector3 amom(ad0c.GetPx(),ad0c.GetPy(),ad0c.GetPz());
555 TVector3 afline(ad0c.GetX()-aprimVtx.GetX(),
556 ad0c.GetY()-aprimVtx.GetY(),
557 ad0c.GetZ()-aprimVtx.GetZ());
558 Double_t apta = amom.Angle(afline);
559 Double_t acosP = TMath::Cos(apta); // cosine of pta for D0bar candidate
561 // calculate P of Pions at Decay Position of D0 and D0bar candidates
562 negKKF.TransportToParticle(d0c);
563 posPiKF.TransportToParticle(d0c);
564 posKKF.TransportToParticle(ad0c);
565 negPiKF.TransportToParticle(ad0c);
567 Double_t pxPi = posPiKF.GetPx();
568 Double_t pyPi = posPiKF.GetPy();
569 Double_t pzPi = posPiKF.GetPz();
570 Double_t ptPi = posPiKF.GetPt();
572 Double_t apxPi = negPiKF.GetPx();
573 Double_t apyPi = negPiKF.GetPy();
574 Double_t apzPi = negPiKF.GetPz();
575 Double_t aptPi = negPiKF.GetPt();
577 // calculate Pt of Kaons at Decay Position of D0 and D0bar candidates
579 Double_t ptK = negKKF.GetPt();
580 Double_t aptK = posKKF.GetPt();
582 //calculate cos(thetastar)
583 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(421)->Mass();
585 massp[0] = TDatabasePDG::Instance()->GetParticle(321)->Mass();
586 massp[1] = TDatabasePDG::Instance()->GetParticle(211)->Mass();
587 Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)
588 -4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
590 // cos(thetastar) for D0 and Pion
592 Double_t d0E = TMath::Sqrt(massvtx*massvtx + d0P*d0P);
593 Double_t beta = d0P/d0E;
594 Double_t gamma = d0E/massvtx;
595 TVector3 momPi(pxPi,pyPi,pzPi);
596 TVector3 momTot(d0Px,d0Py,d0Pz);
597 Double_t q1 = momPi.Dot(momTot)/momTot.Mag();
598 Double_t cts = (q1/gamma-beta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
600 // cos(thetastar) for D0bar and Pion
602 Double_t ad0E = TMath::Sqrt(massvtx*massvtx + ad0P*ad0P);
603 Double_t abeta = ad0P/ad0E;
604 Double_t agamma = ad0E/massvtx;
605 TVector3 amomPi(apxPi,apyPi,apzPi);
606 TVector3 amomTot(ad0Px,ad0Py,ad0Pz);
607 Double_t aq1 = amomPi.Dot(amomTot)/amomTot.Mag();
608 Double_t acts = (aq1/agamma-abeta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
610 // calculate reduced Chi2 for the full D0 fit
611 d0c.SetProductionVertex(primVtx);
612 ad0c.SetProductionVertex(aprimVtx);
613 negKKF.SetProductionVertex(d0c);
614 posPiKF.SetProductionVertex(d0c);
615 posKKF.SetProductionVertex(ad0c);
616 negPiKF.SetProductionVertex(ad0c);
617 d0c.TransportToProductionVertex();
618 ad0c.TransportToProductionVertex();
620 // calculate the decay length
621 Double_t decayLengthD0 = d0c.GetDecayLength();
622 Double_t adecayLengthD0 = ad0c.GetDecayLength();
624 Double_t chi2D0 = 50.;
625 if(d0c.GetNDF() > 0 && d0c.GetChi2() >= 0) {
626 chi2D0 = d0c.GetChi2()/d0c.GetNDF();
629 Double_t achi2D0 = 50.;
630 if(ad0c.GetNDF() > 0 && ad0c.GetChi2() >= 0) {
631 achi2D0 = ad0c.GetChi2()/ad0c.GetNDF();
635 Int_t ptbin=PtBin(pt);
636 Int_t aptbin=PtBin(apt);
638 if(ptbin < 0) okD0 = 0;
639 if(aptbin < 0) okD0bar = 0;
640 if(!okD0 && !okD0bar) returnvalueCuts=0;
642 if(ptK < fCutsRD[GetGlobalIndex(3,ptbin)] || ptPi < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
643 if(aptK < fCutsRD[GetGlobalIndex(3,aptbin)] || aptPi < fCutsRD[GetGlobalIndex(4,aptbin)]) okD0bar = 0;
644 if(!okD0 && !okD0bar) returnvalueCuts=0;
647 if(TMath::Abs(impactKa) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
648 TMath::Abs(impactPi) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
650 if(TMath::Abs(aimpactKa) > fCutsRD[GetGlobalIndex(5,aptbin)] ||
651 TMath::Abs(aimpactPi) > fCutsRD[GetGlobalIndex(6,aptbin)]) okD0bar = 0;
653 if(!okD0 && !okD0bar) returnvalueCuts=0;
655 // for the moment via the standard method due to bug in AliKF
656 if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) okD0 = 0;
657 if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,aptbin)]) okD0bar = 0;
658 if(!okD0 && !okD0bar) returnvalueCuts=0;
661 if(TMath::Abs(d0mass-massvtx) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
662 if(TMath::Abs(ad0mass-massvtx) > fCutsRD[GetGlobalIndex(0,aptbin)]) okD0bar = 0;
663 if(!okD0 && !okD0bar) returnvalueCuts=0;
666 if(TMath::Abs(cts) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
667 if(TMath::Abs(acts) > fCutsRD[GetGlobalIndex(2,aptbin)]) okD0bar = 0;
668 if(!okD0 && !okD0bar) returnvalueCuts=0;
670 if(prodParam > fCutsRD[GetGlobalIndex(7,ptbin)]) okD0 = 0;
671 if(aprodParam > fCutsRD[GetGlobalIndex(7,aptbin)]) okD0bar = 0;
672 if(!okD0 && !okD0bar) returnvalueCuts=0;
674 if(cosP < fCutsRD[GetGlobalIndex(8,ptbin)]) okD0 = 0;
675 if(acosP < fCutsRD[GetGlobalIndex(8,aptbin)]) okD0bar = 0;
676 if(!okD0 && !okD0bar) returnvalueCuts=0;
678 if(chi2D0 > fCutsRD[GetGlobalIndex(10,ptbin)]) okD0 = 0;
679 if(achi2D0 > fCutsRD[GetGlobalIndex(10,aptbin)]) okD0bar = 0;
680 if(!okD0 && !okD0bar) returnvalueCuts=0;
682 if(decayLengthD0 < fCutsRD[GetGlobalIndex(9,ptbin)]) okD0 = 0;
683 if(adecayLengthD0 < fCutsRD[GetGlobalIndex(9,aptbin)]) okD0bar = 0;
684 if(!okD0 && !okD0bar) returnvalueCuts=0;
686 if(returnvalueCuts!=0) {
687 if(okD0) returnvalueCuts=1; //cuts passed as D0
688 if(okD0bar) returnvalueCuts=2; //cuts passed as D0bar
689 if(okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
692 return returnvalueCuts;
695 //---------------------------------------------------------------------------
697 Bool_t AliRDHFCutsD0toKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
700 // Checking if D0 is in fiducial acceptance region
703 if(fMaxRapidityCand>-998.){
704 if(TMath::Abs(y) > fMaxRapidityCand) return kFALSE;
709 // applying cut for pt > 5 GeV
710 AliDebug(2,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt));
711 if (TMath::Abs(y) > 0.8){
715 // appliying smooth cut for pt < 5 GeV
716 Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
717 Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
718 AliDebug(2,Form("pt of D0 = %f (< 5), cutting according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
719 if (y < minFiducialY || y > maxFiducialY){
726 //---------------------------------------------------------------------------
727 Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d)
729 // ############################################################
731 // Apply PID selection
734 // ############################################################
736 if(!fUsePID) return 3;
737 if(fDefaultPID) return IsSelectedPIDdefault(d);
738 if (fCombPID) return IsSelectedCombPID(d); //to use Bayesian method if applicable
740 Int_t isD0D0barPID[2]={1,2};
741 Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
743 // values convention -1 = discarded
744 // 0 = not identified (but compatible) || No PID (->hasPID flag)
746 // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
747 // Initial hypothesis: unknwon (but compatible)
748 combinedPID[0][0]=0; // prima figlia, Kaon
749 combinedPID[0][1]=0; // prima figlia, pione
750 combinedPID[1][0]=0; // seconda figlia, Kaon
751 combinedPID[1][1]=0; // seconda figlia, pion
753 Bool_t checkPIDInfo[2]={kTRUE,kTRUE};
754 Double_t sigma_tmp[3]={fPidHF->GetSigma(0),fPidHF->GetSigma(1),fPidHF->GetSigma(2)};
755 Bool_t isTOFused=fPidHF->GetTOF(),isCompat=fPidHF->GetCompat();
756 for(Int_t daught=0;daught<2;daught++){
758 AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
759 if(fPidHF->IsTOFPiKexcluded(aodtrack,5.)) return 0;
761 if(!(fPidHF->CheckStatus(aodtrack,"TPC")) && !(fPidHF->CheckStatus(aodtrack,"TOF"))) {
762 checkPIDInfo[daught]=kFALSE;
765 Double_t pProng=aodtrack->P();
768 if(pProng<fmaxPtrackForPID){
769 combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
772 if(pProng<fmaxPtrackForPID){
773 if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
774 combinedPID[daught][1]=0;
776 fPidHF->SetTOF(kFALSE);
777 combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
778 if(isTOFused)fPidHF->SetTOF(kTRUE);
779 if(isCompat)fPidHF->SetCompat(kTRUE);
783 if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded
787 else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){
788 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
789 else isD0D0barPID[0]=0;// if K+ D0 excluded
791 /* else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){
796 else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
797 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
798 else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
800 else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
801 if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
802 else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
805 if(fLowPt && d->Pt()<fPtLowPID){
806 Double_t sigmaTPC[3]={3.,2.,0.};
807 fPidHF->SetSigmaForTPC(sigmaTPC);
809 combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
812 fPidHF->SetCompat(kFALSE);
813 combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
814 fPidHF->SetCompat(kTRUE);
817 if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
818 combinedPID[daught][1]=0;
820 fPidHF->SetTOF(kFALSE);
821 Double_t sigmaTPCpi[3]={3.,3.,0.};
822 fPidHF->SetSigmaForTPC(sigmaTPCpi);
823 combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
824 fPidHF->SetTOF(kTRUE);
826 Bool_t isTPCpion=fPidHF->IsPionRaw(aodtrack,"TPC");
828 combinedPID[daught][1]=1;
830 combinedPID[daught][1]=-1;
836 fPidHF->SetSigmaForTPC(sigma_tmp);
837 }// END OF LOOP ON DAUGHTERS
839 if(!checkPIDInfo[0] && !checkPIDInfo[1]) {
840 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
845 // FURTHER PID REQUEST (both daughter info is needed)
846 if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
847 fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
848 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
852 if(fLowPt && d->Pt()<fPtLowPID){
853 if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
854 fWhyRejection=32;// reject cases where the Kaon is not identified
855 fPidHF->SetSigmaForTPC(sigma_tmp);
859 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
861 // cout<<"Why? "<<fWhyRejection<<endl;
862 return isD0D0barPID[0]+isD0D0barPID[1];
864 //---------------------------------------------------------------------------
865 Int_t AliRDHFCutsD0toKpi::IsSelectedPIDdefault(AliAODRecoDecayHF* d)
867 // ############################################################
869 // Apply PID selection
872 // temporary selection: PID AS USED FOR D0 by Andrea Rossi (up to 28/06/2010)
874 // d must be a AliAODRecoDecayHF2Prong object
875 // returns 0 if both D0 and D0bar are rejectecd
876 // 1 if D0 is accepted while D0bar is rejected
877 // 2 if D0bar is accepted while D0 is rejected
878 // 3 if both are accepted
879 // fWhyRejection variable (not returned for the moment, print it if needed)
880 // keeps some information on why a candidate has been
881 // rejected according to the following (unfriendly?) scheme
882 // if more rejection cases are considered interesting, just add numbers
884 // TO BE CONSIDERED WITH A GRAIN OF SALT (the order in which cut are applied is relevant)
885 // from 20 to 30: "detector" selection (PID acceptance)
888 // 28: no (TOF||TPC) pid information (no kTOFpid,kTOFout,kTIME,kTPCpid,...)
890 // from 30 to 40: PID selection
891 // 31: no Kaon compatible tracks found between daughters
892 // 32: no Kaon identified tracks found (strong sel. at low momenta)
893 // 33: both mass hypotheses are rejected
895 // ############################################################
897 if(!fUsePID) return 3;
899 Int_t isD0D0barPID[2]={1,2};
900 Double_t nsigmaTPCpi=-1., nsigmaTPCK=-1.; //used for TPC pid
901 Double_t tofSig,times[5];// used fot TOF pid
902 Int_t hasPID[2]={2,2};// flag to count how many detectors give PID info for the daughters
903 Int_t isKaonPionTOF[2][2],isKaonPionTPC[2][2];
904 Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
906 // values convention -1 = discarded
907 // 0 = not identified (but compatible) || No PID (->hasPID flag)
909 // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
910 // Initial hypothesis: unknwon (but compatible)
911 isKaonPionTOF[0][0]=0;
912 isKaonPionTOF[0][1]=0;
913 isKaonPionTOF[1][0]=0;
914 isKaonPionTOF[1][1]=0;
916 isKaonPionTPC[0][0]=0;
917 isKaonPionTPC[0][1]=0;
918 isKaonPionTPC[1][0]=0;
919 isKaonPionTPC[1][1]=0;
928 for(Int_t daught=0;daught<2;daught++){
931 // ########### Step 0- CHECKING minimal PID "ACCEPTANCE" ####################
933 AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
935 if(!(aodtrack->GetStatus()&AliESDtrack::kTPCrefit)){
939 if(!(aodtrack->GetStatus()&AliESDtrack::kITSrefit)){
944 AliAODPid *pid=aodtrack->GetDetPid();
951 // ########### Step 1- Check of TPC and TOF response ####################
953 Double_t ptrack=aodtrack->P();
954 //#################### TPC PID #######################
955 if (!(aodtrack->GetStatus()&AliESDtrack::kTPCpid )){
956 // NO TPC PID INFO FOR THIS TRACK
960 static AliTPCPIDResponse theTPCpid;
961 AliAODPid *pidObj = aodtrack->GetDetPid();
962 Double_t ptProng=pidObj->GetTPCmomentum();
963 nsigmaTPCpi = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kPion);
964 nsigmaTPCK = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kKaon);
967 if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
968 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
969 if(TMath::Abs(nsigmaTPCpi)<2.)isKaonPionTPC[daught][1]=1;
970 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
972 //else if(ptrack<.8){
974 if(TMath::Abs(nsigmaTPCK)<1.)isKaonPionTPC[daught][0]=1;
975 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
976 if(TMath::Abs(nsigmaTPCpi)<1.)isKaonPionTPC[daught][1]=1;
977 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
980 // if(nsigmaTPCK>-2.&&nsigmaTPCK<1.)isKaonPionTPC[daught][0]=1;
981 if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
982 //if(nsigmaTPCpi>-1.&&nsigmaTPCpi<2.)isKaonPionTPC[daught][1]=1;
983 if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
988 // ##### TOF PID: do not ask nothing for pion/protons ############
989 if(!((aodtrack->GetStatus()&AliESDtrack::kTOFpid)&&(aodtrack->GetStatus()&AliESDtrack::kTOFout)&&(aodtrack->GetStatus()&AliESDtrack::kTIME))){
990 // NO TOF PID INFO FOR THIS TRACK
994 tofSig=pid->GetTOFsignal();
995 pid->GetIntegratedTimes(times);
996 if((tofSig-times[3])>5.*160.)return 0;// PROTON REJECTION
997 if(TMath::Abs(tofSig-times[3])>3.*160.){
998 isKaonPionTOF[daught][0]=-1;
1002 isKaonPionTOF[daught][0]=1;
1007 //######### Step 2: COMBINE TOF and TPC PID ###############
1008 // we apply the following convention: if TPC and TOF disagree (discarded Vs identified) -> unknown
1009 combinedPID[daught][0]=isKaonPionTOF[daught][0]+isKaonPionTPC[daught][0];
1010 combinedPID[daught][1]=isKaonPionTOF[daught][1]+isKaonPionTPC[daught][1];
1013 //######### Step 3: USE PID INFO
1015 if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){// if not a K- and not a pi- both D0 and D0bar excluded
1019 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
1020 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
1021 else isD0D0barPID[0]=0;// if K+ D0 excluded
1023 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
1027 else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
1028 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
1029 else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
1031 else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
1032 if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
1033 else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
1036 // ########## ALSO DIFFERENT TPC PID REQUEST FOR LOW pt D0: request of K identification ###############################
1037 // ########## more tolerant criteria for single particle ID-> more selective criteria for D0 ##############################
1038 // ############### NOT OPTIMIZED YET ###################################
1040 isKaonPionTPC[daught][0]=0;
1041 isKaonPionTPC[daught][1]=0;
1042 AliAODPid *pidObj = aodtrack->GetDetPid();
1043 Double_t ptProng=pidObj->GetTPCmomentum();
1046 if(TMath::Abs(nsigmaTPCK)<3.)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;
1051 //else if(ptrack<.8){
1052 else if(ptProng<.8){
1053 if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
1054 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1055 if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
1056 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1059 if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1060 if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1064 }// END OF LOOP ON DAUGHTERS
1066 // FURTHER PID REQUEST (both daughter info is needed)
1067 if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
1068 fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
1071 else if(hasPID[0]==0&&hasPID[1]==0){
1072 fWhyRejection=28;// reject cases in which no PID info is available
1076 // request of K identification at low D0 pt
1077 combinedPID[0][0]=0;
1078 combinedPID[0][1]=0;
1079 combinedPID[1][0]=0;
1080 combinedPID[1][1]=0;
1082 combinedPID[0][0]=isKaonPionTOF[0][0]+isKaonPionTPC[0][0];
1083 combinedPID[0][1]=isKaonPionTOF[0][1]+isKaonPionTPC[0][1];
1084 combinedPID[1][0]=isKaonPionTOF[1][0]+isKaonPionTPC[1][0];
1085 combinedPID[1][1]=isKaonPionTOF[1][1]+isKaonPionTPC[1][1];
1087 if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
1088 fWhyRejection=32;// reject cases where the Kaon is not identified
1093 // cout<<"Why? "<<fWhyRejection<<endl;
1094 return isD0D0barPID[0]+isD0D0barPID[1];
1099 //---------------------------------------------------------------------------
1100 Int_t AliRDHFCutsD0toKpi::CombineSelectionLevels(Int_t selectionvalTrack,
1101 Int_t selectionvalCand,
1102 Int_t selectionvalPID) const
1105 // This method combines the tracks, PID and cuts selection results
1107 if(selectionvalTrack==0) return 0;
1111 switch(selectionvalPID) {
1116 returnvalue=((selectionvalCand==1 || selectionvalCand==3) ? 1 : 0);
1119 returnvalue=((selectionvalCand==2 || selectionvalCand==3) ? 2 : 0);
1122 returnvalue=selectionvalCand;
1131 //----------------------------------------------------------------------------
1132 Int_t AliRDHFCutsD0toKpi::IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
1135 // Note: this method is temporary
1136 // Additional cuts on decay lenght and lower cut for d0 norm are applied using vertex without candidate's daughters
1141 Float_t normDecLengthCut=1.,decLengthCut=TMath::Min(d->P()*0.0066+0.01,0.06/*cm*/), normd0Cut=0.5;
1142 // "decay length" expo law with tau' = beta*gamma*ctau= p/m*ctau =p*0.0123/1.864~p*0.0066
1143 // decay lenght > ctau' implies to retain (1-1/e) (for signal without considering detector resolution),
1145 Int_t returnvalue=3; //cut passed
1146 for(Int_t i=0;i<2/*prongs*/;i++){
1147 if(TMath::Abs(d->Normalizedd0Prong(i))<normd0Cut) return 0; //normd0Cut not passed
1149 if(d->DecayLength2()<decLengthCut*decLengthCut) return 0; //decLengthCut not passed
1150 if(d->NormalizedDecayLength2()<normDecLengthCut*normDecLengthCut) return 0; //decLengthCut not passed
1155 //----------------------------------------------
1156 void AliRDHFCutsD0toKpi::SetUseKF(Bool_t useKF)
1158 //switch on candidate selection via AliKFparticle
1164 TString varNamesKF[11]={"inv. mass [GeV]",
1175 Bool_t isUpperCutKF[11]={kTRUE,
1186 SetVarNames(nvarsKF,varNamesKF,isUpperCutKF);
1188 Bool_t forOpt[11]={kFALSE,
1199 SetVarsForOpt(4,forOpt);
1205 //---------------------------------------------------------------------------
1206 void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {
1208 //STANDARD CUTS USED FOR 2010 pp analysis
1209 //dca cut will be enlarged soon to 400 micron
1212 SetName("D0toKpiCutsStandard");
1213 SetTitle("Standard Cuts for D0 analysis");
1215 // PILE UP REJECTION
1216 SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1223 // TRACKS ON SINGLE TRACKS
1224 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1225 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1226 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1227 esdTrackCuts->SetRequireITSRefit(kTRUE);
1228 // esdTrackCuts->SetMinNClustersITS(4);
1229 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1230 esdTrackCuts->SetMinDCAToVertexXY(0.);
1231 esdTrackCuts->SetEtaRange(-0.8,0.8);
1232 esdTrackCuts->SetPtRange(0.3,1.e10);
1234 AddTrackCuts(esdTrackCuts);
1235 delete esdTrackCuts;
1238 const Int_t nptbins =14;
1239 const Double_t ptmax = 9999.;
1240 const Int_t nvars=11;
1241 Float_t ptbins[nptbins+1];
1258 SetGlobalIndex(nvars,nptbins);
1259 SetPtBins(nptbins+1,ptbins);
1261 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*/
1262 {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*/
1263 {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 */
1264 {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 */
1265 {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 */
1266 {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 */
1267 {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 */
1268 {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 */
1269 {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 */
1270 {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 */
1271 {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 */
1272 {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 */
1273 {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 */
1274 {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 */
1277 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1278 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1279 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1281 for (Int_t ibin=0;ibin<nptbins;ibin++){
1282 for (Int_t ivar = 0; ivar<nvars; ivar++){
1283 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1287 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1288 SetUseSpecialCuts(kTRUE);
1289 SetRemoveDaughtersFromPrim(kTRUE);
1291 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1292 delete [] cutsMatrixTransposeStand;
1293 cutsMatrixTransposeStand=NULL;
1296 AliAODPidHF* pidObj=new AliAODPidHF();
1297 //pidObj->SetName("pid4D0");
1299 const Int_t nlims=2;
1300 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1301 Bool_t compat=kTRUE; //effective only for this mode
1303 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1304 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1305 pidObj->SetMatch(mode);
1306 pidObj->SetPLimit(plims,nlims);
1307 pidObj->SetSigma(sigmas);
1308 pidObj->SetCompat(compat);
1309 pidObj->SetTPC(kTRUE);
1310 pidObj->SetTOF(kTRUE);
1311 pidObj->SetPCompatTOF(1.5);
1312 pidObj->SetSigmaForTPCCompat(3.);
1313 pidObj->SetSigmaForTOFCompat(3.);
1314 pidObj->SetOldPid(kTRUE);
1318 SetUseDefaultPID(kFALSE);
1331 //___________________________________________________________________________
1332 void AliRDHFCutsD0toKpi::SetStandardCutsPP2010vsMult() {
1334 // Cuts for 2010 pp 7 TeV data analysis in Ntracklets bins (vs multiplicity)
1336 SetName("D0toKpiCuts");
1337 SetTitle("Cuts for D0 analysis in 2010-data pp 7 TeV vs multiplicity");
1342 AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
1343 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1345 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1346 esdTrackCuts->SetRequireITSRefit(kTRUE);
1347 esdTrackCuts->SetEtaRange(-0.8,0.8);
1348 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1349 AliESDtrackCuts::kAny);
1350 // default is kBoth, otherwise kAny
1351 esdTrackCuts->SetMinDCAToVertexXY(0.);
1352 esdTrackCuts->SetPtRange(0.3,1.e10);
1354 AddTrackCuts(esdTrackCuts);
1355 delete esdTrackCuts;
1360 // Cut values per pt bin
1362 const Int_t nvars=11;
1363 const Int_t nptbins=14;
1365 ptbins=new Float_t[nptbins+1];
1382 SetPtBins(nptbins+1,ptbins);
1384 //setting cut values
1385 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*/
1386 {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*/
1387 {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 */
1388 {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 */
1389 {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 */
1390 {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 */
1391 {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 */
1392 {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 */
1393 {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 */
1394 {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 */
1395 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 12<pt<16 */
1396 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 16<pt<20 */
1397 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 20<pt<24 */
1398 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.}};/* pt>24 */
1401 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1402 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1403 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1405 for (Int_t ibin=0;ibin<nptbins;ibin++){
1406 for (Int_t ivar = 0; ivar<nvars; ivar++){
1407 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1411 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1412 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1413 delete [] cutsMatrixTransposeStand;
1415 SetUseSpecialCuts(kTRUE);
1416 SetMaximumPtSpecialCuts(8.);
1418 //Do recalculate the vertex
1419 SetRemoveDaughtersFromPrim(kTRUE);
1424 Bool_t pidflag=kTRUE;
1426 if(pidflag) cout<<"PID is used"<<endl;
1427 else cout<<"PID is not used"<<endl;
1429 AliAODPidHF* pidObj=new AliAODPidHF();
1431 const Int_t nlims=2;
1432 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1433 Bool_t compat=kTRUE; //effective only for this mode
1435 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1436 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1437 pidObj->SetMatch(mode);
1438 pidObj->SetPLimit(plims,nlims);
1439 pidObj->SetSigma(sigmas);
1440 pidObj->SetCompat(compat);
1441 pidObj->SetTPC(kTRUE);
1442 pidObj->SetTOF(kTRUE);
1443 pidObj->SetPCompatTOF(1.5);
1444 pidObj->SetSigmaForTPCCompat(3.);
1445 pidObj->SetSigmaForTOFCompat(3.);
1446 pidObj->SetOldPid(kTRUE);
1447 // pidObj->SetOldPid(kFALSE);
1450 SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
1454 //activate pileup rejection (for pp)
1455 SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1468 //---------------------------------------------------------------------------
1469 void AliRDHFCutsD0toKpi::SetStandardCutsPP2011_276TeV() {
1471 // STANDARD CUTS USED FOR 2011 pp analysis at 2.76TeV
1474 SetName("D0toKpiCutsStandard");
1475 SetTitle("Standard Cuts for D0 analysis in pp2011 at 2.76TeV run");
1480 AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
1481 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1483 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1484 esdTrackCuts->SetRequireITSRefit(kTRUE);
1485 esdTrackCuts->SetEtaRange(-0.8,0.8);
1486 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1487 AliESDtrackCuts::kAny);
1488 // default is kBoth, otherwise kAny
1489 esdTrackCuts->SetMinDCAToVertexXY(0.);
1490 esdTrackCuts->SetPtRange(0.3,1.e10);
1492 esdTrackCuts->SetMaxDCAToVertexXY(1.);
1493 esdTrackCuts->SetMaxDCAToVertexZ(1.);
1494 esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1496 AddTrackCuts(esdTrackCuts);
1497 delete esdTrackCuts;
1501 const Int_t nvars=11;
1502 const Int_t nptbins=13;
1503 Float_t ptbins[nptbins+1];
1519 SetPtBins(nptbins+1,ptbins);
1522 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*/
1523 {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*/
1524 {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 */
1525 {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 */
1526 {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 */
1527 {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 */
1528 {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 */
1529 {0.400,0.09,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* 6<pt<8 */
1530 {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 */
1531 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 12<pt<16 */
1532 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 16<pt<20 */
1533 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 20<pt<24 */
1534 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.}};/* pt>24 */
1536 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1537 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1538 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1539 for (Int_t ibin=0;ibin<nptbins;ibin++){
1540 for (Int_t ivar = 0; ivar<nvars; ivar++){
1541 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1544 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1545 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1546 delete [] cutsMatrixTransposeStand;
1550 AliAODPidHF* pidObj=new AliAODPidHF();
1552 const Int_t nlims=2;
1553 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1554 Bool_t compat=kTRUE; //effective only for this mode
1556 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1557 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1558 pidObj->SetMatch(mode);
1559 pidObj->SetPLimit(plims,nlims);
1560 pidObj->SetSigma(sigmas);
1561 pidObj->SetCompat(compat);
1562 pidObj->SetTPC(kTRUE);
1563 pidObj->SetTOF(kTRUE);
1564 pidObj->SetOldPid(kTRUE);
1567 SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
1571 //activate pileup rejection (for pp)
1572 SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1574 //Do recalculate the vertex
1575 SetRemoveDaughtersFromPrim(kTRUE);
1580 // Use the kFirst selection for tracks with candidate D with pt<2
1581 SetSelectCandTrackSPDFirst(kTRUE,2.);
1583 // Use special cuts for D candidates with pt<2
1584 SetUseSpecialCuts(kTRUE);
1585 SetMaximumPtSpecialCuts(2.);
1596 //---------------------------------------------------------------------------
1597 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
1599 // Final CUTS USED FOR 2010 PbPb analysis of central events
1600 // REMEMBER TO EVENTUALLY SWITCH ON
1601 // SetUseAOD049(kTRUE);
1604 SetName("D0toKpiCutsStandard");
1605 SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
1608 // CENTRALITY SELECTION
1609 SetMinCentrality(0.);
1610 SetMaxCentrality(80.);
1611 SetUseCentrality(AliRDHFCuts::kCentV0M);
1620 SetOptPileup(AliRDHFCuts::kNoPileupSelection);
1621 // PILE UP REJECTION
1622 //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1624 // TRACKS ON SINGLE TRACKS
1625 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1626 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1627 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1628 esdTrackCuts->SetRequireITSRefit(kTRUE);
1629 // esdTrackCuts->SetMinNClustersITS(4);
1630 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1631 esdTrackCuts->SetMinDCAToVertexXY(0.);
1632 esdTrackCuts->SetEtaRange(-0.8,0.8);
1633 esdTrackCuts->SetPtRange(0.7,1.e10);
1635 esdTrackCuts->SetMaxDCAToVertexXY(1.);
1636 esdTrackCuts->SetMaxDCAToVertexZ(1.);
1637 esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1640 AddTrackCuts(esdTrackCuts);
1641 delete esdTrackCuts;
1643 // SPD k FIRST for Pt<3 GeV/c
1644 SetSelectCandTrackSPDFirst(kTRUE, 3);
1647 const Int_t nptbins =13;
1648 const Double_t ptmax = 9999.;
1649 const Int_t nvars=11;
1650 Float_t ptbins[nptbins+1];
1666 SetGlobalIndex(nvars,nptbins);
1667 SetPtBins(nptbins+1,ptbins);
1668 SetMinPtCandidate(2.);
1670 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*/
1671 {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*/
1672 {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 */
1673 {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 */
1674 {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 */
1675 {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 */
1676 {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 */
1677 {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 */
1678 {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 */
1679 {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 */
1680 {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 */
1681 {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 */
1682 {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 */
1685 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1686 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1687 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1689 for (Int_t ibin=0;ibin<nptbins;ibin++){
1690 for (Int_t ivar = 0; ivar<nvars; ivar++){
1691 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1695 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1696 SetUseSpecialCuts(kTRUE);
1697 SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1698 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1699 delete [] cutsMatrixTransposeStand;
1700 cutsMatrixTransposeStand=NULL;
1703 AliAODPidHF* pidObj=new AliAODPidHF();
1704 //pidObj->SetName("pid4D0");
1706 const Int_t nlims=2;
1707 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1708 Bool_t compat=kTRUE; //effective only for this mode
1710 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1711 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1712 pidObj->SetMatch(mode);
1713 pidObj->SetPLimit(plims,nlims);
1714 pidObj->SetSigma(sigmas);
1715 pidObj->SetCompat(compat);
1716 pidObj->SetTPC(kTRUE);
1717 pidObj->SetTOF(kTRUE);
1718 pidObj->SetPCompatTOF(2.);
1719 pidObj->SetSigmaForTPCCompat(3.);
1720 pidObj->SetSigmaForTOFCompat(3.);
1721 pidObj->SetOldPid(kTRUE);
1726 SetUseDefaultPID(kFALSE);
1739 //---------------------------------------------------------------------------
1740 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() {
1741 // CUTS USED FOR D0 RAA for the CENTRALITY RANGE 20-80%
1744 SetName("D0toKpiCutsStandard");
1745 SetTitle("Standard Cuts for D0 analysis in PbPb2010 run, for peripheral events");
1748 // CENTRALITY SELECTION
1749 SetMinCentrality(40.);
1750 SetMaxCentrality(80.);
1751 SetUseCentrality(AliRDHFCuts::kCentV0M);
1758 SetOptPileup(AliRDHFCuts::kNoPileupSelection);
1759 // PILE UP REJECTION
1760 //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1762 // TRACKS ON SINGLE TRACKS
1763 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1764 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1765 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1766 esdTrackCuts->SetRequireITSRefit(kTRUE);
1767 // esdTrackCuts->SetMinNClustersITS(4);
1768 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1769 esdTrackCuts->SetMinDCAToVertexXY(0.);
1770 esdTrackCuts->SetEtaRange(-0.8,0.8);
1771 esdTrackCuts->SetPtRange(0.5,1.e10);
1773 esdTrackCuts->SetMaxDCAToVertexXY(1.);
1774 esdTrackCuts->SetMaxDCAToVertexZ(1.);
1775 esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0025*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1778 AddTrackCuts(esdTrackCuts);
1779 delete esdTrackCuts;
1781 // SPD k FIRST for Pt<3 GeV/c
1782 SetSelectCandTrackSPDFirst(kTRUE, 3);
1785 const Int_t nptbins =13;
1786 const Double_t ptmax = 9999.;
1787 const Int_t nvars=11;
1788 Float_t ptbins[nptbins+1];
1804 SetGlobalIndex(nvars,nptbins);
1805 SetPtBins(nptbins+1,ptbins);
1806 SetMinPtCandidate(0.);
1808 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*/
1809 {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*/
1810 {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 */
1811 {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 */
1812 {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 */
1813 {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 */
1814 {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 */
1815 {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 */
1816 {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 */
1817 {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 */
1818 {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 */
1819 {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 */
1820 {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 */
1823 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1824 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1825 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1827 for (Int_t ibin=0;ibin<nptbins;ibin++){
1828 for (Int_t ivar = 0; ivar<nvars; ivar++){
1829 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1833 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1834 SetUseSpecialCuts(kTRUE);
1835 SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1836 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1837 delete [] cutsMatrixTransposeStand;
1838 cutsMatrixTransposeStand=NULL;
1841 AliAODPidHF* pidObj=new AliAODPidHF();
1842 //pidObj->SetName("pid4D0");
1844 const Int_t nlims=2;
1845 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1846 Bool_t compat=kTRUE; //effective only for this mode
1848 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1849 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1850 pidObj->SetMatch(mode);
1851 pidObj->SetPLimit(plims,nlims);
1852 pidObj->SetSigma(sigmas);
1853 pidObj->SetCompat(compat);
1854 pidObj->SetTPC(kTRUE);
1855 pidObj->SetTOF(kTRUE);
1856 pidObj->SetPCompatTOF(2.);
1857 pidObj->SetSigmaForTPCCompat(3.);
1858 pidObj->SetSigmaForTOFCompat(3.);
1859 pidObj->SetOldPid(kTRUE);
1863 SetUseDefaultPID(kFALSE);
1878 //---------------------------------------------------------------------------
1879 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() {
1881 // Default 2010 PbPb cut object
1882 SetStandardCutsPbPb2010();
1883 AliAODPidHF *pidobj=GetPidHF();
1885 pidobj->SetOldPid(kFALSE);
1888 // Enable all 2011 PbPb run triggers
1890 SetTriggerClass("");
1891 ResetMaskAndEnableMBTrigger();
1892 EnableCentralTrigger();
1893 EnableSemiCentralTrigger();
1897 //---------------------------------------------------------------------------
1898 Int_t AliRDHFCutsD0toKpi::IsSelectedCombPID(AliAODRecoDecayHF* d)
1900 // ############################################################
1902 // Apply Bayesian PID selection
1903 // Implementation of Bayesian PID: Jeremy Wilkinson
1905 // ############################################################
1907 if (!fUsePID || !d) return 3;
1910 if (fBayesianStrategy == kBayesWeightNoFilter) {
1911 //WeightNoFilter: Accept all particles (no PID cut) but fill mass histos with weights in task
1912 CalculateBayesianWeights(d);
1920 Int_t returnvalue = 0;
1922 Int_t isPosKaon = 0, isNegKaon = 0, isPosPion = 0, isNegPion = 0;
1924 //Bayesian methods used here check for ID of kaon, and whether it is positive or negative.
1927 Bool_t checkPIDInfo[2] = {kTRUE, kTRUE};
1928 AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
1930 if ((aodtrack[0]->Charge()*aodtrack[1]->Charge()) != -1) {
1931 return 0; //reject if charges not opposite
1933 Double_t momentumpositive=0., momentumnegative=0.; //Used in "prob > prior" method
1934 for (Int_t daught = 0; daught < 2; daught++) {
1937 if (aodtrack[daught]->Charge() == -1) {
1938 momentumnegative = aodtrack[daught]->P();
1940 if (aodtrack[daught]->Charge() == +1) {
1941 momentumpositive = aodtrack[daught]->P();
1943 if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
1944 checkPIDInfo[daught] = kFALSE;
1950 //Loop over daughters ends here
1952 if (!checkPIDInfo[0] && !checkPIDInfo[1]) {
1953 return 0; //Reject if both daughters lack both TPC and TOF info
1957 CalculateBayesianWeights(d); //Calculates all Bayesian probabilities for both positive and negative tracks
1958 // Double_t prob0[AliPID::kSPECIES];
1959 // fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), prob0);
1960 ///Three possible Bayesian probability cuts: Picked using SetBayesianCondition(int).
1961 switch (fBayesianCondition) {
1962 ///A: Standard max. probability method (accept most likely species)
1964 if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kKaon]) { //If highest probability lies with kaon
1965 isPosKaon = 1; //flag [daught] as a kaon
1968 if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kPion]) { //If highest probability lies with pion
1969 isPosPion = 1; //flag [daught] as a pion
1972 if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kKaon]) { //If highest probability lies with kaon
1973 isNegKaon = 1; //flag [daught] as a kaon
1976 if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kPion]) { //If highest probability lies with kaon
1977 isNegPion = 1; //flag [daught] as a pion
1982 ///B: Accept if probability greater than prior
1985 if (fWeightsNegative[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
1986 GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumnegative)))) { //Retrieves relevant prior, gets value at momentum
1989 if (fWeightsPositive[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
1990 GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumpositive)))) { //Retrieves relevant prior, gets value at momentum
1996 ///C: Accept if probability greater than user-defined threshold
1998 if (fWeightsNegative[AliPID::kKaon] > fProbThreshold) {
2001 if (fWeightsNegative[AliPID::kPion] > fProbThreshold) {
2005 if (fWeightsPositive[AliPID::kKaon] > fProbThreshold) {
2009 if (fWeightsPositive[AliPID::kPion] > fProbThreshold) {
2017 //Momentum-based selection (also applied to filtered weighted method)
2019 if (fBayesianStrategy == kBayesMomentum || fBayesianCondition == kBayesWeight) {
2020 if (isNegKaon && isPosKaon) { // If both are kaons, reject
2023 } else if (isNegKaon && isPosPion) { //If negative kaon present, D0
2025 } else if (isPosKaon && isNegPion) { //If positive kaon present, D0bar
2027 } else { //If neither ID'd as kaon, subject to extra tests
2030 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
2031 if (aodtrack[0]->Charge() == -1) {
2032 isD0 = 0; //Check charge and reject based on pion hypothesis
2034 if (aodtrack[0]->Charge() == 1) {
2038 if (aodtrack[1]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[1], "TOF"))) {
2039 if (aodtrack[1]->Charge() == -1) {
2042 if (aodtrack[1]->Charge() == 1) {
2050 if (isD0 && isD0bar) {
2053 if (isD0 && !isD0bar) {
2056 if (!isD0 && isD0bar) {
2059 if (!isD0 && !isD0bar) {
2065 if (fBayesianStrategy == kBayesSimple) {
2067 if (isPosKaon && isNegKaon) { //If both are ID'd as kaons, accept as possible
2069 } else if (isNegKaon && isPosPion) { //If negative kaon, D0
2071 } else if (isPosKaon && isNegPion) { //If positive kaon, D0-bar
2073 } else if (isPosPion && isNegPion) {
2074 returnvalue = 0; //If neither kaon, reject
2075 } else {returnvalue = 0;} //default
2087 //---------------------------------------------------------------------------
2088 void AliRDHFCutsD0toKpi::CalculateBayesianWeights(AliAODRecoDecayHF* d)
2091 //Function to compute weights for Bayesian method
2093 AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
2094 if ((aodtrack[0]->Charge() * aodtrack[1]->Charge()) != -1) {
2095 return; //Reject if charges do not oppose
2097 for (Int_t daught = 0; daught < 2; daught++) {
2100 if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
2105 // identify kaon, define weights
2106 if (aodtrack[daught]->Charge() == +1) {
2107 fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsPositive);
2110 if (aodtrack[daught]->Charge() == -1) {
2111 fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsNegative);