]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliRDHFCutsD0toKpi.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliRDHFCutsD0toKpi.cxx
CommitLineData
fc051756 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
16/* $Id$ */
17
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"
32#include "AliAODPid.h"
33#include "AliTPCPIDResponse.h"
34#include "AliAODVertex.h"
35#include "AliKFVertex.h"
36#include "AliKFParticle.h"
37
c64cb1f6 38using std::cout;
39using std::endl;
40
fc051756 41ClassImp(AliRDHFCutsD0toKpi)
42
43//--------------------------------------------------------------------------
b8f433ab 44AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const char* name) :
45 AliRDHFCuts(name),
46 fUseSpecialCuts(kFALSE),
47 fLowPt(kTRUE),
48 fDefaultPID(kFALSE),
49 fUseKF(kFALSE),
50 fPtLowPID(2.),
51 fPtMaxSpecialCuts(9999.),
52 fmaxPtrackForPID(9999),
53 fCombPID(kFALSE),
54 fnSpecies(AliPID::kSPECIES),
55 fWeightsPositive(0),
56 fWeightsNegative(0),
57 fProbThreshold(0.5),
58 fBayesianStrategy(0),
59 fBayesianCondition(0)
fc051756 60{
61 //
62 // Default Constructor
63 //
64 Int_t nvars=11;
65 SetNVars(nvars);
66 TString varNames[11]={"inv. mass [GeV]",
67 "dca [cm]",
68 "cosThetaStar",
69 "pTK [GeV/c]",
70 "pTPi [GeV/c]",
71 "d0K [cm]",
72 "d0Pi [cm]",
73 "d0d0 [cm^2]",
74 "cosThetaPoint",
75 "|cosThetaPointXY|",
76 "NormDecayLenghtXY"};
77 Bool_t isUpperCut[11]={kTRUE,
78 kTRUE,
79 kTRUE,
80 kFALSE,
81 kFALSE,
82 kTRUE,
83 kTRUE,
84 kTRUE,
85 kFALSE,
86 kFALSE,
87 kFALSE};
88 SetVarNames(nvars,varNames,isUpperCut);
89 Bool_t forOpt[11]={kFALSE,
90 kTRUE,
91 kTRUE,
92 kFALSE,
93 kFALSE,
94 kFALSE,
95 kFALSE,
96 kTRUE,
97 kTRUE,
98 kFALSE,
99 kFALSE};
100 SetVarsForOpt(4,forOpt);
101 Float_t limits[2]={0,999999999.};
102 SetPtBins(2,limits);
103
c2ec6218 104 fWeightsNegative = new Double_t[AliPID::kSPECIES];
105 fWeightsPositive = new Double_t[AliPID::kSPECIES];
106
17454b62 107 for (Int_t i = 0; i<AliPID::kSPECIES; i++) {
108 fWeightsPositive[i] = 0.;
109 fWeightsNegative[i] = 0.;
110 }
111
fc051756 112}
113//--------------------------------------------------------------------------
114AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :
115 AliRDHFCuts(source),
116 fUseSpecialCuts(source.fUseSpecialCuts),
117 fLowPt(source.fLowPt),
118 fDefaultPID(source.fDefaultPID),
119 fUseKF(source.fUseKF),
120 fPtLowPID(source.fPtLowPID),
1d9bf4ef 121 fPtMaxSpecialCuts(source.fPtMaxSpecialCuts),
e2aa82b6 122 fmaxPtrackForPID(source.fmaxPtrackForPID),
123 fCombPID(source.fCombPID),
c2ec6218 124 fnSpecies(source.fnSpecies),
e2aa82b6 125 fWeightsPositive(source.fWeightsPositive),
126 fWeightsNegative(source.fWeightsNegative),
b8f433ab 127 fProbThreshold(source.fProbThreshold),
128 fBayesianStrategy(source.fBayesianStrategy),
129 fBayesianCondition(source.fBayesianCondition)
fc051756 130{
131 //
132 // Copy constructor
133 //
134
135}
136//--------------------------------------------------------------------------
137AliRDHFCutsD0toKpi &AliRDHFCutsD0toKpi::operator=(const AliRDHFCutsD0toKpi &source)
138{
139 //
140 // assignment operator
141 //
142 if(&source == this) return *this;
143
144 AliRDHFCuts::operator=(source);
145 fUseSpecialCuts=source.fUseSpecialCuts;
146 fLowPt=source.fLowPt;
147 fDefaultPID=source.fDefaultPID;
148 fUseKF=source.fUseKF;
149 fPtLowPID=source.fPtLowPID;
150 fPtMaxSpecialCuts=source.fPtMaxSpecialCuts;
1d9bf4ef 151 fmaxPtrackForPID=source.fmaxPtrackForPID;
e2aa82b6 152 fCombPID = source.fCombPID;
c2ec6218 153 fnSpecies = source.fnSpecies;
e2aa82b6 154 fWeightsPositive = source.fWeightsPositive;
155 fWeightsNegative = source.fWeightsNegative;
b8f433ab 156 fProbThreshold = source.fProbThreshold;
e2aa82b6 157 fBayesianStrategy = source.fBayesianStrategy;
b8f433ab 158 fBayesianCondition = source.fBayesianCondition;
fc051756 159 return *this;
160}
161
162
e2aa82b6 163
164//---------------------------------------------------------------------------
165AliRDHFCutsD0toKpi::~AliRDHFCutsD0toKpi()
166{
167//
168// Destructor
169//
170 if (fWeightsNegative) {
171 delete[] fWeightsNegative;
172 fWeightsNegative = 0;
173 }
174 if (fWeightsPositive) {
175 delete[] fWeightsNegative;
176 fWeightsNegative = 0;
177 }
178}
179
fc051756 180//---------------------------------------------------------------------------
181void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent *aod) {
182 //
183 // Fills in vars the values of the variables
184 //
185
186 if(nvars!=fnVarsForOpt) {
187 printf("AliRDHFCutsD0toKpi::GetCutsVarsForOpt: wrong number of variables\n");
188 return;
189 }
190
191 AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)d;
192
193 //recalculate vertex w/o daughters
194 Bool_t cleanvtx=kFALSE;
195 AliAODVertex *origownvtx=0x0;
196 if(fRemoveDaughtersFromPrimary) {
197 if(dd->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*dd->GetOwnPrimaryVtx());
198 cleanvtx=kTRUE;
199 if(!RecalcOwnPrimaryVtx(dd,aod)) {
200 CleanOwnPrimaryVtx(dd,aod,origownvtx);
201 cleanvtx=kFALSE;
202 }
203 }
204
205 Int_t iter=-1;
206 if(fVarsForOpt[0]){
207 iter++;
208 if(TMath::Abs(pdgdaughters[0])==211) {
209 vars[iter]=dd->InvMassD0();
210 } else {
211 vars[iter]=dd->InvMassD0bar();
212 }
213 }
214 if(fVarsForOpt[1]){
215 iter++;
216 vars[iter]=dd->GetDCA();
217 }
218 if(fVarsForOpt[2]){
219 iter++;
220 if(TMath::Abs(pdgdaughters[0])==211) {
221 vars[iter] = dd->CosThetaStarD0();
222 } else {
223 vars[iter] = dd->CosThetaStarD0bar();
224 }
225 }
226 if(fVarsForOpt[3]){
227 iter++;
228 if(TMath::Abs(pdgdaughters[0])==321) {
229 vars[iter]=dd->PtProng(0);
230 }
231 else{
232 vars[iter]=dd->PtProng(1);
233 }
234 }
235 if(fVarsForOpt[4]){
236 iter++;
237 if(TMath::Abs(pdgdaughters[0])==211) {
238 vars[iter]=dd->PtProng(0);
239 }
240 else{
241 vars[iter]=dd->PtProng(1);
242 }
243 }
244 if(fVarsForOpt[5]){
245 iter++;
246 if(TMath::Abs(pdgdaughters[0])==321) {
247 vars[iter]=dd->Getd0Prong(0);
248 }
249 else{
250 vars[iter]=dd->Getd0Prong(1);
251 }
252 }
253 if(fVarsForOpt[6]){
254 iter++;
255 if(TMath::Abs(pdgdaughters[0])==211) {
256 vars[iter]=dd->Getd0Prong(0);
257 }
258 else{
259 vars[iter]=dd->Getd0Prong(1);
260 }
261 }
262 if(fVarsForOpt[7]){
263 iter++;
264 vars[iter]= dd->Prodd0d0();
265 }
266 if(fVarsForOpt[8]){
267 iter++;
268 vars[iter]=dd->CosPointingAngle();
269 }
270
271 if(fVarsForOpt[9]){
272 iter++;
273 vars[iter]=TMath::Abs(dd->CosPointingAngleXY());
274 }
275
276 if(fVarsForOpt[10]){
277 iter++;
278 vars[iter]=dd->NormalizedDecayLengthXY();
279 }
280
281 if(cleanvtx)CleanOwnPrimaryVtx(dd,aod,origownvtx);
282
283 return;
284}
285//---------------------------------------------------------------------------
286Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod) {
287 //
288 // Apply selection
289 //
290
291
292 fIsSelectedCuts=0;
293 fIsSelectedPID=0;
294
295 if(!fCutsRD){
296 cout<<"Cut matrice not inizialized. Exit..."<<endl;
297 return 0;
298 }
299 //PrintAll();
300 AliAODRecoDecayHF2Prong* d=(AliAODRecoDecayHF2Prong*)obj;
fc051756 301 if(!d){
302 cout<<"AliAODRecoDecayHF2Prong null"<<endl;
303 return 0;
304 }
305
306 if(fKeepSignalMC) if(IsSignalMC(d,aod,421)) return 3;
307
308 Double_t ptD=d->Pt();
309 if(ptD<fMinPtCand) return 0;
310 if(ptD>fMaxPtCand) return 0;
311
3e075b37 312 if(fUseTrackSelectionWithFilterBits && d->HasBadDaughters()) return 0;
fc051756 313
314 // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
315 Int_t returnvaluePID=3;
316 Int_t returnvalueCuts=3;
317
318 // selection on candidate
319 if(selectionLevel==AliRDHFCuts::kAll ||
320 selectionLevel==AliRDHFCuts::kCandidate) {
321
322 if(!fUseKF) {
323
324 //recalculate vertex w/o daughters
325 AliAODVertex *origownvtx=0x0;
326 if(fRemoveDaughtersFromPrimary && !fUseMCVertex) {
327 if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
328 if(!RecalcOwnPrimaryVtx(d,aod)) {
329 CleanOwnPrimaryVtx(d,aod,origownvtx);
330 return 0;
331 }
332 }
333
334 if(fUseMCVertex) {
335 if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
336 if(!SetMCPrimaryVtx(d,aod)) {
337 CleanOwnPrimaryVtx(d,aod,origownvtx);
338 return 0;
339 }
340 }
341
342 Double_t pt=d->Pt();
343
344 Int_t okD0=0,okD0bar=0;
345
346 Int_t ptbin=PtBin(pt);
347 if (ptbin==-1) {
348 CleanOwnPrimaryVtx(d,aod,origownvtx);
349 return 0;
350 }
351
352 Double_t mD0,mD0bar,ctsD0,ctsD0bar;
353 okD0=1; okD0bar=1;
354
355 Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
356
357 d->InvMassD0(mD0,mD0bar);
358 if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
359 if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
360 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
361
362 if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
363
364
365 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;
366 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;
367 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
368
369
370 if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
371 TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
372 if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
373 TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
374 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
375
376 if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
377
378
379 d->CosThetaStarD0(ctsD0,ctsD0bar);
380 if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
381 if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
382 if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
383
384 if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
385
386
387 if(TMath::Abs(d->CosPointingAngleXY()) < fCutsRD[GetGlobalIndex(9,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
388
389 Double_t normalDecayLengXY=d->NormalizedDecayLengthXY();
390 if (normalDecayLengXY < fCutsRD[GetGlobalIndex(10, ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
391
392 if (returnvalueCuts!=0) {
393 if (okD0) returnvalueCuts=1; //cuts passed as D0
394 if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar
395 if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
396 }
397
398 // call special cuts
399 Int_t special=1;
400 if(fUseSpecialCuts && (pt<fPtMaxSpecialCuts)) special=IsSelectedSpecialCuts(d);
401 if(!special) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
402
403 // unset recalculated primary vertex when not needed any more
404 CleanOwnPrimaryVtx(d,aod,origownvtx);
405
406 } else {
407 // go to selection with Kalman vertexing, if requested
408 returnvalueCuts = IsSelectedKF(d,aod);
409 }
410
411 fIsSelectedCuts=returnvalueCuts;
412 if(!returnvalueCuts) return 0;
413 }
414
415
416 // selection on PID
417 if(selectionLevel==AliRDHFCuts::kAll ||
418 selectionLevel==AliRDHFCuts::kCandidate ||
419 selectionLevel==AliRDHFCuts::kPID) {
e2aa82b6 420 if (!fCombPID) {
421 returnvaluePID = IsSelectedPID(d);
422 fIsSelectedPID=returnvaluePID;
423 if(!returnvaluePID) return 0;
424 } else {
775d1f1c 425 returnvaluePID = IsSelectedCombPID(d);
426 if(!returnvaluePID) return 0;
e2aa82b6 427 }
428 }
775d1f1c 429
fc051756 430
e2aa82b6 431
fc051756 432 Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
433
434 if(!returnvalueComb) return 0;
435
436 // selection on daughter tracks
437 if(selectionLevel==AliRDHFCuts::kAll ||
438 selectionLevel==AliRDHFCuts::kTracks) {
439 if(!AreDaughtersSelected(d)) return 0;
440 }
441
442 // cout<<"Pid = "<<returnvaluePID<<endl;
443 return returnvalueComb;
444}
445
446//------------------------------------------------------------------------------------------
447Int_t AliRDHFCutsD0toKpi::IsSelectedKF(AliAODRecoDecayHF2Prong *d,
448 AliAODEvent* aod) const {
449 //
450 // Apply selection using KF-vertexing
451 //
452
453 AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
454 AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
455
456 if(!track0 || !track1) {
457 cout<<"one or two D0 daughters missing!"<<endl;
458 return 0;
459 }
460
461 // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
462 Int_t returnvalueCuts=3;
463
464 // candidate selection with AliKF
465 AliKFParticle::SetField(aod->GetMagneticField()); // set the magnetic field
466
467 Int_t okD0=0,okD0bar=0;
468 okD0=1; okD0bar=1;
469
470 // convert tracks into AliKFParticles
471
472 AliKFParticle negPiKF(*track1,-211); // neg pion kandidate
473 AliKFParticle negKKF(*track1,-321); // neg kaon kandidate
474 AliKFParticle posPiKF(*track0,211); // pos pion kandidate
475 AliKFParticle posKKF(*track0,321); // pos kaon kandidate
476
477 // build D0 candidates
478
479 AliKFParticle d0c(negKKF,posPiKF); // D0 candidate
480 AliKFParticle ad0c(posKKF,negPiKF); // D0bar candidate
481
482 // create kf primary vertices
483
484 AliAODVertex *vtx1 = aod->GetPrimaryVertex();
485 AliKFVertex primVtx(*vtx1);
486 AliKFVertex aprimVtx(*vtx1);
487
488 if(primVtx.GetNContributors()<=0) okD0 = 0;
489 if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
490 if(!okD0 && !okD0bar) returnvalueCuts=0;
491
492 // calculate mass
493
494 Double_t d0mass = d0c.GetMass();
495 Double_t ad0mass = ad0c.GetMass();
496
497 // calculate P of D0 and D0bar
498 Double_t d0P = d0c.GetP();
499 Double_t d0Px = d0c.GetPx();
500 Double_t d0Py = d0c.GetPy();
501 Double_t d0Pz = d0c.GetPz();
502 Double_t ad0P = ad0c.GetP();
503 Double_t ad0Px = ad0c.GetPx();
504 Double_t ad0Py = ad0c.GetPy();
505 Double_t ad0Pz = ad0c.GetPz();
506
507 //calculate Pt of D0 and D0bar
508
509 Double_t pt=d0c.GetPt();
510 Double_t apt=ad0c.GetPt();
511
512 // remove D0 daughters from primary vertices (if used in vertex fit) and add D0-candidates
513
514 if(track0->GetUsedForPrimVtxFit()) {
515 primVtx -= posPiKF;
516 aprimVtx -= posKKF;
517 }
518
519 if(track1->GetUsedForPrimVtxFit()) {
520 primVtx -= negKKF;
521 aprimVtx -= negPiKF;
522 }
523
524 primVtx += d0c;
525 aprimVtx += ad0c;
526
527 if(primVtx.GetNContributors()<=0) okD0 = 0;
528 if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
529 if(!okD0 && !okD0bar) returnvalueCuts=0;
530
531 //calculate cut variables
532
533 // calculate impact params of daughters w.r.t recalculated vertices
534
535 Double_t impactPi = posPiKF.GetDistanceFromVertexXY(primVtx);
536 Double_t aimpactPi = negPiKF.GetDistanceFromVertexXY(aprimVtx);
537 Double_t impactKa = negKKF.GetDistanceFromVertexXY(primVtx);
538 Double_t aimpactKa = posKKF.GetDistanceFromVertexXY(aprimVtx);
539
540 // calculate Product of Impact Params
541
542 Double_t prodParam = impactPi*impactKa;
543 Double_t aprodParam = aimpactPi*aimpactKa;
544
545 // calculate cosine of pointing angles
546
547 TVector3 mom(d0c.GetPx(),d0c.GetPy(),d0c.GetPz());
548 TVector3 fline(d0c.GetX()-primVtx.GetX(),
549 d0c.GetY()-primVtx.GetY(),
550 d0c.GetZ()-primVtx.GetZ());
551 Double_t pta = mom.Angle(fline);
552 Double_t cosP = TMath::Cos(pta); // cosine of pta for D0 candidate
553
554 TVector3 amom(ad0c.GetPx(),ad0c.GetPy(),ad0c.GetPz());
555 TVector3 afline(ad0c.GetX()-aprimVtx.GetX(),
556 ad0c.GetY()-aprimVtx.GetY(),
557 ad0c.GetZ()-aprimVtx.GetZ());
558 Double_t apta = amom.Angle(afline);
559 Double_t acosP = TMath::Cos(apta); // cosine of pta for D0bar candidate
560
561 // calculate P of Pions at Decay Position of D0 and D0bar candidates
562 negKKF.TransportToParticle(d0c);
563 posPiKF.TransportToParticle(d0c);
564 posKKF.TransportToParticle(ad0c);
565 negPiKF.TransportToParticle(ad0c);
566
567 Double_t pxPi = posPiKF.GetPx();
568 Double_t pyPi = posPiKF.GetPy();
569 Double_t pzPi = posPiKF.GetPz();
570 Double_t ptPi = posPiKF.GetPt();
571
572 Double_t apxPi = negPiKF.GetPx();
573 Double_t apyPi = negPiKF.GetPy();
574 Double_t apzPi = negPiKF.GetPz();
575 Double_t aptPi = negPiKF.GetPt();
576
577 // calculate Pt of Kaons at Decay Position of D0 and D0bar candidates
578
579 Double_t ptK = negKKF.GetPt();
580 Double_t aptK = posKKF.GetPt();
581
582 //calculate cos(thetastar)
583 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(421)->Mass();
584 Double_t massp[2];
585 massp[0] = TDatabasePDG::Instance()->GetParticle(321)->Mass();
586 massp[1] = TDatabasePDG::Instance()->GetParticle(211)->Mass();
587 Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)
588 -4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
589
590 // cos(thetastar) for D0 and Pion
591
592 Double_t d0E = TMath::Sqrt(massvtx*massvtx + d0P*d0P);
593 Double_t beta = d0P/d0E;
594 Double_t gamma = d0E/massvtx;
595 TVector3 momPi(pxPi,pyPi,pzPi);
596 TVector3 momTot(d0Px,d0Py,d0Pz);
597 Double_t q1 = momPi.Dot(momTot)/momTot.Mag();
598 Double_t cts = (q1/gamma-beta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
599
600 // cos(thetastar) for D0bar and Pion
601
602 Double_t ad0E = TMath::Sqrt(massvtx*massvtx + ad0P*ad0P);
603 Double_t abeta = ad0P/ad0E;
604 Double_t agamma = ad0E/massvtx;
605 TVector3 amomPi(apxPi,apyPi,apzPi);
606 TVector3 amomTot(ad0Px,ad0Py,ad0Pz);
607 Double_t aq1 = amomPi.Dot(amomTot)/amomTot.Mag();
608 Double_t acts = (aq1/agamma-abeta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
609
610 // calculate reduced Chi2 for the full D0 fit
611 d0c.SetProductionVertex(primVtx);
612 ad0c.SetProductionVertex(aprimVtx);
613 negKKF.SetProductionVertex(d0c);
614 posPiKF.SetProductionVertex(d0c);
615 posKKF.SetProductionVertex(ad0c);
616 negPiKF.SetProductionVertex(ad0c);
617 d0c.TransportToProductionVertex();
618 ad0c.TransportToProductionVertex();
619
620 // calculate the decay length
621 Double_t decayLengthD0 = d0c.GetDecayLength();
622 Double_t adecayLengthD0 = ad0c.GetDecayLength();
623
624 Double_t chi2D0 = 50.;
625 if(d0c.GetNDF() > 0 && d0c.GetChi2() >= 0) {
626 chi2D0 = d0c.GetChi2()/d0c.GetNDF();
627 }
628
629 Double_t achi2D0 = 50.;
630 if(ad0c.GetNDF() > 0 && ad0c.GetChi2() >= 0) {
631 achi2D0 = ad0c.GetChi2()/ad0c.GetNDF();
632 }
633
634 // Get the Pt-bins
635 Int_t ptbin=PtBin(pt);
636 Int_t aptbin=PtBin(apt);
637
638 if(ptbin < 0) okD0 = 0;
639 if(aptbin < 0) okD0bar = 0;
640 if(!okD0 && !okD0bar) returnvalueCuts=0;
641
642 if(ptK < fCutsRD[GetGlobalIndex(3,ptbin)] || ptPi < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
643 if(aptK < fCutsRD[GetGlobalIndex(3,aptbin)] || aptPi < fCutsRD[GetGlobalIndex(4,aptbin)]) okD0bar = 0;
644 if(!okD0 && !okD0bar) returnvalueCuts=0;
645
646
647 if(TMath::Abs(impactKa) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
648 TMath::Abs(impactPi) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
649
650 if(TMath::Abs(aimpactKa) > fCutsRD[GetGlobalIndex(5,aptbin)] ||
651 TMath::Abs(aimpactPi) > fCutsRD[GetGlobalIndex(6,aptbin)]) okD0bar = 0;
652
653 if(!okD0 && !okD0bar) returnvalueCuts=0;
654
655 // for the moment via the standard method due to bug in AliKF
656 if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) okD0 = 0;
657 if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,aptbin)]) okD0bar = 0;
658 if(!okD0 && !okD0bar) returnvalueCuts=0;
659
660
661 if(TMath::Abs(d0mass-massvtx) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
662 if(TMath::Abs(ad0mass-massvtx) > fCutsRD[GetGlobalIndex(0,aptbin)]) okD0bar = 0;
663 if(!okD0 && !okD0bar) returnvalueCuts=0;
664
665
666 if(TMath::Abs(cts) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
667 if(TMath::Abs(acts) > fCutsRD[GetGlobalIndex(2,aptbin)]) okD0bar = 0;
668 if(!okD0 && !okD0bar) returnvalueCuts=0;
669
670 if(prodParam > fCutsRD[GetGlobalIndex(7,ptbin)]) okD0 = 0;
671 if(aprodParam > fCutsRD[GetGlobalIndex(7,aptbin)]) okD0bar = 0;
672 if(!okD0 && !okD0bar) returnvalueCuts=0;
673
674 if(cosP < fCutsRD[GetGlobalIndex(8,ptbin)]) okD0 = 0;
675 if(acosP < fCutsRD[GetGlobalIndex(8,aptbin)]) okD0bar = 0;
676 if(!okD0 && !okD0bar) returnvalueCuts=0;
677
678 if(chi2D0 > fCutsRD[GetGlobalIndex(10,ptbin)]) okD0 = 0;
679 if(achi2D0 > fCutsRD[GetGlobalIndex(10,aptbin)]) okD0bar = 0;
680 if(!okD0 && !okD0bar) returnvalueCuts=0;
681
682 if(decayLengthD0 < fCutsRD[GetGlobalIndex(9,ptbin)]) okD0 = 0;
683 if(adecayLengthD0 < fCutsRD[GetGlobalIndex(9,aptbin)]) okD0bar = 0;
684 if(!okD0 && !okD0bar) returnvalueCuts=0;
685
686 if(returnvalueCuts!=0) {
687 if(okD0) returnvalueCuts=1; //cuts passed as D0
688 if(okD0bar) returnvalueCuts=2; //cuts passed as D0bar
689 if(okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
690 }
691
692 return returnvalueCuts;
693}
694
695//---------------------------------------------------------------------------
696
697Bool_t AliRDHFCutsD0toKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
698{
699 //
700 // Checking if D0 is in fiducial acceptance region
701 //
702
af5d687e 703 if(fMaxRapidityCand>-998.){
704 if(TMath::Abs(y) > fMaxRapidityCand) return kFALSE;
705 else return kTRUE;
706 }
707
fc051756 708 if(pt > 5.) {
709 // applying cut for pt > 5 GeV
710 AliDebug(2,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt));
711 if (TMath::Abs(y) > 0.8){
712 return kFALSE;
713 }
714 } else {
715 // appliying smooth cut for pt < 5 GeV
716 Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
717 Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
718 AliDebug(2,Form("pt of D0 = %f (< 5), cutting according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
719 if (y < minFiducialY || y > maxFiducialY){
720 return kFALSE;
721 }
722 }
723
724 return kTRUE;
725}
726//---------------------------------------------------------------------------
727Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d)
728{
729 // ############################################################
730 //
731 // Apply PID selection
732 //
733 //
734 // ############################################################
735
736 if(!fUsePID) return 3;
737 if(fDefaultPID) return IsSelectedPIDdefault(d);
e2aa82b6 738 if (fCombPID) return IsSelectedCombPID(d); //to use Bayesian method if applicable
fc051756 739 fWhyRejection=0;
740 Int_t isD0D0barPID[2]={1,2};
741 Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
742 // same for prong 2
743 // values convention -1 = discarded
744 // 0 = not identified (but compatible) || No PID (->hasPID flag)
745 // 1 = identified
746 // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
747 // Initial hypothesis: unknwon (but compatible)
748 combinedPID[0][0]=0; // prima figlia, Kaon
749 combinedPID[0][1]=0; // prima figlia, pione
750 combinedPID[1][0]=0; // seconda figlia, Kaon
751 combinedPID[1][1]=0; // seconda figlia, pion
752
753 Bool_t checkPIDInfo[2]={kTRUE,kTRUE};
754 Double_t sigma_tmp[3]={fPidHF->GetSigma(0),fPidHF->GetSigma(1),fPidHF->GetSigma(2)};
dbbe61a2 755 Bool_t isTOFused=fPidHF->GetTOF(),isCompat=fPidHF->GetCompat();
fc051756 756 for(Int_t daught=0;daught<2;daught++){
757 //Loop con prongs
758 AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
759 if(fPidHF->IsTOFPiKexcluded(aodtrack,5.)) return 0;
760
761 if(!(fPidHF->CheckStatus(aodtrack,"TPC")) && !(fPidHF->CheckStatus(aodtrack,"TOF"))) {
762 checkPIDInfo[daught]=kFALSE;
763 continue;
764 }
1d9bf4ef 765 Double_t pProng=aodtrack->P();
fc051756 766
767 // identify kaon
1d9bf4ef 768 if(pProng<fmaxPtrackForPID){
769 combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
770 }
fc051756 771 // identify pion
1d9bf4ef 772 if(pProng<fmaxPtrackForPID){
773 if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
774 combinedPID[daught][1]=0;
775 }else{
776 fPidHF->SetTOF(kFALSE);
777 combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
dbbe61a2 778 if(isTOFused)fPidHF->SetTOF(kTRUE);
779 if(isCompat)fPidHF->SetCompat(kTRUE);
1d9bf4ef 780 }
781 }
dbbe61a2 782
fc051756 783 if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded
784 isD0D0barPID[0]=0;
785 isD0D0barPID[1]=0;
786 }
787 else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){
788 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
789 else isD0D0barPID[0]=0;// if K+ D0 excluded
790 }
791 /* else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){
792 isD0D0barPID[0]=0;
793 isD0D0barPID[1]=0;
794 }
795 */
796 else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
797 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
798 else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
799 }
800 else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
801 if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
802 else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
803 }
804
805 if(fLowPt && d->Pt()<fPtLowPID){
806 Double_t sigmaTPC[3]={3.,2.,0.};
807 fPidHF->SetSigmaForTPC(sigmaTPC);
808 // identify kaon
809 combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
810
1d9bf4ef 811 if(pProng<0.6){
fc051756 812 fPidHF->SetCompat(kFALSE);
813 combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
814 fPidHF->SetCompat(kTRUE);
815 }
816
817 if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
818 combinedPID[daught][1]=0;
819 }else{
820 fPidHF->SetTOF(kFALSE);
821 Double_t sigmaTPCpi[3]={3.,3.,0.};
822 fPidHF->SetSigmaForTPC(sigmaTPCpi);
823 combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
824 fPidHF->SetTOF(kTRUE);
1d9bf4ef 825 if(pProng<0.8){
fc051756 826 Bool_t isTPCpion=fPidHF->IsPionRaw(aodtrack,"TPC");
827 if(isTPCpion){
828 combinedPID[daught][1]=1;
829 }else{
830 combinedPID[daught][1]=-1;
831 }
832 }
833 }
834
835 }
836 fPidHF->SetSigmaForTPC(sigma_tmp);
837 }// END OF LOOP ON DAUGHTERS
838
839 if(!checkPIDInfo[0] && !checkPIDInfo[1]) {
840 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
841 return 0;
842 }
843
844
845 // FURTHER PID REQUEST (both daughter info is needed)
846 if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
847 fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
848 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
849 return 0;
850 }
851
852 if(fLowPt && d->Pt()<fPtLowPID){
853 if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
854 fWhyRejection=32;// reject cases where the Kaon is not identified
855 fPidHF->SetSigmaForTPC(sigma_tmp);
856 return 0;
857 }
858 }
859 if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
860
861 // cout<<"Why? "<<fWhyRejection<<endl;
862 return isD0D0barPID[0]+isD0D0barPID[1];
863}
864//---------------------------------------------------------------------------
865Int_t AliRDHFCutsD0toKpi::IsSelectedPIDdefault(AliAODRecoDecayHF* d)
866{
867 // ############################################################
868 //
869 // Apply PID selection
870 //
871 //
872 // temporary selection: PID AS USED FOR D0 by Andrea Rossi (up to 28/06/2010)
873 //
874 // d must be a AliAODRecoDecayHF2Prong object
875 // returns 0 if both D0 and D0bar are rejectecd
876 // 1 if D0 is accepted while D0bar is rejected
877 // 2 if D0bar is accepted while D0 is rejected
878 // 3 if both are accepted
879 // fWhyRejection variable (not returned for the moment, print it if needed)
880 // keeps some information on why a candidate has been
881 // rejected according to the following (unfriendly?) scheme
882 // if more rejection cases are considered interesting, just add numbers
883 //
884 // TO BE CONSIDERED WITH A GRAIN OF SALT (the order in which cut are applied is relevant)
885 // from 20 to 30: "detector" selection (PID acceptance)
886 // 26: TPC refit
887 // 27: ITS refit
888 // 28: no (TOF||TPC) pid information (no kTOFpid,kTOFout,kTIME,kTPCpid,...)
889 //
890 // from 30 to 40: PID selection
891 // 31: no Kaon compatible tracks found between daughters
892 // 32: no Kaon identified tracks found (strong sel. at low momenta)
893 // 33: both mass hypotheses are rejected
894 //
895 // ############################################################
896
897 if(!fUsePID) return 3;
898 fWhyRejection=0;
899 Int_t isD0D0barPID[2]={1,2};
900 Double_t nsigmaTPCpi=-1., nsigmaTPCK=-1.; //used for TPC pid
901 Double_t tofSig,times[5];// used fot TOF pid
902 Int_t hasPID[2]={2,2};// flag to count how many detectors give PID info for the daughters
903 Int_t isKaonPionTOF[2][2],isKaonPionTPC[2][2];
904 Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
905 // same for prong 2
906 // values convention -1 = discarded
907 // 0 = not identified (but compatible) || No PID (->hasPID flag)
908 // 1 = identified
909 // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
910 // Initial hypothesis: unknwon (but compatible)
911 isKaonPionTOF[0][0]=0;
912 isKaonPionTOF[0][1]=0;
913 isKaonPionTOF[1][0]=0;
914 isKaonPionTOF[1][1]=0;
915
916 isKaonPionTPC[0][0]=0;
917 isKaonPionTPC[0][1]=0;
918 isKaonPionTPC[1][0]=0;
919 isKaonPionTPC[1][1]=0;
920
921 combinedPID[0][0]=0;
922 combinedPID[0][1]=0;
923 combinedPID[1][0]=0;
924 combinedPID[1][1]=0;
925
926
927
928 for(Int_t daught=0;daught<2;daught++){
929 //Loop con prongs
930
931 // ########### Step 0- CHECKING minimal PID "ACCEPTANCE" ####################
932
933 AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
934
935 if(!(aodtrack->GetStatus()&AliESDtrack::kTPCrefit)){
936 fWhyRejection=26;
937 return 0;
938 }
939 if(!(aodtrack->GetStatus()&AliESDtrack::kITSrefit)){
940 fWhyRejection=27;
941 return 0;
942 }
943
944 AliAODPid *pid=aodtrack->GetDetPid();
945 if(!pid) {
946 //delete esdtrack;
947 hasPID[daught]--;
948 continue;
949 }
950
951 // ########### Step 1- Check of TPC and TOF response ####################
952
953 Double_t ptrack=aodtrack->P();
954 //#################### TPC PID #######################
955 if (!(aodtrack->GetStatus()&AliESDtrack::kTPCpid )){
956 // NO TPC PID INFO FOR THIS TRACK
957 hasPID[daught]--;
958 }
959 else {
960 static AliTPCPIDResponse theTPCpid;
961 AliAODPid *pidObj = aodtrack->GetDetPid();
962 Double_t ptProng=pidObj->GetTPCmomentum();
963 nsigmaTPCpi = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kPion);
964 nsigmaTPCK = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kKaon);
965 //if(ptrack<0.6){
966 if(ptProng<0.6){
967 if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
968 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
969 if(TMath::Abs(nsigmaTPCpi)<2.)isKaonPionTPC[daught][1]=1;
970 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
971 }
972 //else if(ptrack<.8){
973 else if(ptProng<.8){
974 if(TMath::Abs(nsigmaTPCK)<1.)isKaonPionTPC[daught][0]=1;
975 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
976 if(TMath::Abs(nsigmaTPCpi)<1.)isKaonPionTPC[daught][1]=1;
977 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
978 }
979 else {
980 // if(nsigmaTPCK>-2.&&nsigmaTPCK<1.)isKaonPionTPC[daught][0]=1;
981 if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
982 //if(nsigmaTPCpi>-1.&&nsigmaTPCpi<2.)isKaonPionTPC[daught][1]=1;
983 if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
984 }
985 }
986
987
988 // ##### TOF PID: do not ask nothing for pion/protons ############
989 if(!((aodtrack->GetStatus()&AliESDtrack::kTOFpid)&&(aodtrack->GetStatus()&AliESDtrack::kTOFout)&&(aodtrack->GetStatus()&AliESDtrack::kTIME))){
990 // NO TOF PID INFO FOR THIS TRACK
991 hasPID[daught]--;
992 }
993 else{
994 tofSig=pid->GetTOFsignal();
995 pid->GetIntegratedTimes(times);
996 if((tofSig-times[3])>5.*160.)return 0;// PROTON REJECTION
997 if(TMath::Abs(tofSig-times[3])>3.*160.){
998 isKaonPionTOF[daught][0]=-1;
999 }
1000 else {
1001 if(ptrack<1.5){
1002 isKaonPionTOF[daught][0]=1;
1003 }
1004 }
1005 }
1006
1007 //######### Step 2: COMBINE TOF and TPC PID ###############
1008 // we apply the following convention: if TPC and TOF disagree (discarded Vs identified) -> unknown
1009 combinedPID[daught][0]=isKaonPionTOF[daught][0]+isKaonPionTPC[daught][0];
1010 combinedPID[daught][1]=isKaonPionTOF[daught][1]+isKaonPionTPC[daught][1];
1011
1012
1013 //######### Step 3: USE PID INFO
1014
1015 if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){// if not a K- and not a pi- both D0 and D0bar excluded
1016 isD0D0barPID[0]=0;
1017 isD0D0barPID[1]=0;
1018 }
1019 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
1020 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
1021 else isD0D0barPID[0]=0;// if K+ D0 excluded
1022 }
1023 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
1024 isD0D0barPID[0]=0;
1025 isD0D0barPID[1]=0;
1026 }
1027 else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
1028 if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
1029 else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
1030 }
1031 else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
1032 if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
1033 else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
1034 }
1035
1036 // ########## ALSO DIFFERENT TPC PID REQUEST FOR LOW pt D0: request of K identification ###############################
1037 // ########## more tolerant criteria for single particle ID-> more selective criteria for D0 ##############################
1038 // ############### NOT OPTIMIZED YET ###################################
1039 if(d->Pt()<2.){
1040 isKaonPionTPC[daught][0]=0;
1041 isKaonPionTPC[daught][1]=0;
1042 AliAODPid *pidObj = aodtrack->GetDetPid();
1043 Double_t ptProng=pidObj->GetTPCmomentum();
1044 //if(ptrack<0.6){
1045 if(ptProng<0.6){
1046 if(TMath::Abs(nsigmaTPCK)<3.)isKaonPionTPC[daught][0]=1;
1047 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1048 if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
1049 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1050 }
1051 //else if(ptrack<.8){
1052 else if(ptProng<.8){
1053 if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
1054 else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1055 if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
1056 else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1057 }
1058 else {
1059 if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1060 if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1061 }
1062 }
1063
1064 }// END OF LOOP ON DAUGHTERS
1065
1066 // FURTHER PID REQUEST (both daughter info is needed)
1067 if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
1068 fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
1069 return 0;
1070 }
1071 else if(hasPID[0]==0&&hasPID[1]==0){
1072 fWhyRejection=28;// reject cases in which no PID info is available
1073 return 0;
1074 }
1075 if(d->Pt()<2.){
1076 // request of K identification at low D0 pt
1077 combinedPID[0][0]=0;
1078 combinedPID[0][1]=0;
1079 combinedPID[1][0]=0;
1080 combinedPID[1][1]=0;
1081
1082 combinedPID[0][0]=isKaonPionTOF[0][0]+isKaonPionTPC[0][0];
1083 combinedPID[0][1]=isKaonPionTOF[0][1]+isKaonPionTPC[0][1];
1084 combinedPID[1][0]=isKaonPionTOF[1][0]+isKaonPionTPC[1][0];
1085 combinedPID[1][1]=isKaonPionTOF[1][1]+isKaonPionTPC[1][1];
1086
1087 if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
1088 fWhyRejection=32;// reject cases where the Kaon is not identified
1089 return 0;
1090 }
1091 }
1092
1093 // cout<<"Why? "<<fWhyRejection<<endl;
1094 return isD0D0barPID[0]+isD0D0barPID[1];
1095}
1096
1097
1098
1099//---------------------------------------------------------------------------
1100Int_t AliRDHFCutsD0toKpi::CombineSelectionLevels(Int_t selectionvalTrack,
1101 Int_t selectionvalCand,
1102 Int_t selectionvalPID) const
1103{
1104 //
1105 // This method combines the tracks, PID and cuts selection results
1106 //
1107 if(selectionvalTrack==0) return 0;
1108
1109 Int_t returnvalue;
1110
1111 switch(selectionvalPID) {
1112 case 0:
1113 returnvalue=0;
1114 break;
1115 case 1:
1116 returnvalue=((selectionvalCand==1 || selectionvalCand==3) ? 1 : 0);
1117 break;
1118 case 2:
1119 returnvalue=((selectionvalCand==2 || selectionvalCand==3) ? 2 : 0);
1120 break;
1121 case 3:
1122 returnvalue=selectionvalCand;
1123 break;
1124 default:
1125 returnvalue=0;
1126 break;
1127 }
1128
1129 return returnvalue;
1130}
1131//----------------------------------------------------------------------------
1132Int_t AliRDHFCutsD0toKpi::IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
1133{
1134 //
1135 // Note: this method is temporary
1136 // Additional cuts on decay lenght and lower cut for d0 norm are applied using vertex without candidate's daughters
1137 //
1138
1139 //apply cuts
1140
1141 Float_t normDecLengthCut=1.,decLengthCut=TMath::Min(d->P()*0.0066+0.01,0.06/*cm*/), normd0Cut=0.5;
1142 // "decay length" expo law with tau' = beta*gamma*ctau= p/m*ctau =p*0.0123/1.864~p*0.0066
1143 // decay lenght > ctau' implies to retain (1-1/e) (for signal without considering detector resolution),
1144
1145 Int_t returnvalue=3; //cut passed
1146 for(Int_t i=0;i<2/*prongs*/;i++){
1147 if(TMath::Abs(d->Normalizedd0Prong(i))<normd0Cut) return 0; //normd0Cut not passed
1148 }
1149 if(d->DecayLength2()<decLengthCut*decLengthCut) return 0; //decLengthCut not passed
1150 if(d->NormalizedDecayLength2()<normDecLengthCut*normDecLengthCut) return 0; //decLengthCut not passed
1151
1152 return returnvalue;
1153}
1154
1155//----------------------------------------------
1156void AliRDHFCutsD0toKpi::SetUseKF(Bool_t useKF)
1157{
1158 //switch on candidate selection via AliKFparticle
1159 if(!useKF) return;
1160 if(useKF){
1161 fUseKF=useKF;
1162 Int_t nvarsKF=11;
1163 SetNVars(nvarsKF);
1164 TString varNamesKF[11]={"inv. mass [GeV]",
1165 "dca [cm]",
1166 "cosThetaStar",
1167 "pTK [GeV/c]",
1168 "pTPi [GeV/c]",
1169 "d0K [cm]",
1170 "d0Pi [cm]",
1171 "d0d0 [cm^2]",
1172 "cosThetaPoint"
1173 "DecayLength[cm]",
1174 "RedChi2"};
1175 Bool_t isUpperCutKF[11]={kTRUE,
1176 kTRUE,
1177 kTRUE,
1178 kFALSE,
1179 kFALSE,
1180 kTRUE,
1181 kTRUE,
1182 kTRUE,
1183 kFALSE,
1184 kFALSE,
1185 kTRUE};
1186 SetVarNames(nvarsKF,varNamesKF,isUpperCutKF);
1187 SetGlobalIndex();
1188 Bool_t forOpt[11]={kFALSE,
1189 kTRUE,
1190 kTRUE,
1191 kFALSE,
1192 kFALSE,
1193 kFALSE,
1194 kFALSE,
1195 kTRUE,
1196 kTRUE,
1197 kFALSE,
1198 kFALSE};
1199 SetVarsForOpt(4,forOpt);
1200 }
1201 return;
1202}
1203
1204
e2aa82b6 1205//---------------------------------------------------------------------------
fc051756 1206void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {
1207 //
1208 //STANDARD CUTS USED FOR 2010 pp analysis
1209 //dca cut will be enlarged soon to 400 micron
1210 //
1211
1212 SetName("D0toKpiCutsStandard");
1213 SetTitle("Standard Cuts for D0 analysis");
1214
1215 // PILE UP REJECTION
1216 SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1217
1218 // EVENT CUTS
1219 SetMinVtxContr(1);
1220 // MAX Z-VERTEX CUT
1221 SetMaxVtxZ(10.);
1222
1223 // TRACKS ON SINGLE TRACKS
1224 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1225 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1226 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1227 esdTrackCuts->SetRequireITSRefit(kTRUE);
1228 // esdTrackCuts->SetMinNClustersITS(4);
1229 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1230 esdTrackCuts->SetMinDCAToVertexXY(0.);
1231 esdTrackCuts->SetEtaRange(-0.8,0.8);
1232 esdTrackCuts->SetPtRange(0.3,1.e10);
1233
1234 AddTrackCuts(esdTrackCuts);
57d15652 1235 delete esdTrackCuts;
1236 esdTrackCuts=NULL;
fc051756 1237
1238 const Int_t nptbins =14;
1239 const Double_t ptmax = 9999.;
1240 const Int_t nvars=11;
1241 Float_t ptbins[nptbins+1];
1242 ptbins[0]=0.;
1243 ptbins[1]=0.5;
1244 ptbins[2]=1.;
1245 ptbins[3]=2.;
1246 ptbins[4]=3.;
1247 ptbins[5]=4.;
1248 ptbins[6]=5.;
1249 ptbins[7]=6.;
1250 ptbins[8]=7.;
1251 ptbins[9]=8.;
1252 ptbins[10]=12.;
1253 ptbins[11]=16.;
1254 ptbins[12]=20.;
1255 ptbins[13]=24.;
1256 ptbins[14]=ptmax;
1257
1258 SetGlobalIndex(nvars,nptbins);
1259 SetPtBins(nptbins+1,ptbins);
1260
1261 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*/
1262 {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*/
1263 {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 */
1264 {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 */
1265 {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 */
1266 {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 */
1267 {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 */
1268 {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 */
1269 {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 */
1270 {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 */
1271 {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 */
1272 {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 */
1273 {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 */
1274 {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 */
1275
1276
1277 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1278 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1279 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1280
1281 for (Int_t ibin=0;ibin<nptbins;ibin++){
1282 for (Int_t ivar = 0; ivar<nvars; ivar++){
1283 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1284 }
1285 }
1286
1287 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1288 SetUseSpecialCuts(kTRUE);
1289 SetRemoveDaughtersFromPrim(kTRUE);
1290
1291 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1292 delete [] cutsMatrixTransposeStand;
1293 cutsMatrixTransposeStand=NULL;
1294
1295 // PID SETTINGS
1296 AliAODPidHF* pidObj=new AliAODPidHF();
1297 //pidObj->SetName("pid4D0");
1298 Int_t mode=1;
1299 const Int_t nlims=2;
1300 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1301 Bool_t compat=kTRUE; //effective only for this mode
1302 Bool_t asym=kTRUE;
1303 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1304 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1305 pidObj->SetMatch(mode);
1306 pidObj->SetPLimit(plims,nlims);
1307 pidObj->SetSigma(sigmas);
1308 pidObj->SetCompat(compat);
1309 pidObj->SetTPC(kTRUE);
1310 pidObj->SetTOF(kTRUE);
1311 pidObj->SetPCompatTOF(1.5);
1312 pidObj->SetSigmaForTPCCompat(3.);
1313 pidObj->SetSigmaForTOFCompat(3.);
ec0f397b 1314 pidObj->SetOldPid(kTRUE);
fc051756 1315
1316 SetPidHF(pidObj);
1317 SetUsePID(kTRUE);
1318 SetUseDefaultPID(kFALSE);
1319
1320 SetLowPt(kFALSE);
1321
1322 PrintAll();
1323
1324 delete pidObj;
1325 pidObj=NULL;
1326
1327 return;
1328
1329}
1330
e0dfdb1a 1331//___________________________________________________________________________
1332void AliRDHFCutsD0toKpi::SetStandardCutsPP2010vsMult() {
1333 //
1334 // Cuts for 2010 pp 7 TeV data analysis in Ntracklets bins (vs multiplicity)
1335 //
1336 SetName("D0toKpiCuts");
1337 SetTitle("Cuts for D0 analysis in 2010-data pp 7 TeV vs multiplicity");
1338
1339 //
1340 // Track cuts
1341 //
1342 AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
1343 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1344 //default
1345 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1346 esdTrackCuts->SetRequireITSRefit(kTRUE);
1347 esdTrackCuts->SetEtaRange(-0.8,0.8);
1348 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1349 AliESDtrackCuts::kAny);
1350 // default is kBoth, otherwise kAny
1351 esdTrackCuts->SetMinDCAToVertexXY(0.);
1352 esdTrackCuts->SetPtRange(0.3,1.e10);
1353
1354 AddTrackCuts(esdTrackCuts);
1355 delete esdTrackCuts;
1356 esdTrackCuts=NULL;
1357
1358
1359 //
1360 // Cut values per pt bin
1361 //
1362 const Int_t nvars=11;
1363 const Int_t nptbins=14;
1364 Float_t* ptbins;
1365 ptbins=new Float_t[nptbins+1];
1366 ptbins[0]=0.;
1367 ptbins[1]=0.5;
1368 ptbins[2]=1.;
1369 ptbins[3]=2.;
1370 ptbins[4]=3.;
1371 ptbins[5]=4.;
1372 ptbins[6]=5.;
1373 ptbins[7]=6.;
1374 ptbins[8]=7.;
1375 ptbins[9]=8.;
1376 ptbins[10]=12.;
1377 ptbins[11]=16.;
1378 ptbins[12]=20.;
1379 ptbins[13]=24.;
1380 ptbins[14]=9999.;
1381
1382 SetPtBins(nptbins+1,ptbins);
1383
1384 //setting cut values
1385 Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,350.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,3.},/* pt<0.5*/
1386 {0.400,350.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,3.},/* 0.5<pt<1*/
1387 {0.400,300.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-30000.*1E-8,0.80,0.,4.},/* 1<pt<2 */
1388 {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-20000.*1E-8,0.85,0.,4.},/* 2<pt<3 */
1389 {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-15000.*1E-8,0.85,0.,4.},/* 3<pt<4 */
1390 {0.400,300.*1E-4,0.75,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.875,0.,0.},/* 4<pt<5 */
1391 {0.400,300.*1E-4,0.75,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.875,0.,0.},/* 5<pt<6 */
1392 {0.400,400.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 6<pt<7 */
1393 {0.400,400.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 7<pt<8 */
1394 {0.400,0.06,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.00001,0.85,0.,0.},/* 8<pt<12 */
1395 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 12<pt<16 */
1396 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 16<pt<20 */
1397 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 20<pt<24 */
1398 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.}};/* pt>24 */
1399
1400
1401 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1402 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1403 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1404
1405 for (Int_t ibin=0;ibin<nptbins;ibin++){
1406 for (Int_t ivar = 0; ivar<nvars; ivar++){
1407 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1408 }
1409 }
1410
1411 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1412 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1413 delete [] cutsMatrixTransposeStand;
1414
1415 SetUseSpecialCuts(kTRUE);
1416 SetMaximumPtSpecialCuts(8.);
1417
1418 //Do recalculate the vertex
1419 SetRemoveDaughtersFromPrim(kTRUE);
1420
1421 //
1422 // Pid settings
1423 //
1424 Bool_t pidflag=kTRUE;
1425 SetUsePID(pidflag);
1426 if(pidflag) cout<<"PID is used"<<endl;
1427 else cout<<"PID is not used"<<endl;
1428 //
1429 AliAODPidHF* pidObj=new AliAODPidHF();
1430 Int_t mode=1;
1431 const Int_t nlims=2;
1432 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1433 Bool_t compat=kTRUE; //effective only for this mode
1434 Bool_t asym=kTRUE;
1435 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1436 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1437 pidObj->SetMatch(mode);
1438 pidObj->SetPLimit(plims,nlims);
1439 pidObj->SetSigma(sigmas);
1440 pidObj->SetCompat(compat);
1441 pidObj->SetTPC(kTRUE);
1442 pidObj->SetTOF(kTRUE);
1443 pidObj->SetPCompatTOF(1.5);
1444 pidObj->SetSigmaForTPCCompat(3.);
1445 pidObj->SetSigmaForTOFCompat(3.);
1446 pidObj->SetOldPid(kTRUE);
1447 // pidObj->SetOldPid(kFALSE);
1448 SetPidHF(pidObj);
1449
1450 SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
1451
1452 SetLowPt(kFALSE);
1453
1454 //activate pileup rejection (for pp)
1455 SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1456
1457
1458 PrintAll();
1459
1460 delete pidObj;
1461 pidObj=NULL;
1462 delete [] ptbins;
1463 ptbins=NULL;
1464
1465 return;
1466}
fc051756 1467
e2aa82b6 1468//---------------------------------------------------------------------------
fc051756 1469void AliRDHFCutsD0toKpi::SetStandardCutsPP2011_276TeV() {
1470 //
1471 // STANDARD CUTS USED FOR 2011 pp analysis at 2.76TeV
1472 //
1473
1474 SetName("D0toKpiCutsStandard");
1475 SetTitle("Standard Cuts for D0 analysis in pp2011 at 2.76TeV run");
1476
1477 //
1478 // Track cuts
1479 //
1480 AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
1481 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1482 //default
1483 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1484 esdTrackCuts->SetRequireITSRefit(kTRUE);
1485 esdTrackCuts->SetEtaRange(-0.8,0.8);
1486 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1487 AliESDtrackCuts::kAny);
1488 // default is kBoth, otherwise kAny
1489 esdTrackCuts->SetMinDCAToVertexXY(0.);
1490 esdTrackCuts->SetPtRange(0.3,1.e10);
1491
1492 esdTrackCuts->SetMaxDCAToVertexXY(1.);
1493 esdTrackCuts->SetMaxDCAToVertexZ(1.);
1494 esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1495
1496 AddTrackCuts(esdTrackCuts);
57d15652 1497 delete esdTrackCuts;
1498 esdTrackCuts=NULL;
fc051756 1499
1500
1501 const Int_t nvars=11;
1502 const Int_t nptbins=13;
d1029f2b 1503 Float_t ptbins[nptbins+1];
fc051756 1504 ptbins[0]=0.;
1505 ptbins[1]=0.5;
1506 ptbins[2]=1.;
1507 ptbins[3]=2.;
1508 ptbins[4]=3.;
1509 ptbins[5]=4.;
1510 ptbins[6]=5.;
1511 ptbins[7]=6.;
1512 ptbins[8]=8.;
1513 ptbins[9]=12.;
1514 ptbins[10]=16.;
1515 ptbins[11]=20.;
1516 ptbins[12]=24.;
1517 ptbins[13]=9999.;
1518
1519 SetPtBins(nptbins+1,ptbins);
1520
fc051756 1521
1522 Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,0.04,0.75,0.3,0.3,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* pt<0.5*/
1523 {0.400,0.04,0.75,0.3,0.3,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* 0.5<pt<1*/
1524 {0.400,0.03,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.8,0.,0.},/* 1<pt<2 */
1525 {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.0003,0.9,0.,0.},/* 2<pt<3 */
1526 {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.0002,0.9,0.,0.},/* 3<pt<4 */
1527 {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.00015,0.9,0.,0.},/* 4<pt<5 */
1528 {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.0001,0.9,0.,0.},/* 5<pt<6 */
1529 {0.400,0.09,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* 6<pt<8 */
1530 {0.400,0.06,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.00001,0.85,0.,0.},/* 8<pt<12 */
1531 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 12<pt<16 */
1532 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 16<pt<20 */
1533 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 20<pt<24 */
1534 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.}};/* pt>24 */
1535
1536 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1537 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1538 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1539 for (Int_t ibin=0;ibin<nptbins;ibin++){
1540 for (Int_t ivar = 0; ivar<nvars; ivar++){
1541 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1542 }
1543 }
1544 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
d1029f2b 1545 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1546 delete [] cutsMatrixTransposeStand;
fc051756 1547
1548
1549 //pid settings
1550 AliAODPidHF* pidObj=new AliAODPidHF();
1551 Int_t mode=1;
1552 const Int_t nlims=2;
1553 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1554 Bool_t compat=kTRUE; //effective only for this mode
1555 Bool_t asym=kTRUE;
1556 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1557 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1558 pidObj->SetMatch(mode);
1559 pidObj->SetPLimit(plims,nlims);
1560 pidObj->SetSigma(sigmas);
1561 pidObj->SetCompat(compat);
1562 pidObj->SetTPC(kTRUE);
1563 pidObj->SetTOF(kTRUE);
ec0f397b 1564 pidObj->SetOldPid(kTRUE);
fc051756 1565 SetPidHF(pidObj);
1566
1567 SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
1568
1569 SetUsePID(kTRUE);
1570
1571 //activate pileup rejection (for pp)
1572 SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1573
1574 //Do recalculate the vertex
1575 SetRemoveDaughtersFromPrim(kTRUE);
1576
1577 // Cut on the zvtx
1578 SetMaxVtxZ(10.);
1579
1580 // Use the kFirst selection for tracks with candidate D with pt<2
1581 SetSelectCandTrackSPDFirst(kTRUE,2.);
1582
1583 // Use special cuts for D candidates with pt<2
1584 SetUseSpecialCuts(kTRUE);
1585 SetMaximumPtSpecialCuts(2.);
1586
1587 PrintAll();
1588
1589 delete pidObj;
1590 pidObj=NULL;
1591
1592 return;
1593}
1594
1595
e2aa82b6 1596//---------------------------------------------------------------------------
fc051756 1597void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
1598 //
1599 // Final CUTS USED FOR 2010 PbPb analysis of central events
1600 // REMEMBER TO EVENTUALLY SWITCH ON
1601 // SetUseAOD049(kTRUE);
1602 //
1603
1604 SetName("D0toKpiCutsStandard");
1605 SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
1606
1607
1608 // CENTRALITY SELECTION
1609 SetMinCentrality(0.);
1610 SetMaxCentrality(80.);
1611 SetUseCentrality(AliRDHFCuts::kCentV0M);
1612
1613
1614
1615 // EVENT CUTS
1616 SetMinVtxContr(1);
1617 // MAX Z-VERTEX CUT
1618 SetMaxVtxZ(10.);
1619 // PILE UP
1620 SetOptPileup(AliRDHFCuts::kNoPileupSelection);
1621 // PILE UP REJECTION
1622 //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1623
1624 // TRACKS ON SINGLE TRACKS
1625 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1626 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1627 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1628 esdTrackCuts->SetRequireITSRefit(kTRUE);
1629 // esdTrackCuts->SetMinNClustersITS(4);
1630 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1631 esdTrackCuts->SetMinDCAToVertexXY(0.);
1632 esdTrackCuts->SetEtaRange(-0.8,0.8);
1633 esdTrackCuts->SetPtRange(0.7,1.e10);
1634
1635 esdTrackCuts->SetMaxDCAToVertexXY(1.);
1636 esdTrackCuts->SetMaxDCAToVertexZ(1.);
1637 esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1638
1639
1640 AddTrackCuts(esdTrackCuts);
57d15652 1641 delete esdTrackCuts;
1642 esdTrackCuts=NULL;
fc051756 1643 // SPD k FIRST for Pt<3 GeV/c
1644 SetSelectCandTrackSPDFirst(kTRUE, 3);
1645
1646 // CANDIDATE CUTS
1647 const Int_t nptbins =13;
1648 const Double_t ptmax = 9999.;
1649 const Int_t nvars=11;
1650 Float_t ptbins[nptbins+1];
1651 ptbins[0]=0.;
1652 ptbins[1]=0.5;
1653 ptbins[2]=1.;
1654 ptbins[3]=2.;
1655 ptbins[4]=3.;
1656 ptbins[5]=4.;
1657 ptbins[6]=5.;
1658 ptbins[7]=6.;
1659 ptbins[8]=8.;
1660 ptbins[9]=12.;
1661 ptbins[10]=16.;
1662 ptbins[11]=20.;
1663 ptbins[12]=24.;
1664 ptbins[13]=ptmax;
1665
1666 SetGlobalIndex(nvars,nptbins);
1667 SetPtBins(nptbins+1,ptbins);
1668 SetMinPtCandidate(2.);
1669
1670 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.99,2.},/* pt<0.5*/
1671 {0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.9,0.99,2.},/* 0.5<pt<1*/
dab01ae7 1672 {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-43000.*1E-8,0.85,0.995,8.},/* 1<pt<2 */
fc051756 1673 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-45000.*1E-8,0.95,0.998,7.},/* 2<pt<3 */
1674 {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 */
1675 {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 */
1676 {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 */
1677 {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 */
1678 {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 */
1679 {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.83,0.998,8.},/* 12<pt<16 */
1680 {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.82,0.995,6.},/* 16<pt<20 */
1681 {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.81,0.995,6.},/* 20<pt<24 */
1682 {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.8,0.99,2.}};/* pt>24 */
1683
1684
1685 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1686 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1687 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1688
1689 for (Int_t ibin=0;ibin<nptbins;ibin++){
1690 for (Int_t ivar = 0; ivar<nvars; ivar++){
1691 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1692 }
1693 }
1694
1695 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1696 SetUseSpecialCuts(kTRUE);
1697 SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1698 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1699 delete [] cutsMatrixTransposeStand;
1700 cutsMatrixTransposeStand=NULL;
1701
1702 // PID SETTINGS
1703 AliAODPidHF* pidObj=new AliAODPidHF();
1704 //pidObj->SetName("pid4D0");
1705 Int_t mode=1;
1706 const Int_t nlims=2;
1707 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1708 Bool_t compat=kTRUE; //effective only for this mode
1709 Bool_t asym=kTRUE;
1710 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1711 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1712 pidObj->SetMatch(mode);
1713 pidObj->SetPLimit(plims,nlims);
1714 pidObj->SetSigma(sigmas);
1715 pidObj->SetCompat(compat);
1716 pidObj->SetTPC(kTRUE);
1717 pidObj->SetTOF(kTRUE);
1718 pidObj->SetPCompatTOF(2.);
1719 pidObj->SetSigmaForTPCCompat(3.);
1720 pidObj->SetSigmaForTOFCompat(3.);
ec0f397b 1721 pidObj->SetOldPid(kTRUE);
fc051756 1722
1723
1724 SetPidHF(pidObj);
1725 SetUsePID(kTRUE);
1726 SetUseDefaultPID(kFALSE);
1727
1728
1729 PrintAll();
1730
1731
1732 delete pidObj;
1733 pidObj=NULL;
1734
1735 return;
1736
1737}
1738
e2aa82b6 1739//---------------------------------------------------------------------------
fc051756 1740void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() {
1741 // CUTS USED FOR D0 RAA for the CENTRALITY RANGE 20-80%
1742
1743
1744 SetName("D0toKpiCutsStandard");
1745 SetTitle("Standard Cuts for D0 analysis in PbPb2010 run, for peripheral events");
1746
1747
1748 // CENTRALITY SELECTION
1749 SetMinCentrality(40.);
1750 SetMaxCentrality(80.);
1751 SetUseCentrality(AliRDHFCuts::kCentV0M);
1752
1753 // EVENT CUTS
1754 SetMinVtxContr(1);
1755 // MAX Z-VERTEX CUT
1756 SetMaxVtxZ(10.);
1757 // PILE UP
1758 SetOptPileup(AliRDHFCuts::kNoPileupSelection);
1759 // PILE UP REJECTION
1760 //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1761
1762 // TRACKS ON SINGLE TRACKS
1763 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1764 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1765 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1766 esdTrackCuts->SetRequireITSRefit(kTRUE);
1767 // esdTrackCuts->SetMinNClustersITS(4);
1768 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1769 esdTrackCuts->SetMinDCAToVertexXY(0.);
1770 esdTrackCuts->SetEtaRange(-0.8,0.8);
1771 esdTrackCuts->SetPtRange(0.5,1.e10);
1772
1773 esdTrackCuts->SetMaxDCAToVertexXY(1.);
1774 esdTrackCuts->SetMaxDCAToVertexZ(1.);
1775 esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0025*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1776
1777
1778 AddTrackCuts(esdTrackCuts);
57d15652 1779 delete esdTrackCuts;
1780 esdTrackCuts=NULL;
fc051756 1781 // SPD k FIRST for Pt<3 GeV/c
1782 SetSelectCandTrackSPDFirst(kTRUE, 3);
1783
1784 // CANDIDATE CUTS
1785 const Int_t nptbins =13;
1786 const Double_t ptmax = 9999.;
1787 const Int_t nvars=11;
1788 Float_t ptbins[nptbins+1];
1789 ptbins[0]=0.;
1790 ptbins[1]=0.5;
1791 ptbins[2]=1.;
1792 ptbins[3]=2.;
1793 ptbins[4]=3.;
1794 ptbins[5]=4.;
1795 ptbins[6]=5.;
1796 ptbins[7]=6.;
1797 ptbins[8]=8.;
1798 ptbins[9]=12.;
1799 ptbins[10]=16.;
1800 ptbins[11]=20.;
1801 ptbins[12]=24.;
1802 ptbins[13]=ptmax;
1803
1804 SetGlobalIndex(nvars,nptbins);
1805 SetPtBins(nptbins+1,ptbins);
1806 SetMinPtCandidate(0.);
1807
1808 Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-20000.*1E-8,0.7,0.993,2.},/* pt<0.5*/
1809 {0.400,400.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.85,0.993,2.},/* 0.5<pt<1*/
1810 {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.85,0.995,6.},/* 1<pt<2 */
1811 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.95,0.991,5.},/* 2<pt<3 */
1812 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-36000.*1E-8,0.95,0.993,5.},/* 3<pt<4 */
1813 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-27000.*1E-8,0.95,0.995,5.},/* 4<pt<5 */
1814 {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 */
1815 {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 */
1816 {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.995,5.},/* 8<pt<12 */
1817 {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.83,0.995,4.},/* 12<pt<16 */
1818 {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.82,0.995,4.},/* 16<pt<20 */
1819 {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.81,0.995,4.},/* 20<pt<24 */
1820 {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.8,0.995,4.}};/* pt>24 */
1821
1822
1823 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1824 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1825 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1826
1827 for (Int_t ibin=0;ibin<nptbins;ibin++){
1828 for (Int_t ivar = 0; ivar<nvars; ivar++){
1829 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1830 }
1831 }
1832
1833 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1834 SetUseSpecialCuts(kTRUE);
1835 SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1836 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1837 delete [] cutsMatrixTransposeStand;
1838 cutsMatrixTransposeStand=NULL;
1839
1840 // PID SETTINGS
1841 AliAODPidHF* pidObj=new AliAODPidHF();
1842 //pidObj->SetName("pid4D0");
1843 Int_t mode=1;
1844 const Int_t nlims=2;
1845 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1846 Bool_t compat=kTRUE; //effective only for this mode
1847 Bool_t asym=kTRUE;
1848 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1849 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1850 pidObj->SetMatch(mode);
1851 pidObj->SetPLimit(plims,nlims);
1852 pidObj->SetSigma(sigmas);
1853 pidObj->SetCompat(compat);
1854 pidObj->SetTPC(kTRUE);
1855 pidObj->SetTOF(kTRUE);
1856 pidObj->SetPCompatTOF(2.);
1857 pidObj->SetSigmaForTPCCompat(3.);
1858 pidObj->SetSigmaForTOFCompat(3.);
ec0f397b 1859 pidObj->SetOldPid(kTRUE);
fc051756 1860
1861 SetPidHF(pidObj);
1862 SetUsePID(kTRUE);
1863 SetUseDefaultPID(kFALSE);
1864
1865 SetLowPt(kTRUE,2.);
1866
1867 PrintAll();
1868
1869
1870 delete pidObj;
1871 pidObj=NULL;
1872
1873 return;
1874
1875}
1876
1877
e2aa82b6 1878//---------------------------------------------------------------------------
fc051756 1879void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() {
1880
1881 // Default 2010 PbPb cut object
1882 SetStandardCutsPbPb2010();
0ebb13ca 1883 AliAODPidHF *pidobj=GetPidHF();
1884
1885 pidobj->SetOldPid(kFALSE);
fc051756 1886
1887 //
1888 // Enable all 2011 PbPb run triggers
1889 //
1890 SetTriggerClass("");
1891 ResetMaskAndEnableMBTrigger();
1892 EnableCentralTrigger();
1893 EnableSemiCentralTrigger();
1894
1895}
e2aa82b6 1896
1897//---------------------------------------------------------------------------
1898Int_t AliRDHFCutsD0toKpi::IsSelectedCombPID(AliAODRecoDecayHF* d)
1899{
1900 // ############################################################
1901 //
1902 // Apply Bayesian PID selection
1903 // Implementation of Bayesian PID: Jeremy Wilkinson
1904 //
1905 // ############################################################
1906
1907 if (!fUsePID || !d) return 3;
1908
775d1f1c 1909
1910 if (fBayesianStrategy == kBayesWeightNoFilter) {
1911 //WeightNoFilter: Accept all particles (no PID cut) but fill mass histos with weights in task
1912 CalculateBayesianWeights(d);
1913 return 3;
6bb29c8c 1914 }
e2aa82b6 1915
1916
1917
1918 Int_t isD0 = 0;
1919 Int_t isD0bar = 0;
1920 Int_t returnvalue = 0;
1921
6bb29c8c 1922 Int_t isPosKaon = 0, isNegKaon = 0, isPosPion = 0, isNegPion = 0;
e2aa82b6 1923
775d1f1c 1924 //Bayesian methods used here check for ID of kaon, and whether it is positive or negative.
1925
1926
e2aa82b6 1927 Bool_t checkPIDInfo[2] = {kTRUE, kTRUE};
1928 AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
1929
1930 if ((aodtrack[0]->Charge()*aodtrack[1]->Charge()) != -1) {
1931 return 0; //reject if charges not opposite
1932 }
b8f433ab 1933 Double_t momentumpositive=0., momentumnegative=0.; //Used in "prob > prior" method
e2aa82b6 1934 for (Int_t daught = 0; daught < 2; daught++) {
1935 //Loop over prongs
1936
b8f433ab 1937 if (aodtrack[daught]->Charge() == -1) {
1938 momentumnegative = aodtrack[daught]->P();
1939 }
1940 if (aodtrack[daught]->Charge() == +1) {
1941 momentumpositive = aodtrack[daught]->P();
1942 }
1943 if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
1944 checkPIDInfo[daught] = kFALSE;
1945 continue;
1946 }
e2aa82b6 1947
b8f433ab 1948 }
e2aa82b6 1949
b8f433ab 1950 //Loop over daughters ends here
e2aa82b6 1951
b8f433ab 1952 if (!checkPIDInfo[0] && !checkPIDInfo[1]) {
1953 return 0; //Reject if both daughters lack both TPC and TOF info
1954 }
e2aa82b6 1955
e2aa82b6 1956
b8f433ab 1957 CalculateBayesianWeights(d); //Calculates all Bayesian probabilities for both positive and negative tracks
1958 // Double_t prob0[AliPID::kSPECIES];
1959 // fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), prob0);
775d1f1c 1960 ///Three possible Bayesian probability cuts: Picked using SetBayesianCondition(int).
b8f433ab 1961 switch (fBayesianCondition) {
775d1f1c 1962 ///A: Standard max. probability method (accept most likely species)
b8f433ab 1963 case kMaxProb:
6bb29c8c 1964 if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kKaon]) { //If highest probability lies with kaon
1965 isPosKaon = 1; //flag [daught] as a kaon
1966 }
1967
1968 if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kPion]) { //If highest probability lies with pion
1969 isPosPion = 1; //flag [daught] as a pion
1970 }
1971
1972 if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kKaon]) { //If highest probability lies with kaon
1973 isNegKaon = 1; //flag [daught] as a kaon
1974 }
1975
1976 if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kPion]) { //If highest probability lies with kaon
1977 isNegPion = 1; //flag [daught] as a pion
1978 }
1979
1980
1981 break;
1982 ///B: Accept if probability greater than prior
1983 case kAbovePrior:
1984
1985 if (fWeightsNegative[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
1986 GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumnegative)))) { //Retrieves relevant prior, gets value at momentum
1987 isNegKaon = 1;
1988 }
1989 if (fWeightsPositive[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
1990 GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumpositive)))) { //Retrieves relevant prior, gets value at momentum
1991 isPosKaon = 1;
1992 }
1993
1994 break;
1995
1996 ///C: Accept if probability greater than user-defined threshold
b8f433ab 1997 case kThreshold:
1998 if (fWeightsNegative[AliPID::kKaon] > fProbThreshold) {
1999 isNegKaon = 1;
2000 }
6bb29c8c 2001 if (fWeightsNegative[AliPID::kPion] > fProbThreshold) {
2002 isNegPion = 1;
2003 }
2004
b8f433ab 2005 if (fWeightsPositive[AliPID::kKaon] > fProbThreshold) {
2006 isPosKaon = 1;
2007 }
6bb29c8c 2008
2009 if (fWeightsPositive[AliPID::kPion] > fProbThreshold) {
2010 isPosPion = 1;
2011 }
b8f433ab 2012
2013 break;
6bb29c8c 2014 }
e2aa82b6 2015
775d1f1c 2016
2017 //Momentum-based selection (also applied to filtered weighted method)
2018
2019 if (fBayesianStrategy == kBayesMomentum || fBayesianCondition == kBayesWeight) {
2020 if (isNegKaon && isPosKaon) { // If both are kaons, reject
775d1f1c 2021 isD0 = 1;
6bb29c8c 2022 isD0bar = 1;
2023 } else if (isNegKaon && isPosPion) { //If negative kaon present, D0
2024 isD0 = 1;
2025 } else if (isPosKaon && isNegPion) { //If positive kaon present, D0bar
775d1f1c 2026 isD0bar = 1;
2027 } else { //If neither ID'd as kaon, subject to extra tests
2028 isD0 = 1;
2029 isD0bar = 1;
2030 if (aodtrack[0]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[0], "TOF"))) { //If it's low momentum and not ID'd as a kaon, we assume it must be a pion
2031 if (aodtrack[0]->Charge() == -1) {
2032 isD0 = 0; //Check charge and reject based on pion hypothesis
2033 }
2034 if (aodtrack[0]->Charge() == 1) {
2035 isD0bar = 0;
2036 }
2037 }
2038 if (aodtrack[1]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[1], "TOF"))) {
2039 if (aodtrack[1]->Charge() == -1) {
2040 isD0 = 0;
2041 }
2042 if (aodtrack[1]->Charge() == 1) {
2043 isD0bar = 0;
2044 }
2045 }
2046 }
e2aa82b6 2047
775d1f1c 2048
e2aa82b6 2049
775d1f1c 2050 if (isD0 && isD0bar) {
2051 returnvalue = 3;
b8f433ab 2052 }
775d1f1c 2053 if (isD0 && !isD0bar) {
2054 returnvalue = 1;
b8f433ab 2055 }
775d1f1c 2056 if (!isD0 && isD0bar) {
2057 returnvalue = 2;
b8f433ab 2058 }
775d1f1c 2059 if (!isD0 && !isD0bar) {
2060 returnvalue = 0;
b8f433ab 2061 }
775d1f1c 2062 }
b8f433ab 2063
775d1f1c 2064 //Simple Bayesian
2065 if (fBayesianStrategy == kBayesSimple) {
2066
731b87ba 2067 if (isPosKaon && isNegKaon) { //If both are ID'd as kaons, accept as possible
2068 returnvalue = 3;
6bb29c8c 2069 } else if (isNegKaon && isPosPion) { //If negative kaon, D0
775d1f1c 2070 returnvalue = 1;
6bb29c8c 2071 } else if (isPosKaon && isNegPion) { //If positive kaon, D0-bar
775d1f1c 2072 returnvalue = 2;
6bb29c8c 2073 } else if (isPosPion && isNegPion) {
775d1f1c 2074 returnvalue = 0; //If neither kaon, reject
6bb29c8c 2075 } else {returnvalue = 0;} //default
2076
775d1f1c 2077 }
2078
2079 return returnvalue;
e2aa82b6 2080
2081
e2aa82b6 2082
e2aa82b6 2083}
2084
2085
2086
e2aa82b6 2087//---------------------------------------------------------------------------
2088void AliRDHFCutsD0toKpi::CalculateBayesianWeights(AliAODRecoDecayHF* d)
2089
2090{
2091 //Function to compute weights for Bayesian method
e2aa82b6 2092
2093 AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
2094 if ((aodtrack[0]->Charge() * aodtrack[1]->Charge()) != -1) {
2095 return; //Reject if charges do not oppose
2096 }
2097 for (Int_t daught = 0; daught < 2; daught++) {
2098 //Loop over prongs
2099
2100 if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
2101 continue;
2102 }
2103
2104
c2ec6218 2105 // identify kaon, define weights
e2aa82b6 2106 if (aodtrack[daught]->Charge() == +1) {
c2ec6218 2107 fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsPositive);
e2aa82b6 2108 }
c2ec6218 2109
e2aa82b6 2110 if (aodtrack[daught]->Charge() == -1) {
c2ec6218 2111 fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsNegative);
e2aa82b6 2112 }
e2aa82b6 2113 }
2114}
2115
2116
2117