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"
38 ClassImp(AliRDHFCutsD0toKpi)
40 //--------------------------------------------------------------------------
41 AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const char* name) :
43 fUseSpecialCuts(kFALSE),
49 // Default Constructor
53 TString varNames[9]={"inv. mass [GeV]",
62 Bool_t isUpperCut[9]={kTRUE,
71 SetVarNames(nvars,varNames,isUpperCut);
72 Bool_t forOpt[9]={kFALSE,
81 SetVarsForOpt(4,forOpt);
82 Float_t limits[2]={0,999999999.};
86 //--------------------------------------------------------------------------
87 AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :
89 fUseSpecialCuts(source.fUseSpecialCuts),
90 fLowPt(source.fLowPt),
91 fDefaultPID(source.fDefaultPID),
99 //--------------------------------------------------------------------------
100 AliRDHFCutsD0toKpi &AliRDHFCutsD0toKpi::operator=(const AliRDHFCutsD0toKpi &source)
103 // assignment operator
105 if(&source == this) return *this;
107 AliRDHFCuts::operator=(source);
108 fUseSpecialCuts=source.fUseSpecialCuts;
109 fLowPt=source.fLowPt;
110 fDefaultPID=source.fDefaultPID;
116 //---------------------------------------------------------------------------
117 void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
119 // Fills in vars the values of the variables
122 if(nvars!=fnVarsForOpt) {
123 printf("AliRDHFCutsD0toKpi::GetCutsVarsForOpt: wrong number of variables\n");
127 AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)d;
132 if(TMath::Abs(pdgdaughters[0])==211) {
133 vars[iter]=dd->InvMassD0();
135 vars[iter]=dd->InvMassD0bar();
140 vars[iter]=dd->GetDCA();
144 if(TMath::Abs(pdgdaughters[0])==211) {
145 vars[iter] = dd->CosThetaStarD0();
147 vars[iter] = dd->CosThetaStarD0bar();
152 if(TMath::Abs(pdgdaughters[0])==321) {
153 vars[iter]=dd->PtProng(0);
156 vars[iter]=dd->PtProng(1);
161 if(TMath::Abs(pdgdaughters[0])==211) {
162 vars[iter]=dd->PtProng(0);
165 vars[iter]=dd->PtProng(1);
170 if(TMath::Abs(pdgdaughters[0])==321) {
171 vars[iter]=dd->Getd0Prong(0);
174 vars[iter]=dd->Getd0Prong(1);
179 if(TMath::Abs(pdgdaughters[0])==211) {
180 vars[iter]=dd->Getd0Prong(0);
183 vars[iter]=dd->Getd0Prong(1);
188 vars[iter]= dd->Prodd0d0();
192 vars[iter]=dd->CosPointingAngle();
197 //---------------------------------------------------------------------------
198 Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod) {
208 cout<<"Cut matrice not inizialized. Exit..."<<endl;
212 AliAODRecoDecayHF2Prong* d=(AliAODRecoDecayHF2Prong*)obj;
215 cout<<"AliAODRecoDecayHF2Prong null"<<endl;
219 Double_t ptD=d->Pt();
220 if(ptD<fMinPtCand) return 0;
221 if(ptD>fMaxPtCand) return 0;
223 // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
224 Int_t returnvaluePID=3;
225 Int_t returnvalueCuts=3;
227 // selection on candidate
228 if(selectionLevel==AliRDHFCuts::kAll ||
229 selectionLevel==AliRDHFCuts::kCandidate) {
233 //recalculate vertex w/o daughters
234 AliAODVertex *origownvtx=0x0;
235 if(fRemoveDaughtersFromPrimary && !fUseMCVertex) {
236 if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
237 if(!RecalcOwnPrimaryVtx(d,aod)) {
238 CleanOwnPrimaryVtx(d,aod,origownvtx);
243 if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
244 if(!SetMCPrimaryVtx(d,aod)) {
245 CleanOwnPrimaryVtx(d,aod,origownvtx);
252 Int_t okD0=0,okD0bar=0;
254 Int_t ptbin=PtBin(pt);
256 CleanOwnPrimaryVtx(d,aod,origownvtx);
260 Double_t mD0,mD0bar,ctsD0,ctsD0bar;
263 Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
265 d->InvMassD0(mD0,mD0bar);
266 if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
267 if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
268 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
270 if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
273 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;
274 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;
275 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
278 if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
279 TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
280 if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
281 TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
282 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
284 if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
287 d->CosThetaStarD0(ctsD0,ctsD0bar);
288 if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
289 if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
290 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
292 if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
294 if (returnvalueCuts!=0) {
295 if (okD0) returnvalueCuts=1; //cuts passed as D0
296 if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar
297 if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
302 if(fUseSpecialCuts) special=IsSelectedSpecialCuts(d);
303 if(!special) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
305 // unset recalculated primary vertex when not needed any more
306 CleanOwnPrimaryVtx(d,aod,origownvtx);
309 // go to selection with Kalman vertexing, if requested
310 returnvalueCuts = IsSelectedKF(d,aod);
313 fIsSelectedCuts=returnvalueCuts;
314 if(!returnvalueCuts) return 0;
319 if(selectionLevel==AliRDHFCuts::kAll ||
320 selectionLevel==AliRDHFCuts::kCandidate ||
321 selectionLevel==AliRDHFCuts::kPID) {
322 returnvaluePID = IsSelectedPID(d);
323 fIsSelectedPID=returnvaluePID;
324 if(!returnvaluePID) return 0;
327 Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
329 if(!returnvalueComb) return 0;
331 // selection on daughter tracks
332 if(selectionLevel==AliRDHFCuts::kAll ||
333 selectionLevel==AliRDHFCuts::kTracks) {
334 if(!AreDaughtersSelected(d)) return 0;
337 // cout<<"Pid = "<<returnvaluePID<<endl;
338 return returnvalueComb;
341 //------------------------------------------------------------------------------------------
342 Int_t AliRDHFCutsD0toKpi::IsSelectedKF(AliAODRecoDecayHF2Prong *d,
343 AliAODEvent* aod) const {
345 // Apply selection using KF-vertexing
348 AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
349 AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
351 if(!track0 || !track1) {
352 cout<<"one or two D0 daughters missing!"<<endl;
356 // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
357 Int_t returnvalueCuts=3;
359 // candidate selection with AliKF
360 AliKFParticle::SetField(aod->GetMagneticField()); // set the magnetic field
362 Int_t okD0=0,okD0bar=0;
365 // convert tracks into AliKFParticles
367 AliKFParticle negPiKF(*track1,-211); // neg pion kandidate
368 AliKFParticle negKKF(*track1,-321); // neg kaon kandidate
369 AliKFParticle posPiKF(*track0,211); // pos pion kandidate
370 AliKFParticle posKKF(*track0,321); // pos kaon kandidate
372 // build D0 candidates
374 AliKFParticle d0c(negKKF,posPiKF); // D0 candidate
375 AliKFParticle ad0c(posKKF,negPiKF); // D0bar candidate
377 // create kf primary vertices
379 AliAODVertex *vtx1 = aod->GetPrimaryVertex();
380 AliKFVertex primVtx(*vtx1);
381 AliKFVertex aprimVtx(*vtx1);
383 if(primVtx.GetNContributors()<=0) okD0 = 0;
384 if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
385 if(!okD0 && !okD0bar) returnvalueCuts=0;
389 Double_t d0mass = d0c.GetMass();
390 Double_t ad0mass = ad0c.GetMass();
392 // calculate P of D0 and D0bar
393 Double_t d0P = d0c.GetP();
394 Double_t d0Px = d0c.GetPx();
395 Double_t d0Py = d0c.GetPy();
396 Double_t d0Pz = d0c.GetPz();
397 Double_t ad0P = ad0c.GetP();
398 Double_t ad0Px = ad0c.GetPx();
399 Double_t ad0Py = ad0c.GetPy();
400 Double_t ad0Pz = ad0c.GetPz();
402 //calculate Pt of D0 and D0bar
404 Double_t pt=d0c.GetPt();
405 Double_t apt=ad0c.GetPt();
407 // remove D0 daughters from primary vertices (if used in vertex fit) and add D0-candidates
409 if(track0->GetUsedForPrimVtxFit()) {
414 if(track1->GetUsedForPrimVtxFit()) {
422 if(primVtx.GetNContributors()<=0) okD0 = 0;
423 if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
424 if(!okD0 && !okD0bar) returnvalueCuts=0;
426 //calculate cut variables
428 // calculate impact params of daughters w.r.t recalculated vertices
430 Double_t impactPi = posPiKF.GetDistanceFromVertexXY(primVtx);
431 Double_t aimpactPi = negPiKF.GetDistanceFromVertexXY(aprimVtx);
432 Double_t impactKa = negKKF.GetDistanceFromVertexXY(primVtx);
433 Double_t aimpactKa = posKKF.GetDistanceFromVertexXY(aprimVtx);
435 // calculate Product of Impact Params
437 Double_t prodParam = impactPi*impactKa;
438 Double_t aprodParam = aimpactPi*aimpactKa;
440 // calculate cosine of pointing angles
442 TVector3 mom(d0c.GetPx(),d0c.GetPy(),d0c.GetPz());
443 TVector3 fline(d0c.GetX()-primVtx.GetX(),
444 d0c.GetY()-primVtx.GetY(),
445 d0c.GetZ()-primVtx.GetZ());
446 Double_t pta = mom.Angle(fline);
447 Double_t cosP = TMath::Cos(pta); // cosine of pta for D0 candidate
449 TVector3 amom(ad0c.GetPx(),ad0c.GetPy(),ad0c.GetPz());
450 TVector3 afline(ad0c.GetX()-aprimVtx.GetX(),
451 ad0c.GetY()-aprimVtx.GetY(),
452 ad0c.GetZ()-aprimVtx.GetZ());
453 Double_t apta = amom.Angle(afline);
454 Double_t acosP = TMath::Cos(apta); // cosine of pta for D0bar candidate
456 // calculate P of Pions at Decay Position of D0 and D0bar candidates
457 negKKF.TransportToParticle(d0c);
458 posPiKF.TransportToParticle(d0c);
459 posKKF.TransportToParticle(ad0c);
460 negPiKF.TransportToParticle(ad0c);
462 Double_t pxPi = posPiKF.GetPx();
463 Double_t pyPi = posPiKF.GetPy();
464 Double_t pzPi = posPiKF.GetPz();
465 Double_t ptPi = posPiKF.GetPt();
467 Double_t apxPi = negPiKF.GetPx();
468 Double_t apyPi = negPiKF.GetPy();
469 Double_t apzPi = negPiKF.GetPz();
470 Double_t aptPi = negPiKF.GetPt();
472 // calculate Pt of Kaons at Decay Position of D0 and D0bar candidates
474 Double_t ptK = negKKF.GetPt();
475 Double_t aptK = posKKF.GetPt();
477 //calculate cos(thetastar)
478 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(421)->Mass();
480 massp[0] = TDatabasePDG::Instance()->GetParticle(321)->Mass();
481 massp[1] = TDatabasePDG::Instance()->GetParticle(211)->Mass();
482 Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)
483 -4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
485 // cos(thetastar) for D0 and Pion
487 Double_t d0E = TMath::Sqrt(massvtx*massvtx + d0P*d0P);
488 Double_t beta = d0P/d0E;
489 Double_t gamma = d0E/massvtx;
490 TVector3 momPi(pxPi,pyPi,pzPi);
491 TVector3 momTot(d0Px,d0Py,d0Pz);
492 Double_t q1 = momPi.Dot(momTot)/momTot.Mag();
493 Double_t cts = (q1/gamma-beta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
495 // cos(thetastar) for D0bar and Pion
497 Double_t ad0E = TMath::Sqrt(massvtx*massvtx + ad0P*ad0P);
498 Double_t abeta = ad0P/ad0E;
499 Double_t agamma = ad0E/massvtx;
500 TVector3 amomPi(apxPi,apyPi,apzPi);
501 TVector3 amomTot(ad0Px,ad0Py,ad0Pz);
502 Double_t aq1 = amomPi.Dot(amomTot)/amomTot.Mag();
503 Double_t acts = (aq1/agamma-abeta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
505 // calculate reduced Chi2 for the full D0 fit
506 d0c.SetProductionVertex(primVtx);
507 ad0c.SetProductionVertex(aprimVtx);
508 negKKF.SetProductionVertex(d0c);
509 posPiKF.SetProductionVertex(d0c);
510 posKKF.SetProductionVertex(ad0c);
511 negPiKF.SetProductionVertex(ad0c);
512 d0c.TransportToProductionVertex();
513 ad0c.TransportToProductionVertex();
515 // calculate the decay length
516 Double_t decayLengthD0 = d0c.GetDecayLength();
517 Double_t adecayLengthD0 = ad0c.GetDecayLength();
519 Double_t chi2D0 = 50.;
520 if(d0c.GetNDF() > 0 && d0c.GetChi2() >= 0) {
521 chi2D0 = d0c.GetChi2()/d0c.GetNDF();
524 Double_t achi2D0 = 50.;
525 if(ad0c.GetNDF() > 0 && ad0c.GetChi2() >= 0) {
526 achi2D0 = ad0c.GetChi2()/ad0c.GetNDF();
530 Int_t ptbin=PtBin(pt);
531 Int_t aptbin=PtBin(apt);
533 if(ptbin < 0) okD0 = 0;
534 if(aptbin < 0) okD0bar = 0;
535 if(!okD0 && !okD0bar) returnvalueCuts=0;
537 if(ptK < fCutsRD[GetGlobalIndex(3,ptbin)] || ptPi < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
538 if(aptK < fCutsRD[GetGlobalIndex(3,aptbin)] || aptPi < fCutsRD[GetGlobalIndex(4,aptbin)]) okD0bar = 0;
539 if(!okD0 && !okD0bar) returnvalueCuts=0;
542 if(TMath::Abs(impactKa) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
543 TMath::Abs(impactPi) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
545 if(TMath::Abs(aimpactKa) > fCutsRD[GetGlobalIndex(5,aptbin)] ||
546 TMath::Abs(aimpactPi) > fCutsRD[GetGlobalIndex(6,aptbin)]) okD0bar = 0;
548 if(!okD0 && !okD0bar) returnvalueCuts=0;
550 // for the moment via the standard method due to bug in AliKF
551 if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) okD0 = 0;
552 if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,aptbin)]) okD0bar = 0;
553 if(!okD0 && !okD0bar) returnvalueCuts=0;
556 if(TMath::Abs(d0mass-massvtx) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
557 if(TMath::Abs(ad0mass-massvtx) > fCutsRD[GetGlobalIndex(0,aptbin)]) okD0bar = 0;
558 if(!okD0 && !okD0bar) returnvalueCuts=0;
561 if(TMath::Abs(cts) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
562 if(TMath::Abs(acts) > fCutsRD[GetGlobalIndex(2,aptbin)]) okD0bar = 0;
563 if(!okD0 && !okD0bar) returnvalueCuts=0;
565 if(prodParam > fCutsRD[GetGlobalIndex(7,ptbin)]) okD0 = 0;
566 if(aprodParam > fCutsRD[GetGlobalIndex(7,aptbin)]) okD0bar = 0;
567 if(!okD0 && !okD0bar) returnvalueCuts=0;
569 if(cosP < fCutsRD[GetGlobalIndex(8,ptbin)]) okD0 = 0;
570 if(acosP < fCutsRD[GetGlobalIndex(8,aptbin)]) okD0bar = 0;
571 if(!okD0 && !okD0bar) returnvalueCuts=0;
573 if(chi2D0 > fCutsRD[GetGlobalIndex(10,ptbin)]) okD0 = 0;
574 if(achi2D0 > fCutsRD[GetGlobalIndex(10,aptbin)]) okD0bar = 0;
575 if(!okD0 && !okD0bar) returnvalueCuts=0;
577 if(decayLengthD0 < fCutsRD[GetGlobalIndex(9,ptbin)]) okD0 = 0;
578 if(adecayLengthD0 < fCutsRD[GetGlobalIndex(9,aptbin)]) okD0bar = 0;
579 if(!okD0 && !okD0bar) returnvalueCuts=0;
581 if(returnvalueCuts!=0) {
582 if(okD0) returnvalueCuts=1; //cuts passed as D0
583 if(okD0bar) returnvalueCuts=2; //cuts passed as D0bar
584 if(okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
587 return returnvalueCuts;
590 //---------------------------------------------------------------------------
592 Bool_t AliRDHFCutsD0toKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
595 // Checking if D0 is in fiducial acceptance region
599 // applying cut for pt > 5 GeV
600 AliDebug(2,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt));
601 if (TMath::Abs(y) > 0.8){
605 // appliying smooth cut for pt < 5 GeV
606 Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
607 Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
608 AliDebug(2,Form("pt of D0 = %f (< 5), cutting according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
609 if (y < minFiducialY || y > maxFiducialY){
616 //---------------------------------------------------------------------------
617 Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d)
619 // ############################################################
621 // Apply PID selection
624 // ############################################################
626 if(!fUsePID) return 3;
627 if(fDefaultPID) return IsSelectedPIDdefault(d);
629 Int_t isD0D0barPID[2]={1,2};
630 Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
632 // values convention -1 = discarded
633 // 0 = not identified (but compatible) || No PID (->hasPID flag)
635 // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
636 // Initial hypothesis: unknwon (but compatible)
637 combinedPID[0][0]=0; // prima figlia, Kaon
638 combinedPID[0][1]=0; // prima figlia, pione
639 combinedPID[1][0]=0; // seconda figlia, Kaon
640 combinedPID[1][1]=0; // seconda figlia, pion
642 Bool_t checkPIDInfo[2]={kTRUE,kTRUE};
643 Double_t sigma_tmp[3]={fPidHF->GetSigma(0),fPidHF->GetSigma(1),fPidHF->GetSigma(2)};
644 for(Int_t daught=0;daught<2;daught++){
646 AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
648 if(!(fPidHF->CheckStatus(aodtrack,"TPC")) && !(fPidHF->CheckStatus(aodtrack,"TOF"))) {
649 checkPIDInfo[daught]=kFALSE;
654 combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
658 if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
659 combinedPID[daught][1]=0;
661 fPidHF->SetTOF(kFALSE);
662 combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
663 fPidHF->SetTOF(kTRUE);
664 fPidHF->SetCompat(kTRUE);
668 if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded
672 else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){
673 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
674 else isD0D0barPID[0]=0;// if K+ D0 excluded
676 else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){
680 else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
681 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
682 else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
684 else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
685 if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
686 else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
689 if(fLowPt && d->Pt()<2.){
690 Double_t sigmaTPC[3]={3.,2.,0.};
691 fPidHF->SetSigmaForTPC(sigmaTPC);
693 combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
695 Double_t ptProng=aodtrack->P();
698 fPidHF->SetCompat(kFALSE);
699 combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
700 fPidHF->SetCompat(kTRUE);
703 if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
704 combinedPID[daught][1]=0;
706 fPidHF->SetTOF(kFALSE);
707 Double_t sigmaTPCpi[3]={3.,3.,0.};
708 fPidHF->SetSigmaForTPC(sigmaTPCpi);
709 combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
710 fPidHF->SetTOF(kTRUE);
712 Bool_t isTPCpion=fPidHF->IsPionRaw(aodtrack,"TPC");
714 combinedPID[daught][1]=1;
716 combinedPID[daught][1]=-1;
722 fPidHF->SetSigmaForTPC(sigma_tmp);
723 }// END OF LOOP ON DAUGHTERS
725 if(!checkPIDInfo[0] && !checkPIDInfo[1]) {
726 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
731 // FURTHER PID REQUEST (both daughter info is needed)
732 if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
733 fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
734 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
739 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
740 if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
741 fWhyRejection=32;// reject cases where the Kaon is not identified
745 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
747 // cout<<"Why? "<<fWhyRejection<<endl;
748 return isD0D0barPID[0]+isD0D0barPID[1];
750 //---------------------------------------------------------------------------
751 Int_t AliRDHFCutsD0toKpi::IsSelectedPIDdefault(AliAODRecoDecayHF* d)
753 // ############################################################
755 // Apply PID selection
758 // temporary selection: PID AS USED FOR D0 by Andrea Rossi (up to 28/06/2010)
760 // d must be a AliAODRecoDecayHF2Prong object
761 // returns 0 if both D0 and D0bar are rejectecd
762 // 1 if D0 is accepted while D0bar is rejected
763 // 2 if D0bar is accepted while D0 is rejected
764 // 3 if both are accepted
765 // fWhyRejection variable (not returned for the moment, print it if needed)
766 // keeps some information on why a candidate has been
767 // rejected according to the following (unfriendly?) scheme
768 // if more rejection cases are considered interesting, just add numbers
770 // TO BE CONSIDERED WITH A GRAIN OF SALT (the order in which cut are applied is relevant)
771 // from 20 to 30: "detector" selection (PID acceptance)
774 // 28: no (TOF||TPC) pid information (no kTOFpid,kTOFout,kTIME,kTPCpid,...)
776 // from 30 to 40: PID selection
777 // 31: no Kaon compatible tracks found between daughters
778 // 32: no Kaon identified tracks found (strong sel. at low momenta)
779 // 33: both mass hypotheses are rejected
781 // ############################################################
783 if(!fUsePID) return 3;
785 Int_t isD0D0barPID[2]={1,2};
786 Double_t nsigmaTPCpi=-1., nsigmaTPCK=-1.; //used for TPC pid
787 Double_t tofSig,times[5];// used fot TOF pid
788 Int_t hasPID[2]={2,2};// flag to count how many detectors give PID info for the daughters
789 Int_t isKaonPionTOF[2][2],isKaonPionTPC[2][2];
790 Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
792 // values convention -1 = discarded
793 // 0 = not identified (but compatible) || No PID (->hasPID flag)
795 // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
796 // Initial hypothesis: unknwon (but compatible)
797 isKaonPionTOF[0][0]=0;
798 isKaonPionTOF[0][1]=0;
799 isKaonPionTOF[1][0]=0;
800 isKaonPionTOF[1][1]=0;
802 isKaonPionTPC[0][0]=0;
803 isKaonPionTPC[0][1]=0;
804 isKaonPionTPC[1][0]=0;
805 isKaonPionTPC[1][1]=0;
814 for(Int_t daught=0;daught<2;daught++){
817 // ########### Step 0- CHECKING minimal PID "ACCEPTANCE" ####################
819 AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
821 if(!(aodtrack->GetStatus()&AliESDtrack::kTPCrefit)){
825 if(!(aodtrack->GetStatus()&AliESDtrack::kITSrefit)){
830 AliAODPid *pid=aodtrack->GetDetPid();
837 // ########### Step 1- Check of TPC and TOF response ####################
839 Double_t ptrack=aodtrack->P();
840 //#################### TPC PID #######################
841 if (!(aodtrack->GetStatus()&AliESDtrack::kTPCpid )){
842 // NO TPC PID INFO FOR THIS TRACK
846 static AliTPCPIDResponse theTPCpid;
847 AliAODPid *pidObj = aodtrack->GetDetPid();
848 Double_t ptProng=pidObj->GetTPCmomentum();
849 nsigmaTPCpi = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kPion);
850 nsigmaTPCK = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kKaon);
853 if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
854 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
855 if(TMath::Abs(nsigmaTPCpi)<2.)isKaonPionTPC[daught][1]=1;
856 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
858 //else if(ptrack<.8){
860 if(TMath::Abs(nsigmaTPCK)<1.)isKaonPionTPC[daught][0]=1;
861 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
862 if(TMath::Abs(nsigmaTPCpi)<1.)isKaonPionTPC[daught][1]=1;
863 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
866 // if(nsigmaTPCK>-2.&&nsigmaTPCK<1.)isKaonPionTPC[daught][0]=1;
867 if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
868 //if(nsigmaTPCpi>-1.&&nsigmaTPCpi<2.)isKaonPionTPC[daught][1]=1;
869 if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
874 // ##### TOF PID: do not ask nothing for pion/protons ############
875 if(!((aodtrack->GetStatus()&AliESDtrack::kTOFpid)&&(aodtrack->GetStatus()&AliESDtrack::kTOFout)&&(aodtrack->GetStatus()&AliESDtrack::kTIME))){
876 // NO TOF PID INFO FOR THIS TRACK
880 tofSig=pid->GetTOFsignal();
881 pid->GetIntegratedTimes(times);
882 if((tofSig-times[3])>5.*160.)return 0;// PROTON REJECTION
883 if(TMath::Abs(tofSig-times[3])>3.*160.){
884 isKaonPionTOF[daught][0]=-1;
888 isKaonPionTOF[daught][0]=1;
893 //######### Step 2: COMBINE TOF and TPC PID ###############
894 // we apply the following convention: if TPC and TOF disagree (discarded Vs identified) -> unknown
895 combinedPID[daught][0]=isKaonPionTOF[daught][0]+isKaonPionTPC[daught][0];
896 combinedPID[daught][1]=isKaonPionTOF[daught][1]+isKaonPionTPC[daught][1];
899 //######### Step 3: USE PID INFO
901 if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){// if not a K- and not a pi- both D0 and D0bar excluded
905 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
906 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
907 else isD0D0barPID[0]=0;// if K+ D0 excluded
909 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
913 else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
914 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
915 else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
917 else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
918 if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
919 else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
922 // ########## ALSO DIFFERENT TPC PID REQUEST FOR LOW pt D0: request of K identification ###############################
923 // ########## more tolerant criteria for single particle ID-> more selective criteria for D0 ##############################
924 // ############### NOT OPTIMIZED YET ###################################
926 isKaonPionTPC[daught][0]=0;
927 isKaonPionTPC[daught][1]=0;
928 AliAODPid *pidObj = aodtrack->GetDetPid();
929 Double_t ptProng=pidObj->GetTPCmomentum();
932 if(TMath::Abs(nsigmaTPCK)<3.)isKaonPionTPC[daught][0]=1;
933 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
934 if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
935 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
937 //else if(ptrack<.8){
939 if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
940 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
941 if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
942 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
945 if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
946 if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
950 }// END OF LOOP ON DAUGHTERS
952 // FURTHER PID REQUEST (both daughter info is needed)
953 if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
954 fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
957 else if(hasPID[0]==0&&hasPID[1]==0){
958 fWhyRejection=28;// reject cases in which no PID info is available
962 // request of K identification at low D0 pt
968 combinedPID[0][0]=isKaonPionTOF[0][0]+isKaonPionTPC[0][0];
969 combinedPID[0][1]=isKaonPionTOF[0][1]+isKaonPionTPC[0][1];
970 combinedPID[1][0]=isKaonPionTOF[1][0]+isKaonPionTPC[1][0];
971 combinedPID[1][1]=isKaonPionTOF[1][1]+isKaonPionTPC[1][1];
973 if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
974 fWhyRejection=32;// reject cases where the Kaon is not identified
979 // cout<<"Why? "<<fWhyRejection<<endl;
980 return isD0D0barPID[0]+isD0D0barPID[1];
985 //---------------------------------------------------------------------------
986 Int_t AliRDHFCutsD0toKpi::CombineSelectionLevels(Int_t selectionvalTrack,
987 Int_t selectionvalCand,
988 Int_t selectionvalPID) const
991 // This method combines the tracks, PID and cuts selection results
993 if(selectionvalTrack==0) return 0;
997 switch(selectionvalPID) {
1002 returnvalue=((selectionvalCand==1 || selectionvalCand==3) ? 1 : 0);
1005 returnvalue=((selectionvalCand==2 || selectionvalCand==3) ? 2 : 0);
1008 returnvalue=selectionvalCand;
1017 //----------------------------------------------------------------------------
1018 Int_t AliRDHFCutsD0toKpi::IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
1021 // Note: this method is temporary
1022 // Additional cuts on decay lenght and lower cut for d0 norm are applied using vertex without candidate's daughters
1027 Float_t normDecLengthCut=1.,decLengthCut=TMath::Min(d->P()*0.0066+0.01,0.06/*cm*/), normd0Cut=0.5;
1028 // "decay length" expo law with tau' = beta*gamma*ctau= p/m*ctau =p*0.0123/1.864~p*0.0066
1029 // decay lenght > ctau' implies to retain (1-1/e) (for signal without considering detector resolution),
1031 Int_t returnvalue=3; //cut passed
1032 for(Int_t i=0;i<2/*prongs*/;i++){
1033 if(TMath::Abs(d->Normalizedd0Prong(i))<normd0Cut) return 0; //normd0Cut not passed
1035 if(d->DecayLength2()<decLengthCut*decLengthCut) return 0; //decLengthCut not passed
1036 if(d->NormalizedDecayLength2()<normDecLengthCut*normDecLengthCut) return 0; //decLengthCut not passed
1042 //----------------------------------------------
1043 void AliRDHFCutsD0toKpi::SetUseKF(Bool_t useKF)
1045 //switch on candidate selection via AliKFparticle
1051 TString varNamesKF[11]={"inv. mass [GeV]",
1062 Bool_t isUpperCutKF[11]={kTRUE,
1073 SetVarNames(nvarsKF,varNamesKF,isUpperCutKF);
1075 Bool_t forOpt[11]={kFALSE,
1086 SetVarsForOpt(4,forOpt);
1092 void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {
1094 //STANDARD CUTS USED FOR 2010 pp analysis
1095 //dca cut will be enlarged soon to 400 micron
1098 SetName("D0toKpiCutsStandard");
1099 SetTitle("Standard Cuts for D0 analysis");
1101 // PILE UP REJECTION
1102 SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1108 // TRACKS ON SINGLE TRACKS
1109 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1110 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1111 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1112 esdTrackCuts->SetRequireITSRefit(kTRUE);
1113 // esdTrackCuts->SetMinNClustersITS(4);
1114 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1115 esdTrackCuts->SetMinDCAToVertexXY(0.);
1116 esdTrackCuts->SetEtaRange(-0.8,0.8);
1117 esdTrackCuts->SetPtRange(0.3,1.e10);
1119 AddTrackCuts(esdTrackCuts);
1121 const Int_t nptbins =13;
1122 const Double_t ptmax = 9999.;
1123 const Int_t nvars=9;
1124 Float_t ptbins[nptbins+1];
1140 SetGlobalIndex(nvars,nptbins);
1141 SetPtBins(nptbins+1,ptbins);
1143 Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,300.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.73},/* pt<0.5*/
1144 {0.400,300.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.73},/* 0.5<pt<1*/
1145 {0.400,200.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.75},/* 1<pt<2 */
1146 {0.400,200.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-15000.*1E-8,0.8},/* 2<pt<3 */
1147 {0.400,200.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85},/* 3<pt<4 */
1148 {0.400,200.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85},/* 4<pt<5 */
1149 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85},/* 5<pt<6 */
1150 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85},/* 6<pt<8 */
1151 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85},/* 8<pt<12 */
1152 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.*1E-8,0.85},/* 12<pt<16 */
1153 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.85},/* 16<pt<20 */
1154 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.85},/* 20<pt<24 */
1155 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.85}};/* pt>24 */
1158 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1159 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1160 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1162 for (Int_t ibin=0;ibin<nptbins;ibin++){
1163 for (Int_t ivar = 0; ivar<nvars; ivar++){
1164 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1168 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1169 SetUseSpecialCuts(kTRUE);
1170 SetRemoveDaughtersFromPrim(kTRUE);
1172 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1173 delete [] cutsMatrixTransposeStand;
1174 cutsMatrixTransposeStand=NULL;
1177 AliAODPidHF* pidObj=new AliAODPidHF();
1178 //pidObj->SetName("pid4D0");
1180 const Int_t nlims=2;
1181 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1182 Bool_t compat=kTRUE; //effective only for this mode
1184 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1185 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1186 pidObj->SetMatch(mode);
1187 pidObj->SetPLimit(plims,nlims);
1188 pidObj->SetSigma(sigmas);
1189 pidObj->SetCompat(compat);
1190 pidObj->SetTPC(kTRUE);
1191 pidObj->SetTOF(kTRUE);
1195 SetUseDefaultPID(kFALSE);
1208 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
1210 //PRELIMINARY CUTS USED FOR 2010 PbPb analysis
1214 SetName("D0toKpiCutsStandard");
1215 SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
1217 // PILE UP REJECTION
1218 //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1224 // TRACKS ON SINGLE TRACKS
1225 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1226 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1227 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1228 esdTrackCuts->SetRequireITSRefit(kTRUE);
1229 // esdTrackCuts->SetMinNClustersITS(4);
1230 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1231 esdTrackCuts->SetMinDCAToVertexXY(0.);
1232 esdTrackCuts->SetEtaRange(-0.8,0.8);
1233 esdTrackCuts->SetPtRange(0.3,1.e10);
1235 AddTrackCuts(esdTrackCuts);
1237 const Int_t nptbins =13;
1238 const Double_t ptmax = 9999.;
1239 const Int_t nvars=9;
1240 Float_t ptbins[nptbins+1];
1256 SetGlobalIndex(nvars,nptbins);
1257 SetPtBins(nptbins+1,ptbins);
1259 Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,300.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.8},/* pt<0.5*/
1260 {0.400,300.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.8},/* 0.5<pt<1*/
1261 {0.400,250.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-32000.*1E-8,0.8},/* 1<pt<2 */
1262 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-26000.*1E-8,0.94},/* 2<pt<3 */
1263 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-1500.*1E-8,0.88},/* 3<pt<4 */
1264 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-1500.*1E-8,0.88},/* 4<pt<5 */
1265 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.90},/* 5<pt<6 */
1266 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.90},/* 6<pt<8 */
1267 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.90},/* 8<pt<12 */
1268 {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.90},/* 12<pt<16 */
1269 {0.400,350.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85},/* 16<pt<20 */
1270 {0.400,350.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.*1E-8,0.85},/* 20<pt<24 */
1271 {0.400,350.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.*1E-8,0.82}};/* pt>24 */
1274 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1275 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1276 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1278 for (Int_t ibin=0;ibin<nptbins;ibin++){
1279 for (Int_t ivar = 0; ivar<nvars; ivar++){
1280 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1284 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1285 SetUseSpecialCuts(kTRUE);
1286 SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1287 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1288 delete [] cutsMatrixTransposeStand;
1289 cutsMatrixTransposeStand=NULL;
1292 AliAODPidHF* pidObj=new AliAODPidHF();
1293 //pidObj->SetName("pid4D0");
1295 const Int_t nlims=2;
1296 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1297 Bool_t compat=kTRUE; //effective only for this mode
1299 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1300 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1301 pidObj->SetMatch(mode);
1302 pidObj->SetPLimit(plims,nlims);
1303 pidObj->SetSigma(sigmas);
1304 pidObj->SetCompat(compat);
1305 pidObj->SetTPC(kTRUE);
1306 pidObj->SetTOF(kTRUE);
1310 SetUseDefaultPID(kTRUE);// TEMPORARY: PROTON EXCLUSION SET ONLY IN DEFAULT PID