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