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