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