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