]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliRDHFCutsD0toKpi.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added proton rejection (Andrea)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliRDHFCutsD0toKpi.cxx
CommitLineData
650b3ced 1/**************************************************************************
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
27de2dfb 16/* $Id$ */
17
650b3ced 18/////////////////////////////////////////////////////////////
19//
20// Class for cuts on AOD reconstructed D0->Kpi
21//
22// Author: A.Dainese, andrea.dainese@pd.infn.it
23/////////////////////////////////////////////////////////////
24
25#include <TDatabasePDG.h>
26#include <Riostream.h>
27
28#include "AliRDHFCutsD0toKpi.h"
29#include "AliAODRecoDecayHF2Prong.h"
30#include "AliAODTrack.h"
31#include "AliESDtrack.h"
937e1290 32#include "AliAODPid.h"
33#include "AliTPCPIDResponse.h"
34#include "AliAODVertex.h"
5f25117b 35#include "AliKFVertex.h"
36#include "AliKFParticle.h"
650b3ced 37
38ClassImp(AliRDHFCutsD0toKpi)
39
40//--------------------------------------------------------------------------
a9b75906 41AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const char* name) :
937e1290 42AliRDHFCuts(name),
08ea5c32 43fUseSpecialCuts(kFALSE),
44fLowPt(kTRUE),
5f25117b 45fDefaultPID(kTRUE),
46fUseKF(kFALSE)
650b3ced 47{
48 //
49 // Default Constructor
50 //
51 Int_t nvars=9;
52 SetNVars(nvars);
53 TString varNames[9]={"inv. mass [GeV]",
54 "dca [cm]",
55 "cosThetaStar",
56 "pTK [GeV/c]",
57 "pTPi [GeV/c]",
58 "d0K [cm]",
59 "d0Pi [cm]",
60 "d0d0 [cm^2]",
61 "cosThetaPoint"};
62 Bool_t isUpperCut[9]={kTRUE,
63 kTRUE,
64 kTRUE,
65 kFALSE,
66 kFALSE,
67 kTRUE,
68 kTRUE,
69 kTRUE,
70 kFALSE};
71 SetVarNames(nvars,varNames,isUpperCut);
72 Bool_t forOpt[9]={kFALSE,
73 kTRUE,
74 kTRUE,
75 kFALSE,
76 kFALSE,
77 kFALSE,
78 kFALSE,
79 kTRUE,
80 kTRUE};
81 SetVarsForOpt(4,forOpt);
82 Float_t limits[2]={0,999999999.};
83 SetPtBins(2,limits);
937e1290 84
650b3ced 85}
86//--------------------------------------------------------------------------
87AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :
08ea5c32 88 AliRDHFCuts(source),
89 fUseSpecialCuts(source.fUseSpecialCuts),
90 fLowPt(source.fLowPt),
5f25117b 91 fDefaultPID(source.fDefaultPID),
92 fUseKF(source.fUseKF)
650b3ced 93{
94 //
95 // Copy constructor
96 //
97
98}
99//--------------------------------------------------------------------------
100AliRDHFCutsD0toKpi &AliRDHFCutsD0toKpi::operator=(const AliRDHFCutsD0toKpi &source)
101{
102 //
103 // assignment operator
104 //
105 if(&source == this) return *this;
106
08ea5c32 107 AliRDHFCuts::operator=(source);
937e1290 108 fUseSpecialCuts=source.fUseSpecialCuts;
08ea5c32 109 fLowPt=source.fLowPt;
110 fDefaultPID=source.fDefaultPID;
650b3ced 111
112 return *this;
113}
114
115
116//---------------------------------------------------------------------------
117void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
118 //
119 // Fills in vars the values of the variables
120 //
121
650b3ced 122 if(nvars!=fnVarsForOpt) {
123 printf("AliRDHFCutsD0toKpi::GetCutsVarsForOpt: wrong number of variables\n");
124 return;
125 }
126
127 AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)d;
650b3ced 128
650b3ced 129 Int_t iter=-1;
130 if(fVarsForOpt[0]){
131 iter++;
132 if(TMath::Abs(pdgdaughters[0])==211) {
133 vars[iter]=dd->InvMassD0();
134 } else {
135 vars[iter]=dd->InvMassD0bar();
136 }
137 }
138 if(fVarsForOpt[1]){
139 iter++;
140 vars[iter]=dd->GetDCA();
141 }
142 if(fVarsForOpt[2]){
143 iter++;
144 if(TMath::Abs(pdgdaughters[0])==211) {
145 vars[iter] = dd->CosThetaStarD0();
146 } else {
147 vars[iter] = dd->CosThetaStarD0bar();
148 }
149 }
150 if(fVarsForOpt[3]){
151 iter++;
152 if(TMath::Abs(pdgdaughters[0])==321) {
153 vars[iter]=dd->PtProng(0);
154 }
155 else{
156 vars[iter]=dd->PtProng(1);
157 }
158 }
159 if(fVarsForOpt[4]){
160 iter++;
161 if(TMath::Abs(pdgdaughters[0])==211) {
162 vars[iter]=dd->PtProng(0);
163 }
164 else{
165 vars[iter]=dd->PtProng(1);
166 }
167 }
168 if(fVarsForOpt[5]){
169 iter++;
170 if(TMath::Abs(pdgdaughters[0])==321) {
171 vars[iter]=dd->Getd0Prong(0);
172 }
173 else{
174 vars[iter]=dd->Getd0Prong(1);
175 }
176 }
177 if(fVarsForOpt[6]){
178 iter++;
179 if(TMath::Abs(pdgdaughters[0])==211) {
180 vars[iter]=dd->Getd0Prong(0);
181 }
182 else{
183 vars[iter]=dd->Getd0Prong(1);
184 }
185 }
186 if(fVarsForOpt[7]){
187 iter++;
188 vars[iter]= dd->Prodd0d0();
189 }
190 if(fVarsForOpt[8]){
191 iter++;
192 vars[iter]=dd->CosPointingAngle();
193 }
194
195 return;
196}
197//---------------------------------------------------------------------------
937e1290 198Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod) {
650b3ced 199 //
200 // Apply selection
201 //
47aa3d55 202
203
d7688946 204 fIsSelectedCuts=0;
205 fIsSelectedPID=0;
206
650b3ced 207 if(!fCutsRD){
208 cout<<"Cut matrice not inizialized. Exit..."<<endl;
3cbc09cd 209 return 0;
650b3ced 210 }
211 //PrintAll();
3cbc09cd 212 AliAODRecoDecayHF2Prong* d=(AliAODRecoDecayHF2Prong*)obj;
650b3ced 213
650b3ced 214 if(!d){
215 cout<<"AliAODRecoDecayHF2Prong null"<<endl;
216 return 0;
217 }
218
47aa3d55 219 Double_t ptD=d->Pt();
220 if(ptD<fMinPtCand) return 0;
221 if(ptD>fMaxPtCand) return 0;
937e1290 222
c96634a2 223 // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
224 Int_t returnvaluePID=3;
225 Int_t returnvalueCuts=3;
c96634a2 226
650b3ced 227 // selection on candidate
228 if(selectionLevel==AliRDHFCuts::kAll ||
229 selectionLevel==AliRDHFCuts::kCandidate) {
937e1290 230
d7688946 231 if(!fUseKF) {
5f25117b 232
233 //recalculate vertex w/o daughters
234 AliAODVertex *origownvtx=0x0;
748e77de 235 if(fRemoveDaughtersFromPrimary && !fUseMCVertex) {
793fdbec 236 if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
237 if(!RecalcOwnPrimaryVtx(d,aod)) {
238 CleanOwnPrimaryVtx(d,aod,origownvtx);
239 return 0;
240 }
5f25117b 241 }
87baf03d 242
70841bf2 243 if(fUseMCVertex) {
793fdbec 244 if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
245 if(!SetMCPrimaryVtx(d,aod)) {
246 CleanOwnPrimaryVtx(d,aod,origownvtx);
247 return 0;
248 }
70841bf2 249 }
5f25117b 250
5f25117b 251 Double_t pt=d->Pt();
252
253 Int_t okD0=0,okD0bar=0;
254
255 Int_t ptbin=PtBin(pt);
256 if (ptbin==-1) {
793fdbec 257 CleanOwnPrimaryVtx(d,aod,origownvtx);
937e1290 258 return 0;
259 }
d7688946 260
5f25117b 261 Double_t mD0,mD0bar,ctsD0,ctsD0bar;
262 okD0=1; okD0bar=1;
263
264 Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
d7688946 265
266 d->InvMassD0(mD0,mD0bar);
267 if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
268 if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
793fdbec 269 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
d7688946 270
793fdbec 271 if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
d7688946 272
5f25117b 273
be9cad0c 274 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;
275 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;
793fdbec 276 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
5f25117b 277
278
279 if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
280 TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
281 if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
282 TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
793fdbec 283 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
5f25117b 284
793fdbec 285 if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
5f25117b 286
650b3ced 287
5f25117b 288 d->CosThetaStarD0(ctsD0,ctsD0bar);
289 if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
290 if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
793fdbec 291 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
5f25117b 292
793fdbec 293 if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
5f25117b 294
295 if (returnvalueCuts!=0) {
296 if (okD0) returnvalueCuts=1; //cuts passed as D0
297 if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar
298 if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
299 }
300
301 // call special cuts
302 Int_t special=1;
303 if(fUseSpecialCuts) special=IsSelectedSpecialCuts(d);
5a4e6774 304 if(!special) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
5f25117b 305
306 // unset recalculated primary vertex when not needed any more
793fdbec 307 CleanOwnPrimaryVtx(d,aod,origownvtx);
650b3ced 308
d7688946 309 } else {
310 // go to selection with Kalman vertexing, if requested
311 returnvalueCuts = IsSelectedKF(d,aod);
312 }
dcc2ade0 313
d7688946 314 fIsSelectedCuts=returnvalueCuts;
dcc2ade0 315 if(!returnvalueCuts) return 0;
5f25117b 316 }
317
650b3ced 318
dcc2ade0 319 // selection on PID
320 if(selectionLevel==AliRDHFCuts::kAll ||
321 selectionLevel==AliRDHFCuts::kCandidate ||
322 selectionLevel==AliRDHFCuts::kPID) {
323 returnvaluePID = IsSelectedPID(d);
d7688946 324 fIsSelectedPID=returnvaluePID;
dcc2ade0 325 if(!returnvaluePID) return 0;
326 }
327
328 Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
329
330 if(!returnvalueComb) return 0;
331
332 // selection on daughter tracks
333 if(selectionLevel==AliRDHFCuts::kAll ||
334 selectionLevel==AliRDHFCuts::kTracks) {
335 if(!AreDaughtersSelected(d)) return 0;
336 }
337
5f25117b 338 // cout<<"Pid = "<<returnvaluePID<<endl;
e2b51a19 339 return returnvalueComb;
5f25117b 340}
341
342//------------------------------------------------------------------------------------------
343Int_t AliRDHFCutsD0toKpi::IsSelectedKF(AliAODRecoDecayHF2Prong *d,
344 AliAODEvent* aod) const {
345 //
346 // Apply selection using KF-vertexing
347 //
348
349 AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
350 AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
351
352 if(!track0 || !track1) {
353 cout<<"one or two D0 daughters missing!"<<endl;
354 return 0;
355 }
356
357 // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
358 Int_t returnvalueCuts=3;
359
360 // candidate selection with AliKF
361 AliKFParticle::SetField(aod->GetMagneticField()); // set the magnetic field
937e1290 362
5f25117b 363 Int_t okD0=0,okD0bar=0;
364 okD0=1; okD0bar=1;
365
366 // convert tracks into AliKFParticles
367
368 AliKFParticle negPiKF(*track1,-211); // neg pion kandidate
369 AliKFParticle negKKF(*track1,-321); // neg kaon kandidate
370 AliKFParticle posPiKF(*track0,211); // pos pion kandidate
371 AliKFParticle posKKF(*track0,321); // pos kaon kandidate
372
373 // build D0 candidates
374
375 AliKFParticle d0c(negKKF,posPiKF); // D0 candidate
376 AliKFParticle ad0c(posKKF,negPiKF); // D0bar candidate
377
378 // create kf primary vertices
379
380 AliAODVertex *vtx1 = aod->GetPrimaryVertex();
381 AliKFVertex primVtx(*vtx1);
382 AliKFVertex aprimVtx(*vtx1);
383
384 if(primVtx.GetNContributors()<=0) okD0 = 0;
385 if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
386 if(!okD0 && !okD0bar) returnvalueCuts=0;
387
388 // calculate mass
389
390 Double_t d0mass = d0c.GetMass();
391 Double_t ad0mass = ad0c.GetMass();
392
393 // calculate P of D0 and D0bar
394 Double_t d0P = d0c.GetP();
395 Double_t d0Px = d0c.GetPx();
396 Double_t d0Py = d0c.GetPy();
397 Double_t d0Pz = d0c.GetPz();
398 Double_t ad0P = ad0c.GetP();
399 Double_t ad0Px = ad0c.GetPx();
400 Double_t ad0Py = ad0c.GetPy();
401 Double_t ad0Pz = ad0c.GetPz();
402
403 //calculate Pt of D0 and D0bar
404
405 Double_t pt=d0c.GetPt();
406 Double_t apt=ad0c.GetPt();
407
408 // remove D0 daughters from primary vertices (if used in vertex fit) and add D0-candidates
409
410 if(track0->GetUsedForPrimVtxFit()) {
411 primVtx -= posPiKF;
412 aprimVtx -= posKKF;
413 }
414
415 if(track1->GetUsedForPrimVtxFit()) {
416 primVtx -= negKKF;
417 aprimVtx -= negPiKF;
418 }
419
420 primVtx += d0c;
421 aprimVtx += ad0c;
422
423 if(primVtx.GetNContributors()<=0) okD0 = 0;
424 if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
425 if(!okD0 && !okD0bar) returnvalueCuts=0;
426
427 //calculate cut variables
428
429 // calculate impact params of daughters w.r.t recalculated vertices
430
431 Double_t impactPi = posPiKF.GetDistanceFromVertexXY(primVtx);
432 Double_t aimpactPi = negPiKF.GetDistanceFromVertexXY(aprimVtx);
433 Double_t impactKa = negKKF.GetDistanceFromVertexXY(primVtx);
434 Double_t aimpactKa = posKKF.GetDistanceFromVertexXY(aprimVtx);
435
436 // calculate Product of Impact Params
437
438 Double_t prodParam = impactPi*impactKa;
439 Double_t aprodParam = aimpactPi*aimpactKa;
440
441 // calculate cosine of pointing angles
442
443 TVector3 mom(d0c.GetPx(),d0c.GetPy(),d0c.GetPz());
444 TVector3 fline(d0c.GetX()-primVtx.GetX(),
445 d0c.GetY()-primVtx.GetY(),
446 d0c.GetZ()-primVtx.GetZ());
447 Double_t pta = mom.Angle(fline);
448 Double_t cosP = TMath::Cos(pta); // cosine of pta for D0 candidate
449
450 TVector3 amom(ad0c.GetPx(),ad0c.GetPy(),ad0c.GetPz());
451 TVector3 afline(ad0c.GetX()-aprimVtx.GetX(),
452 ad0c.GetY()-aprimVtx.GetY(),
453 ad0c.GetZ()-aprimVtx.GetZ());
454 Double_t apta = amom.Angle(afline);
455 Double_t acosP = TMath::Cos(apta); // cosine of pta for D0bar candidate
456
457 // calculate P of Pions at Decay Position of D0 and D0bar candidates
458 negKKF.TransportToParticle(d0c);
459 posPiKF.TransportToParticle(d0c);
460 posKKF.TransportToParticle(ad0c);
461 negPiKF.TransportToParticle(ad0c);
462
463 Double_t pxPi = posPiKF.GetPx();
464 Double_t pyPi = posPiKF.GetPy();
465 Double_t pzPi = posPiKF.GetPz();
466 Double_t ptPi = posPiKF.GetPt();
467
468 Double_t apxPi = negPiKF.GetPx();
469 Double_t apyPi = negPiKF.GetPy();
470 Double_t apzPi = negPiKF.GetPz();
471 Double_t aptPi = negPiKF.GetPt();
472
473 // calculate Pt of Kaons at Decay Position of D0 and D0bar candidates
474
475 Double_t ptK = negKKF.GetPt();
476 Double_t aptK = posKKF.GetPt();
477
478 //calculate cos(thetastar)
479 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(421)->Mass();
480 Double_t massp[2];
481 massp[0] = TDatabasePDG::Instance()->GetParticle(321)->Mass();
482 massp[1] = TDatabasePDG::Instance()->GetParticle(211)->Mass();
483 Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)
484 -4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
485
486 // cos(thetastar) for D0 and Pion
487
488 Double_t d0E = TMath::Sqrt(massvtx*massvtx + d0P*d0P);
489 Double_t beta = d0P/d0E;
490 Double_t gamma = d0E/massvtx;
491 TVector3 momPi(pxPi,pyPi,pzPi);
492 TVector3 momTot(d0Px,d0Py,d0Pz);
493 Double_t q1 = momPi.Dot(momTot)/momTot.Mag();
494 Double_t cts = (q1/gamma-beta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
495
496 // cos(thetastar) for D0bar and Pion
497
498 Double_t ad0E = TMath::Sqrt(massvtx*massvtx + ad0P*ad0P);
499 Double_t abeta = ad0P/ad0E;
500 Double_t agamma = ad0E/massvtx;
501 TVector3 amomPi(apxPi,apyPi,apzPi);
502 TVector3 amomTot(ad0Px,ad0Py,ad0Pz);
503 Double_t aq1 = amomPi.Dot(amomTot)/amomTot.Mag();
504 Double_t acts = (aq1/agamma-abeta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
505
506 // calculate reduced Chi2 for the full D0 fit
507 d0c.SetProductionVertex(primVtx);
508 ad0c.SetProductionVertex(aprimVtx);
509 negKKF.SetProductionVertex(d0c);
510 posPiKF.SetProductionVertex(d0c);
511 posKKF.SetProductionVertex(ad0c);
512 negPiKF.SetProductionVertex(ad0c);
513 d0c.TransportToProductionVertex();
514 ad0c.TransportToProductionVertex();
515
516 // calculate the decay length
517 Double_t decayLengthD0 = d0c.GetDecayLength();
518 Double_t adecayLengthD0 = ad0c.GetDecayLength();
519
520 Double_t chi2D0 = 50.;
521 if(d0c.GetNDF() > 0 && d0c.GetChi2() >= 0) {
522 chi2D0 = d0c.GetChi2()/d0c.GetNDF();
523 }
524
525 Double_t achi2D0 = 50.;
526 if(ad0c.GetNDF() > 0 && ad0c.GetChi2() >= 0) {
527 achi2D0 = ad0c.GetChi2()/ad0c.GetNDF();
528 }
529
530 // Get the Pt-bins
531 Int_t ptbin=PtBin(pt);
532 Int_t aptbin=PtBin(apt);
533
534 if(ptbin < 0) okD0 = 0;
535 if(aptbin < 0) okD0bar = 0;
536 if(!okD0 && !okD0bar) returnvalueCuts=0;
537
538 if(ptK < fCutsRD[GetGlobalIndex(3,ptbin)] || ptPi < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
539 if(aptK < fCutsRD[GetGlobalIndex(3,aptbin)] || aptPi < fCutsRD[GetGlobalIndex(4,aptbin)]) okD0bar = 0;
540 if(!okD0 && !okD0bar) returnvalueCuts=0;
541
542
543 if(TMath::Abs(impactKa) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
544 TMath::Abs(impactPi) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
545
546 if(TMath::Abs(aimpactKa) > fCutsRD[GetGlobalIndex(5,aptbin)] ||
547 TMath::Abs(aimpactPi) > fCutsRD[GetGlobalIndex(6,aptbin)]) okD0bar = 0;
548
549 if(!okD0 && !okD0bar) returnvalueCuts=0;
550
551 // for the moment via the standard method due to bug in AliKF
552 if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) okD0 = 0;
553 if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,aptbin)]) okD0bar = 0;
554 if(!okD0 && !okD0bar) returnvalueCuts=0;
650b3ced 555
650b3ced 556
5f25117b 557 if(TMath::Abs(d0mass-massvtx) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
558 if(TMath::Abs(ad0mass-massvtx) > fCutsRD[GetGlobalIndex(0,aptbin)]) okD0bar = 0;
559 if(!okD0 && !okD0bar) returnvalueCuts=0;
560
561
562 if(TMath::Abs(cts) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
563 if(TMath::Abs(acts) > fCutsRD[GetGlobalIndex(2,aptbin)]) okD0bar = 0;
564 if(!okD0 && !okD0bar) returnvalueCuts=0;
565
566 if(prodParam > fCutsRD[GetGlobalIndex(7,ptbin)]) okD0 = 0;
567 if(aprodParam > fCutsRD[GetGlobalIndex(7,aptbin)]) okD0bar = 0;
568 if(!okD0 && !okD0bar) returnvalueCuts=0;
650b3ced 569
5f25117b 570 if(cosP < fCutsRD[GetGlobalIndex(8,ptbin)]) okD0 = 0;
571 if(acosP < fCutsRD[GetGlobalIndex(8,aptbin)]) okD0bar = 0;
572 if(!okD0 && !okD0bar) returnvalueCuts=0;
573
574 if(chi2D0 > fCutsRD[GetGlobalIndex(10,ptbin)]) okD0 = 0;
575 if(achi2D0 > fCutsRD[GetGlobalIndex(10,aptbin)]) okD0bar = 0;
576 if(!okD0 && !okD0bar) returnvalueCuts=0;
577
578 if(decayLengthD0 < fCutsRD[GetGlobalIndex(9,ptbin)]) okD0 = 0;
579 if(adecayLengthD0 < fCutsRD[GetGlobalIndex(9,aptbin)]) okD0bar = 0;
580 if(!okD0 && !okD0bar) returnvalueCuts=0;
650b3ced 581
5f25117b 582 if(returnvalueCuts!=0) {
583 if(okD0) returnvalueCuts=1; //cuts passed as D0
584 if(okD0bar) returnvalueCuts=2; //cuts passed as D0bar
585 if(okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
650b3ced 586 }
587
5f25117b 588 return returnvalueCuts;
650b3ced 589}
5f25117b 590
650b3ced 591//---------------------------------------------------------------------------
5f25117b 592
aa8d25fe 593Bool_t AliRDHFCutsD0toKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
594{
595 //
596 // Checking if D0 is in fiducial acceptance region
597 //
598
c96634a2 599 if(pt > 5.) {
aa8d25fe 600 // applying cut for pt > 5 GeV
937e1290 601 AliDebug(2,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt));
aa8d25fe 602 if (TMath::Abs(y) > 0.8){
603 return kFALSE;
604 }
605 } else {
606 // appliying smooth cut for pt < 5 GeV
607 Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
608 Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
937e1290 609 AliDebug(2,Form("pt of D0 = %f (< 5), cutting according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
aa8d25fe 610 if (y < minFiducialY || y > maxFiducialY){
611 return kFALSE;
612 }
613 }
c96634a2 614
aa8d25fe 615 return kTRUE;
aa8d25fe 616}
c96634a2 617//---------------------------------------------------------------------------
937e1290 618Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d)
08ea5c32 619{
620 // ############################################################
621 //
622 // Apply PID selection
623 //
624 //
625 // ############################################################
626
627 if(!fUsePID) return 3;
628 if(fDefaultPID) return IsSelectedPIDdefault(d);
629 fWhyRejection=0;
630 Int_t isD0D0barPID[2]={1,2};
631 Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
632 // same for prong 2
633 // values convention -1 = discarded
634 // 0 = not identified (but compatible) || No PID (->hasPID flag)
635 // 1 = identified
636 // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
637 // Initial hypothesis: unknwon (but compatible)
638 combinedPID[0][0]=0; // prima figlia, Kaon
639 combinedPID[0][1]=0; // prima figlia, pione
640 combinedPID[1][0]=0; // seconda figlia, Kaon
641 combinedPID[1][1]=0; // seconda figlia, pion
642
643 Bool_t checkPIDInfo[2]={kTRUE,kTRUE};
644 Double_t sigma_tmp[3]={fPidHF->GetSigma(0),fPidHF->GetSigma(1),fPidHF->GetSigma(2)};
645 for(Int_t daught=0;daught<2;daught++){
646 //Loop con prongs
647 AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
87baf03d 648 if(fPidHF->IsTOFPiKexcluded(aodtrack,5.)) return 0;
08ea5c32 649
650 if(!(fPidHF->CheckStatus(aodtrack,"TPC")) && !(fPidHF->CheckStatus(aodtrack,"TOF"))) {
87baf03d 651 checkPIDInfo[daught]=kFALSE;
652 continue;
08ea5c32 653 }
654
655 // identify kaon
656 combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
657
658 // identify pion
659
660 if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
661 combinedPID[daught][1]=0;
662 }else{
663 fPidHF->SetTOF(kFALSE);
664 combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
665 fPidHF->SetTOF(kTRUE);
666 fPidHF->SetCompat(kTRUE);
667 }
668
669
670 if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded
671 isD0D0barPID[0]=0;
672 isD0D0barPID[1]=0;
673 }
674 else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){
675 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
676 else isD0D0barPID[0]=0;// if K+ D0 excluded
677 }
678 else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){
679 isD0D0barPID[0]=0;
680 isD0D0barPID[1]=0;
681 }
682 else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
683 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
684 else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
685 }
686 else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
687 if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
688 else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
689 }
690
691 if(fLowPt && d->Pt()<2.){
692 Double_t sigmaTPC[3]={3.,2.,0.};
693 fPidHF->SetSigmaForTPC(sigmaTPC);
694 // identify kaon
695 combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
696
697 Double_t ptProng=aodtrack->P();
698
699 if(ptProng<0.6){
700 fPidHF->SetCompat(kFALSE);
701 combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
702 fPidHF->SetCompat(kTRUE);
703 }
704
705 if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
706 combinedPID[daught][1]=0;
707 }else{
708 fPidHF->SetTOF(kFALSE);
709 Double_t sigmaTPCpi[3]={3.,3.,0.};
710 fPidHF->SetSigmaForTPC(sigmaTPCpi);
711 combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
712 fPidHF->SetTOF(kTRUE);
713 if(ptProng<0.8){
714 Bool_t isTPCpion=fPidHF->IsPionRaw(aodtrack,"TPC");
715 if(isTPCpion){
716 combinedPID[daught][1]=1;
717 }else{
718 combinedPID[daught][1]=-1;
719 }
720 }
721 }
722
723 }
724 fPidHF->SetSigmaForTPC(sigma_tmp);
725 }// END OF LOOP ON DAUGHTERS
726
727 if(!checkPIDInfo[0] && !checkPIDInfo[1]) {
728 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
729 return 0;
730 }
731
732
733 // FURTHER PID REQUEST (both daughter info is needed)
734 if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
735 fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
736 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
737 return 0;
738 }
739
740 if(d->Pt()<2.){
741 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
742 if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
743 fWhyRejection=32;// reject cases where the Kaon is not identified
744 return 0;
745 }
746 }
747 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
748
749 // cout<<"Why? "<<fWhyRejection<<endl;
750 return isD0D0barPID[0]+isD0D0barPID[1];
751}
752//---------------------------------------------------------------------------
753Int_t AliRDHFCutsD0toKpi::IsSelectedPIDdefault(AliAODRecoDecayHF* d)
937e1290 754{
755 // ############################################################
c96634a2 756 //
757 // Apply PID selection
c96634a2 758 //
937e1290 759 //
760 // temporary selection: PID AS USED FOR D0 by Andrea Rossi (up to 28/06/2010)
761 //
762 // d must be a AliAODRecoDecayHF2Prong object
763 // returns 0 if both D0 and D0bar are rejectecd
764 // 1 if D0 is accepted while D0bar is rejected
765 // 2 if D0bar is accepted while D0 is rejected
766 // 3 if both are accepted
767 // fWhyRejection variable (not returned for the moment, print it if needed)
768 // keeps some information on why a candidate has been
769 // rejected according to the following (unfriendly?) scheme
770 // if more rejection cases are considered interesting, just add numbers
771 //
772 // TO BE CONSIDERED WITH A GRAIN OF SALT (the order in which cut are applied is relevant)
773 // from 20 to 30: "detector" selection (PID acceptance)
774 // 26: TPC refit
775 // 27: ITS refit
776 // 28: no (TOF||TPC) pid information (no kTOFpid,kTOFout,kTIME,kTPCpid,...)
777 //
778 // from 30 to 40: PID selection
779 // 31: no Kaon compatible tracks found between daughters
780 // 32: no Kaon identified tracks found (strong sel. at low momenta)
781 // 33: both mass hypotheses are rejected
782 //
783 // ############################################################
784
785 if(!fUsePID) return 3;
786 fWhyRejection=0;
787 Int_t isD0D0barPID[2]={1,2};
788 Double_t nsigmaTPCpi=-1., nsigmaTPCK=-1.; //used for TPC pid
789 Double_t tofSig,times[5];// used fot TOF pid
790 Int_t hasPID[2]={2,2};// flag to count how many detectors give PID info for the daughters
791 Int_t isKaonPionTOF[2][2],isKaonPionTPC[2][2];
792 Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
793 // same for prong 2
794 // values convention -1 = discarded
795 // 0 = not identified (but compatible) || No PID (->hasPID flag)
796 // 1 = identified
797 // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
798 // Initial hypothesis: unknwon (but compatible)
799 isKaonPionTOF[0][0]=0;
800 isKaonPionTOF[0][1]=0;
801 isKaonPionTOF[1][0]=0;
802 isKaonPionTOF[1][1]=0;
803
804 isKaonPionTPC[0][0]=0;
805 isKaonPionTPC[0][1]=0;
806 isKaonPionTPC[1][0]=0;
807 isKaonPionTPC[1][1]=0;
808
809 combinedPID[0][0]=0;
810 combinedPID[0][1]=0;
811 combinedPID[1][0]=0;
812 combinedPID[1][1]=0;
813
814
815
816 for(Int_t daught=0;daught<2;daught++){
817 //Loop con prongs
818
819 // ########### Step 0- CHECKING minimal PID "ACCEPTANCE" ####################
820
821 AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
822
823 if(!(aodtrack->GetStatus()&AliESDtrack::kTPCrefit)){
824 fWhyRejection=26;
825 return 0;
826 }
827 if(!(aodtrack->GetStatus()&AliESDtrack::kITSrefit)){
828 fWhyRejection=27;
829 return 0;
830 }
831
832 AliAODPid *pid=aodtrack->GetDetPid();
833 if(!pid) {
834 //delete esdtrack;
835 hasPID[daught]--;
836 continue;
837 }
838
839 // ########### Step 1- Check of TPC and TOF response ####################
840
841 Double_t ptrack=aodtrack->P();
842 //#################### TPC PID #######################
843 if (!(aodtrack->GetStatus()&AliESDtrack::kTPCpid )){
844 // NO TPC PID INFO FOR THIS TRACK
845 hasPID[daught]--;
846 }
847 else {
848 static AliTPCPIDResponse theTPCpid;
08ea5c32 849 AliAODPid *pidObj = aodtrack->GetDetPid();
850 Double_t ptProng=pidObj->GetTPCmomentum();
851 nsigmaTPCpi = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kPion);
852 nsigmaTPCK = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kKaon);
853 //if(ptrack<0.6){
854 if(ptProng<0.6){
937e1290 855 if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
856 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
857 if(TMath::Abs(nsigmaTPCpi)<2.)isKaonPionTPC[daught][1]=1;
858 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
859 }
08ea5c32 860 //else if(ptrack<.8){
861 else if(ptProng<.8){
937e1290 862 if(TMath::Abs(nsigmaTPCK)<1.)isKaonPionTPC[daught][0]=1;
863 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
864 if(TMath::Abs(nsigmaTPCpi)<1.)isKaonPionTPC[daught][1]=1;
865 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
866 }
867 else {
868 // if(nsigmaTPCK>-2.&&nsigmaTPCK<1.)isKaonPionTPC[daught][0]=1;
869 if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
870 //if(nsigmaTPCpi>-1.&&nsigmaTPCpi<2.)isKaonPionTPC[daught][1]=1;
871 if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
872 }
873 }
874
875
876 // ##### TOF PID: do not ask nothing for pion/protons ############
877 if(!((aodtrack->GetStatus()&AliESDtrack::kTOFpid)&&(aodtrack->GetStatus()&AliESDtrack::kTOFout)&&(aodtrack->GetStatus()&AliESDtrack::kTIME))){
878 // NO TOF PID INFO FOR THIS TRACK
879 hasPID[daught]--;
880 }
881 else{
882 tofSig=pid->GetTOFsignal();
883 pid->GetIntegratedTimes(times);
b14fe7f1 884 if((tofSig-times[3])>5.*160.)return 0;// PROTON REJECTION
937e1290 885 if(TMath::Abs(tofSig-times[3])>3.*160.){
886 isKaonPionTOF[daught][0]=-1;
887 }
888 else {
889 if(ptrack<1.5){
890 isKaonPionTOF[daught][0]=1;
891 }
892 }
893 }
894
895 //######### Step 2: COMBINE TOF and TPC PID ###############
896 // we apply the following convention: if TPC and TOF disagree (discarded Vs identified) -> unknown
897 combinedPID[daught][0]=isKaonPionTOF[daught][0]+isKaonPionTPC[daught][0];
898 combinedPID[daught][1]=isKaonPionTOF[daught][1]+isKaonPionTPC[daught][1];
899
900
901 //######### Step 3: USE PID INFO
902
903 if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){// if not a K- and not a pi- both D0 and D0bar excluded
904 isD0D0barPID[0]=0;
905 isD0D0barPID[1]=0;
906 }
907 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
908 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
909 else isD0D0barPID[0]=0;// if K+ D0 excluded
910 }
911 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
912 isD0D0barPID[0]=0;
913 isD0D0barPID[1]=0;
914 }
915 else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
916 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
917 else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
918 }
919 else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
920 if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
921 else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
922 }
923
924 // ########## ALSO DIFFERENT TPC PID REQUEST FOR LOW pt D0: request of K identification ###############################
925 // ########## more tolerant criteria for single particle ID-> more selective criteria for D0 ##############################
926 // ############### NOT OPTIMIZED YET ###################################
927 if(d->Pt()<2.){
928 isKaonPionTPC[daught][0]=0;
929 isKaonPionTPC[daught][1]=0;
08ea5c32 930 AliAODPid *pidObj = aodtrack->GetDetPid();
931 Double_t ptProng=pidObj->GetTPCmomentum();
932 //if(ptrack<0.6){
933 if(ptProng<0.6){
937e1290 934 if(TMath::Abs(nsigmaTPCK)<3.)isKaonPionTPC[daught][0]=1;
935 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
936 if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
937 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
938 }
08ea5c32 939 //else if(ptrack<.8){
940 else if(ptProng<.8){
937e1290 941 if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
942 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
943 if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
944 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
945 }
946 else {
947 if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
948 if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
949 }
950 }
951
952 }// END OF LOOP ON DAUGHTERS
953
954 // FURTHER PID REQUEST (both daughter info is needed)
955 if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
956 fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
957 return 0;
958 }
959 else if(hasPID[0]==0&&hasPID[1]==0){
960 fWhyRejection=28;// reject cases in which no PID info is available
961 return 0;
962 }
963 if(d->Pt()<2.){
964 // request of K identification at low D0 pt
965 combinedPID[0][0]=0;
966 combinedPID[0][1]=0;
967 combinedPID[1][0]=0;
968 combinedPID[1][1]=0;
969
970 combinedPID[0][0]=isKaonPionTOF[0][0]+isKaonPionTPC[0][0];
971 combinedPID[0][1]=isKaonPionTOF[0][1]+isKaonPionTPC[0][1];
aaf00f06 972 combinedPID[1][0]=isKaonPionTOF[1][0]+isKaonPionTPC[1][0];
973 combinedPID[1][1]=isKaonPionTOF[1][1]+isKaonPionTPC[1][1];
937e1290 974
975 if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
976 fWhyRejection=32;// reject cases where the Kaon is not identified
977 return 0;
978 }
979 }
980
981 // cout<<"Why? "<<fWhyRejection<<endl;
982 return isD0D0barPID[0]+isD0D0barPID[1];
983}
984
08ea5c32 985
986
937e1290 987//---------------------------------------------------------------------------
988Int_t AliRDHFCutsD0toKpi::CombineSelectionLevels(Int_t selectionvalTrack,
989 Int_t selectionvalCand,
990 Int_t selectionvalPID) const
991{
992 //
993 // This method combines the tracks, PID and cuts selection results
994 //
995 if(selectionvalTrack==0) return 0;
996
997 Int_t returnvalue;
c96634a2 998
937e1290 999 switch(selectionvalPID) {
1000 case 0:
1001 returnvalue=0;
1002 break;
1003 case 1:
1004 returnvalue=((selectionvalCand==1 || selectionvalCand==3) ? 1 : 0);
1005 break;
1006 case 2:
1007 returnvalue=((selectionvalCand==2 || selectionvalCand==3) ? 2 : 0);
1008 break;
1009 case 3:
1010 returnvalue=selectionvalCand;
1011 break;
1012 default:
1013 returnvalue=0;
1014 break;
1015 }
1016
1017 return returnvalue;
1018}
1019//----------------------------------------------------------------------------
1020Int_t AliRDHFCutsD0toKpi::IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
1021{
1022 //
1023 // Note: this method is temporary
1024 // Additional cuts on decay lenght and lower cut for d0 norm are applied using vertex without candidate's daughters
1025 //
c96634a2 1026
937e1290 1027 //apply cuts
aa8d25fe 1028
937e1290 1029 Float_t normDecLengthCut=1.,decLengthCut=TMath::Min(d->P()*0.0066+0.01,0.06/*cm*/), normd0Cut=0.5;
1030 // "decay length" expo law with tau' = beta*gamma*ctau= p/m*ctau =p*0.0123/1.864~p*0.0066
1031 // decay lenght > ctau' implies to retain (1-1/e) (for signal without considering detector resolution),
c96634a2 1032
937e1290 1033 Int_t returnvalue=3; //cut passed
1034 for(Int_t i=0;i<2/*prongs*/;i++){
1035 if(TMath::Abs(d->Normalizedd0Prong(i))<normd0Cut) return 0; //normd0Cut not passed
1036 }
be9cad0c 1037 if(d->DecayLength2()<decLengthCut*decLengthCut) return 0; //decLengthCut not passed
1038 if(d->NormalizedDecayLength2()<normDecLengthCut*normDecLengthCut) return 0; //decLengthCut not passed
937e1290 1039
c96634a2 1040
1041 return returnvalue;
1042}
08ea5c32 1043
5f25117b 1044//----------------------------------------------
5f25117b 1045void AliRDHFCutsD0toKpi::SetUseKF(Bool_t useKF)
1046{
1047 //switch on candidate selection via AliKFparticle
1048 if(!useKF) return;
1049 if(useKF){
1050 fUseKF=useKF;
1051 Int_t nvarsKF=11;
1052 SetNVars(nvarsKF);
1053 TString varNamesKF[11]={"inv. mass [GeV]",
1054 "dca [cm]",
1055 "cosThetaStar",
1056 "pTK [GeV/c]",
1057 "pTPi [GeV/c]",
1058 "d0K [cm]",
1059 "d0Pi [cm]",
1060 "d0d0 [cm^2]",
1061 "cosThetaPoint"
1062 "DecayLength[cm]",
1063 "RedChi2"};
1064 Bool_t isUpperCutKF[11]={kTRUE,
1065 kTRUE,
1066 kTRUE,
1067 kFALSE,
1068 kFALSE,
1069 kTRUE,
1070 kTRUE,
1071 kTRUE,
1072 kFALSE,
1073 kFALSE,
1074 kTRUE};
1075 SetVarNames(nvarsKF,varNamesKF,isUpperCutKF);
1076 SetGlobalIndex();
1077 Bool_t forOpt[11]={kFALSE,
1078 kTRUE,
1079 kTRUE,
1080 kFALSE,
1081 kFALSE,
1082 kFALSE,
1083 kFALSE,
1084 kTRUE,
1085 kTRUE,
1086 kFALSE,
1087 kFALSE};
1088 SetVarsForOpt(4,forOpt);
1089 }
1090 return;
1091}
b14fe7f1 1092
1093
73f4ccee 1094void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {
1095 //
1096 //STANDARD CUTS USED FOR 2010 pp analysis
1097 //dca cut will be enlarged soon to 400 micron
1098 //
b14fe7f1 1099
1100 SetName("D0toKpiCutsStandard");
1101 SetTitle("Standard Cuts for D0 analysis");
1102
1103 // PILE UP REJECTION
1104 SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1105
1106 // EVENT CUTS
1107 SetMinVtxContr(1);
1108
1109
1110 // TRACKS ON SINGLE TRACKS
1111 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1112 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1113 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1114 esdTrackCuts->SetRequireITSRefit(kTRUE);
1115 // esdTrackCuts->SetMinNClustersITS(4);
1116 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1117 esdTrackCuts->SetMinDCAToVertexXY(0.);
1118 esdTrackCuts->SetEtaRange(-0.8,0.8);
1119 esdTrackCuts->SetPtRange(0.3,1.e10);
1120
1121 AddTrackCuts(esdTrackCuts);
1122
1123 const Int_t nptbins =13;
1124 const Double_t ptmax = 9999.;
1125 const Int_t nvars=9;
1126 Float_t ptbins[nptbins+1];
1127 ptbins[0]=0.;
1128 ptbins[1]=0.5;
1129 ptbins[2]=1.;
1130 ptbins[3]=2.;
1131 ptbins[4]=3.;
1132 ptbins[5]=4.;
1133 ptbins[6]=5.;
1134 ptbins[7]=6.;
1135 ptbins[8]=8.;
1136 ptbins[9]=12.;
1137 ptbins[10]=16.;
1138 ptbins[11]=20.;
1139 ptbins[12]=24.;
1140 ptbins[13]=ptmax;
1141
1142 SetGlobalIndex(nvars,nptbins);
1143 SetPtBins(nptbins+1,ptbins);
1144
1145 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*/
1146 {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*/
1147 {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 */
1148 {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 */
1149 {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 */
1150 {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 */
1151 {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 */
1152 {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 */
1153 {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 */
1154 {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 */
1155 {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 */
1156 {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 */
1157 {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
1159
1160 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1161 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1162 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1163
1164 for (Int_t ibin=0;ibin<nptbins;ibin++){
1165 for (Int_t ivar = 0; ivar<nvars; ivar++){
1166 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1167 }
1168 }
1169
1170 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1171 SetUseSpecialCuts(kTRUE);
1172 SetRemoveDaughtersFromPrim(kTRUE);
1173
e11ae259 1174 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1175 delete [] cutsMatrixTransposeStand;
1176 cutsMatrixTransposeStand=NULL;
1177
b14fe7f1 1178 // PID SETTINGS
1179 AliAODPidHF* pidObj=new AliAODPidHF();
1180 //pidObj->SetName("pid4D0");
1181 Int_t mode=1;
1182 const Int_t nlims=2;
1183 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1184 Bool_t compat=kTRUE; //effective only for this mode
1185 Bool_t asym=kTRUE;
1186 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1187 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1188 pidObj->SetMatch(mode);
1189 pidObj->SetPLimit(plims,nlims);
1190 pidObj->SetSigma(sigmas);
1191 pidObj->SetCompat(compat);
1192 pidObj->SetTPC(kTRUE);
1193 pidObj->SetTOF(kTRUE);
1194
1195 SetPidHF(pidObj);
1196 SetUsePID(kTRUE);
1197 SetUseDefaultPID(kFALSE);
1198
1199
1200 PrintAll();
1201
e11ae259 1202 delete pidObj;
1203 pidObj=NULL;
1204
b14fe7f1 1205 return;
1206
1207}
1208
1209
73f4ccee 1210void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
1211 //
1212 //PRELIMINARY CUTS USED FOR 2010 PbPb analysis
1213 //... EVOLVING SOON
1214 //
b14fe7f1 1215
1216 SetName("D0toKpiCutsStandard");
1217 SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
1218
1219 // PILE UP REJECTION
1220 //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1221
1222 // EVENT CUTS
1223 SetMinVtxContr(1);
1224
1225
1226 // TRACKS ON SINGLE TRACKS
1227 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1228 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1229 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1230 esdTrackCuts->SetRequireITSRefit(kTRUE);
1231 // esdTrackCuts->SetMinNClustersITS(4);
1232 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1233 esdTrackCuts->SetMinDCAToVertexXY(0.);
1234 esdTrackCuts->SetEtaRange(-0.8,0.8);
1235 esdTrackCuts->SetPtRange(0.3,1.e10);
1236
1237 AddTrackCuts(esdTrackCuts);
1238
1239 const Int_t nptbins =13;
1240 const Double_t ptmax = 9999.;
1241 const Int_t nvars=9;
1242 Float_t ptbins[nptbins+1];
1243 ptbins[0]=0.;
1244 ptbins[1]=0.5;
1245 ptbins[2]=1.;
1246 ptbins[3]=2.;
1247 ptbins[4]=3.;
1248 ptbins[5]=4.;
1249 ptbins[6]=5.;
1250 ptbins[7]=6.;
1251 ptbins[8]=8.;
1252 ptbins[9]=12.;
1253 ptbins[10]=16.;
1254 ptbins[11]=20.;
1255 ptbins[12]=24.;
1256 ptbins[13]=ptmax;
1257
1258 SetGlobalIndex(nvars,nptbins);
1259 SetPtBins(nptbins+1,ptbins);
1260
1261 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*/
1262 {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*/
1263 {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 */
1264 {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 */
1265 {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 */
1266 {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 */
1267 {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 */
1268 {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 */
1269 {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 */
1270 {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 */
1271 {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 */
1272 {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 */
1273 {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
1275
1276 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1277 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1278 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1279
1280 for (Int_t ibin=0;ibin<nptbins;ibin++){
1281 for (Int_t ivar = 0; ivar<nvars; ivar++){
1282 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1283 }
1284 }
1285
1286 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1287 SetUseSpecialCuts(kTRUE);
1288 SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
e11ae259 1289 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1290 delete [] cutsMatrixTransposeStand;
1291 cutsMatrixTransposeStand=NULL;
b14fe7f1 1292
1293 // PID SETTINGS
1294 AliAODPidHF* pidObj=new AliAODPidHF();
1295 //pidObj->SetName("pid4D0");
1296 Int_t mode=1;
1297 const Int_t nlims=2;
1298 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1299 Bool_t compat=kTRUE; //effective only for this mode
1300 Bool_t asym=kTRUE;
1301 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1302 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1303 pidObj->SetMatch(mode);
1304 pidObj->SetPLimit(plims,nlims);
1305 pidObj->SetSigma(sigmas);
1306 pidObj->SetCompat(compat);
1307 pidObj->SetTPC(kTRUE);
1308 pidObj->SetTOF(kTRUE);
1309
1310 SetPidHF(pidObj);
1311 SetUsePID(kTRUE);
1312 SetUseDefaultPID(kTRUE);// TEMPORARY: PROTON EXCLUSION SET ONLY IN DEFAULT PID
1313
1314
1315 PrintAll();
1316
e11ae259 1317
1318 delete pidObj;
1319 pidObj=NULL;
1320
b14fe7f1 1321 return;
1322
1323}