]> 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 'master' of https://git.cern.ch/reps/AliRoot
[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);
1235
1236
1237 const Int_t nptbins =14;
1238 const Double_t ptmax = 9999.;
1239 const Int_t nvars=11;
1240 Float_t ptbins[nptbins+1];
1241 ptbins[0]=0.;
1242 ptbins[1]=0.5;
1243 ptbins[2]=1.;
1244 ptbins[3]=2.;
1245 ptbins[4]=3.;
1246 ptbins[5]=4.;
1247 ptbins[6]=5.;
1248 ptbins[7]=6.;
1249 ptbins[8]=7.;
1250 ptbins[9]=8.;
1251 ptbins[10]=12.;
1252 ptbins[11]=16.;
1253 ptbins[12]=20.;
1254 ptbins[13]=24.;
1255 ptbins[14]=ptmax;
1256
1257 SetGlobalIndex(nvars,nptbins);
1258 SetPtBins(nptbins+1,ptbins);
1259
1260 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*/
1261 {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*/
1262 {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 */
1263 {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 */
1264 {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 */
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.},/* 4<pt<5 */
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.},/* 5<pt<6 */
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.},/* 6<pt<7 */
1268 {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 */
1269 {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 */
1270 {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 */
1271 {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 */
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.},/* 20<pt<24 */
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.}};/* pt>24 */
1274
1275
1276 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1277 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1278 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1279
1280 for (Int_t ibin=0;ibin<nptbins;ibin++){
1281 for (Int_t ivar = 0; ivar<nvars; ivar++){
1282 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1283 }
1284 }
1285
1286 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1287 SetUseSpecialCuts(kTRUE);
1288 SetRemoveDaughtersFromPrim(kTRUE);
1289
1290 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1291 delete [] cutsMatrixTransposeStand;
1292 cutsMatrixTransposeStand=NULL;
1293
1294 // PID SETTINGS
1295 AliAODPidHF* pidObj=new AliAODPidHF();
1296 //pidObj->SetName("pid4D0");
1297 Int_t mode=1;
1298 const Int_t nlims=2;
1299 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1300 Bool_t compat=kTRUE; //effective only for this mode
1301 Bool_t asym=kTRUE;
1302 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1303 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1304 pidObj->SetMatch(mode);
1305 pidObj->SetPLimit(plims,nlims);
1306 pidObj->SetSigma(sigmas);
1307 pidObj->SetCompat(compat);
1308 pidObj->SetTPC(kTRUE);
1309 pidObj->SetTOF(kTRUE);
1310 pidObj->SetPCompatTOF(1.5);
1311 pidObj->SetSigmaForTPCCompat(3.);
1312 pidObj->SetSigmaForTOFCompat(3.);
ec0f397b 1313 pidObj->SetOldPid(kTRUE);
fc051756 1314
1315 SetPidHF(pidObj);
1316 SetUsePID(kTRUE);
1317 SetUseDefaultPID(kFALSE);
1318
1319 SetLowPt(kFALSE);
1320
1321 PrintAll();
1322
1323 delete pidObj;
1324 pidObj=NULL;
1325
1326 return;
1327
1328}
1329
1330
e2aa82b6 1331//---------------------------------------------------------------------------
fc051756 1332void AliRDHFCutsD0toKpi::SetStandardCutsPP2011_276TeV() {
1333 //
1334 // STANDARD CUTS USED FOR 2011 pp analysis at 2.76TeV
1335 //
1336
1337 SetName("D0toKpiCutsStandard");
1338 SetTitle("Standard Cuts for D0 analysis in pp2011 at 2.76TeV run");
1339
1340 //
1341 // Track cuts
1342 //
1343 AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
1344 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1345 //default
1346 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1347 esdTrackCuts->SetRequireITSRefit(kTRUE);
1348 esdTrackCuts->SetEtaRange(-0.8,0.8);
1349 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1350 AliESDtrackCuts::kAny);
1351 // default is kBoth, otherwise kAny
1352 esdTrackCuts->SetMinDCAToVertexXY(0.);
1353 esdTrackCuts->SetPtRange(0.3,1.e10);
1354
1355 esdTrackCuts->SetMaxDCAToVertexXY(1.);
1356 esdTrackCuts->SetMaxDCAToVertexZ(1.);
1357 esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1358
1359 AddTrackCuts(esdTrackCuts);
1360
1361
1362 const Int_t nvars=11;
1363 const Int_t nptbins=13;
d1029f2b 1364 Float_t ptbins[nptbins+1];
fc051756 1365 ptbins[0]=0.;
1366 ptbins[1]=0.5;
1367 ptbins[2]=1.;
1368 ptbins[3]=2.;
1369 ptbins[4]=3.;
1370 ptbins[5]=4.;
1371 ptbins[6]=5.;
1372 ptbins[7]=6.;
1373 ptbins[8]=8.;
1374 ptbins[9]=12.;
1375 ptbins[10]=16.;
1376 ptbins[11]=20.;
1377 ptbins[12]=24.;
1378 ptbins[13]=9999.;
1379
1380 SetPtBins(nptbins+1,ptbins);
1381
fc051756 1382
1383 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*/
1384 {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*/
1385 {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 */
1386 {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 */
1387 {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 */
1388 {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 */
1389 {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 */
1390 {0.400,0.09,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* 6<pt<8 */
1391 {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 */
1392 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 12<pt<16 */
1393 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 16<pt<20 */
1394 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 20<pt<24 */
1395 {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.}};/* pt>24 */
1396
1397 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1398 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1399 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1400 for (Int_t ibin=0;ibin<nptbins;ibin++){
1401 for (Int_t ivar = 0; ivar<nvars; ivar++){
1402 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1403 }
1404 }
1405 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
d1029f2b 1406 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1407 delete [] cutsMatrixTransposeStand;
fc051756 1408
1409
1410 //pid settings
1411 AliAODPidHF* pidObj=new AliAODPidHF();
1412 Int_t mode=1;
1413 const Int_t nlims=2;
1414 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1415 Bool_t compat=kTRUE; //effective only for this mode
1416 Bool_t asym=kTRUE;
1417 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1418 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1419 pidObj->SetMatch(mode);
1420 pidObj->SetPLimit(plims,nlims);
1421 pidObj->SetSigma(sigmas);
1422 pidObj->SetCompat(compat);
1423 pidObj->SetTPC(kTRUE);
1424 pidObj->SetTOF(kTRUE);
ec0f397b 1425 pidObj->SetOldPid(kTRUE);
fc051756 1426 SetPidHF(pidObj);
1427
1428 SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
1429
1430 SetUsePID(kTRUE);
1431
1432 //activate pileup rejection (for pp)
1433 SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1434
1435 //Do recalculate the vertex
1436 SetRemoveDaughtersFromPrim(kTRUE);
1437
1438 // Cut on the zvtx
1439 SetMaxVtxZ(10.);
1440
1441 // Use the kFirst selection for tracks with candidate D with pt<2
1442 SetSelectCandTrackSPDFirst(kTRUE,2.);
1443
1444 // Use special cuts for D candidates with pt<2
1445 SetUseSpecialCuts(kTRUE);
1446 SetMaximumPtSpecialCuts(2.);
1447
1448 PrintAll();
1449
1450 delete pidObj;
1451 pidObj=NULL;
1452
1453 return;
1454}
1455
1456
e2aa82b6 1457//---------------------------------------------------------------------------
fc051756 1458void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
1459 //
1460 // Final CUTS USED FOR 2010 PbPb analysis of central events
1461 // REMEMBER TO EVENTUALLY SWITCH ON
1462 // SetUseAOD049(kTRUE);
1463 //
1464
1465 SetName("D0toKpiCutsStandard");
1466 SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
1467
1468
1469 // CENTRALITY SELECTION
1470 SetMinCentrality(0.);
1471 SetMaxCentrality(80.);
1472 SetUseCentrality(AliRDHFCuts::kCentV0M);
1473
1474
1475
1476 // EVENT CUTS
1477 SetMinVtxContr(1);
1478 // MAX Z-VERTEX CUT
1479 SetMaxVtxZ(10.);
1480 // PILE UP
1481 SetOptPileup(AliRDHFCuts::kNoPileupSelection);
1482 // PILE UP REJECTION
1483 //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1484
1485 // TRACKS ON SINGLE TRACKS
1486 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1487 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1488 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1489 esdTrackCuts->SetRequireITSRefit(kTRUE);
1490 // esdTrackCuts->SetMinNClustersITS(4);
1491 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1492 esdTrackCuts->SetMinDCAToVertexXY(0.);
1493 esdTrackCuts->SetEtaRange(-0.8,0.8);
1494 esdTrackCuts->SetPtRange(0.7,1.e10);
1495
1496 esdTrackCuts->SetMaxDCAToVertexXY(1.);
1497 esdTrackCuts->SetMaxDCAToVertexZ(1.);
1498 esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1499
1500
1501 AddTrackCuts(esdTrackCuts);
1502 // SPD k FIRST for Pt<3 GeV/c
1503 SetSelectCandTrackSPDFirst(kTRUE, 3);
1504
1505 // CANDIDATE CUTS
1506 const Int_t nptbins =13;
1507 const Double_t ptmax = 9999.;
1508 const Int_t nvars=11;
1509 Float_t ptbins[nptbins+1];
1510 ptbins[0]=0.;
1511 ptbins[1]=0.5;
1512 ptbins[2]=1.;
1513 ptbins[3]=2.;
1514 ptbins[4]=3.;
1515 ptbins[5]=4.;
1516 ptbins[6]=5.;
1517 ptbins[7]=6.;
1518 ptbins[8]=8.;
1519 ptbins[9]=12.;
1520 ptbins[10]=16.;
1521 ptbins[11]=20.;
1522 ptbins[12]=24.;
1523 ptbins[13]=ptmax;
1524
1525 SetGlobalIndex(nvars,nptbins);
1526 SetPtBins(nptbins+1,ptbins);
1527 SetMinPtCandidate(2.);
1528
1529 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*/
1530 {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 1531 {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 1532 {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 */
1533 {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 */
1534 {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 */
1535 {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 */
1536 {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 */
1537 {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 */
1538 {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 */
1539 {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 */
1540 {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 */
1541 {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 */
1542
1543
1544 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1545 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1546 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1547
1548 for (Int_t ibin=0;ibin<nptbins;ibin++){
1549 for (Int_t ivar = 0; ivar<nvars; ivar++){
1550 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1551 }
1552 }
1553
1554 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1555 SetUseSpecialCuts(kTRUE);
1556 SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1557 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1558 delete [] cutsMatrixTransposeStand;
1559 cutsMatrixTransposeStand=NULL;
1560
1561 // PID SETTINGS
1562 AliAODPidHF* pidObj=new AliAODPidHF();
1563 //pidObj->SetName("pid4D0");
1564 Int_t mode=1;
1565 const Int_t nlims=2;
1566 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1567 Bool_t compat=kTRUE; //effective only for this mode
1568 Bool_t asym=kTRUE;
1569 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1570 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1571 pidObj->SetMatch(mode);
1572 pidObj->SetPLimit(plims,nlims);
1573 pidObj->SetSigma(sigmas);
1574 pidObj->SetCompat(compat);
1575 pidObj->SetTPC(kTRUE);
1576 pidObj->SetTOF(kTRUE);
1577 pidObj->SetPCompatTOF(2.);
1578 pidObj->SetSigmaForTPCCompat(3.);
1579 pidObj->SetSigmaForTOFCompat(3.);
ec0f397b 1580 pidObj->SetOldPid(kTRUE);
fc051756 1581
1582
1583 SetPidHF(pidObj);
1584 SetUsePID(kTRUE);
1585 SetUseDefaultPID(kFALSE);
1586
1587
1588 PrintAll();
1589
1590
1591 delete pidObj;
1592 pidObj=NULL;
1593
1594 return;
1595
1596}
1597
e2aa82b6 1598//---------------------------------------------------------------------------
fc051756 1599void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() {
1600 // CUTS USED FOR D0 RAA for the CENTRALITY RANGE 20-80%
1601
1602
1603 SetName("D0toKpiCutsStandard");
1604 SetTitle("Standard Cuts for D0 analysis in PbPb2010 run, for peripheral events");
1605
1606
1607 // CENTRALITY SELECTION
1608 SetMinCentrality(40.);
1609 SetMaxCentrality(80.);
1610 SetUseCentrality(AliRDHFCuts::kCentV0M);
1611
1612 // EVENT CUTS
1613 SetMinVtxContr(1);
1614 // MAX Z-VERTEX CUT
1615 SetMaxVtxZ(10.);
1616 // PILE UP
1617 SetOptPileup(AliRDHFCuts::kNoPileupSelection);
1618 // PILE UP REJECTION
1619 //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1620
1621 // TRACKS ON SINGLE TRACKS
1622 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1623 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1624 esdTrackCuts->SetRequireTPCRefit(kTRUE);
1625 esdTrackCuts->SetRequireITSRefit(kTRUE);
1626 // esdTrackCuts->SetMinNClustersITS(4);
1627 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1628 esdTrackCuts->SetMinDCAToVertexXY(0.);
1629 esdTrackCuts->SetEtaRange(-0.8,0.8);
1630 esdTrackCuts->SetPtRange(0.5,1.e10);
1631
1632 esdTrackCuts->SetMaxDCAToVertexXY(1.);
1633 esdTrackCuts->SetMaxDCAToVertexZ(1.);
1634 esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0025*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1635
1636
1637 AddTrackCuts(esdTrackCuts);
1638 // SPD k FIRST for Pt<3 GeV/c
1639 SetSelectCandTrackSPDFirst(kTRUE, 3);
1640
1641 // CANDIDATE CUTS
1642 const Int_t nptbins =13;
1643 const Double_t ptmax = 9999.;
1644 const Int_t nvars=11;
1645 Float_t ptbins[nptbins+1];
1646 ptbins[0]=0.;
1647 ptbins[1]=0.5;
1648 ptbins[2]=1.;
1649 ptbins[3]=2.;
1650 ptbins[4]=3.;
1651 ptbins[5]=4.;
1652 ptbins[6]=5.;
1653 ptbins[7]=6.;
1654 ptbins[8]=8.;
1655 ptbins[9]=12.;
1656 ptbins[10]=16.;
1657 ptbins[11]=20.;
1658 ptbins[12]=24.;
1659 ptbins[13]=ptmax;
1660
1661 SetGlobalIndex(nvars,nptbins);
1662 SetPtBins(nptbins+1,ptbins);
1663 SetMinPtCandidate(0.);
1664
1665 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*/
1666 {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*/
1667 {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 */
1668 {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 */
1669 {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 */
1670 {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 */
1671 {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 */
1672 {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 */
1673 {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 */
1674 {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 */
1675 {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 */
1676 {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 */
1677 {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 */
1678
1679
1680 //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1681 Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1682 for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1683
1684 for (Int_t ibin=0;ibin<nptbins;ibin++){
1685 for (Int_t ivar = 0; ivar<nvars; ivar++){
1686 cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1687 }
1688 }
1689
1690 SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1691 SetUseSpecialCuts(kTRUE);
1692 SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1693 for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1694 delete [] cutsMatrixTransposeStand;
1695 cutsMatrixTransposeStand=NULL;
1696
1697 // PID SETTINGS
1698 AliAODPidHF* pidObj=new AliAODPidHF();
1699 //pidObj->SetName("pid4D0");
1700 Int_t mode=1;
1701 const Int_t nlims=2;
1702 Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1703 Bool_t compat=kTRUE; //effective only for this mode
1704 Bool_t asym=kTRUE;
1705 Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1706 pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1707 pidObj->SetMatch(mode);
1708 pidObj->SetPLimit(plims,nlims);
1709 pidObj->SetSigma(sigmas);
1710 pidObj->SetCompat(compat);
1711 pidObj->SetTPC(kTRUE);
1712 pidObj->SetTOF(kTRUE);
1713 pidObj->SetPCompatTOF(2.);
1714 pidObj->SetSigmaForTPCCompat(3.);
1715 pidObj->SetSigmaForTOFCompat(3.);
ec0f397b 1716 pidObj->SetOldPid(kTRUE);
fc051756 1717
1718 SetPidHF(pidObj);
1719 SetUsePID(kTRUE);
1720 SetUseDefaultPID(kFALSE);
1721
1722 SetLowPt(kTRUE,2.);
1723
1724 PrintAll();
1725
1726
1727 delete pidObj;
1728 pidObj=NULL;
1729
1730 return;
1731
1732}
1733
1734
e2aa82b6 1735//---------------------------------------------------------------------------
fc051756 1736void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() {
1737
1738 // Default 2010 PbPb cut object
1739 SetStandardCutsPbPb2010();
0ebb13ca 1740 AliAODPidHF *pidobj=GetPidHF();
1741
1742 pidobj->SetOldPid(kFALSE);
fc051756 1743
1744 //
1745 // Enable all 2011 PbPb run triggers
1746 //
1747 SetTriggerClass("");
1748 ResetMaskAndEnableMBTrigger();
1749 EnableCentralTrigger();
1750 EnableSemiCentralTrigger();
1751
1752}
e2aa82b6 1753
1754//---------------------------------------------------------------------------
1755Int_t AliRDHFCutsD0toKpi::IsSelectedCombPID(AliAODRecoDecayHF* d)
1756{
1757 // ############################################################
1758 //
1759 // Apply Bayesian PID selection
1760 // Implementation of Bayesian PID: Jeremy Wilkinson
1761 //
1762 // ############################################################
1763
1764 if (!fUsePID || !d) return 3;
1765
775d1f1c 1766
1767 if (fBayesianStrategy == kBayesWeightNoFilter) {
1768 //WeightNoFilter: Accept all particles (no PID cut) but fill mass histos with weights in task
1769 CalculateBayesianWeights(d);
1770 return 3;
6bb29c8c 1771 }
e2aa82b6 1772
1773
1774
1775 Int_t isD0 = 0;
1776 Int_t isD0bar = 0;
1777 Int_t returnvalue = 0;
1778
6bb29c8c 1779 Int_t isPosKaon = 0, isNegKaon = 0, isPosPion = 0, isNegPion = 0;
e2aa82b6 1780
775d1f1c 1781 //Bayesian methods used here check for ID of kaon, and whether it is positive or negative.
1782
1783
e2aa82b6 1784 Bool_t checkPIDInfo[2] = {kTRUE, kTRUE};
1785 AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
1786
1787 if ((aodtrack[0]->Charge()*aodtrack[1]->Charge()) != -1) {
1788 return 0; //reject if charges not opposite
1789 }
b8f433ab 1790 Double_t momentumpositive=0., momentumnegative=0.; //Used in "prob > prior" method
e2aa82b6 1791 for (Int_t daught = 0; daught < 2; daught++) {
1792 //Loop over prongs
1793
b8f433ab 1794 if (aodtrack[daught]->Charge() == -1) {
1795 momentumnegative = aodtrack[daught]->P();
1796 }
1797 if (aodtrack[daught]->Charge() == +1) {
1798 momentumpositive = aodtrack[daught]->P();
1799 }
1800 if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
1801 checkPIDInfo[daught] = kFALSE;
1802 continue;
1803 }
e2aa82b6 1804
b8f433ab 1805 }
e2aa82b6 1806
b8f433ab 1807 //Loop over daughters ends here
e2aa82b6 1808
b8f433ab 1809 if (!checkPIDInfo[0] && !checkPIDInfo[1]) {
1810 return 0; //Reject if both daughters lack both TPC and TOF info
1811 }
e2aa82b6 1812
e2aa82b6 1813
b8f433ab 1814 CalculateBayesianWeights(d); //Calculates all Bayesian probabilities for both positive and negative tracks
1815 // Double_t prob0[AliPID::kSPECIES];
1816 // fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), prob0);
775d1f1c 1817 ///Three possible Bayesian probability cuts: Picked using SetBayesianCondition(int).
b8f433ab 1818 switch (fBayesianCondition) {
775d1f1c 1819 ///A: Standard max. probability method (accept most likely species)
b8f433ab 1820 case kMaxProb:
6bb29c8c 1821 if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kKaon]) { //If highest probability lies with kaon
1822 isPosKaon = 1; //flag [daught] as a kaon
1823 }
1824
1825 if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kPion]) { //If highest probability lies with pion
1826 isPosPion = 1; //flag [daught] as a pion
1827 }
1828
1829 if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kKaon]) { //If highest probability lies with kaon
1830 isNegKaon = 1; //flag [daught] as a kaon
1831 }
1832
1833 if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kPion]) { //If highest probability lies with kaon
1834 isNegPion = 1; //flag [daught] as a pion
1835 }
1836
1837
1838 break;
1839 ///B: Accept if probability greater than prior
1840 case kAbovePrior:
1841
1842 if (fWeightsNegative[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
1843 GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumnegative)))) { //Retrieves relevant prior, gets value at momentum
1844 isNegKaon = 1;
1845 }
1846 if (fWeightsPositive[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
1847 GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumpositive)))) { //Retrieves relevant prior, gets value at momentum
1848 isPosKaon = 1;
1849 }
1850
1851 break;
1852
1853 ///C: Accept if probability greater than user-defined threshold
b8f433ab 1854 case kThreshold:
1855 if (fWeightsNegative[AliPID::kKaon] > fProbThreshold) {
1856 isNegKaon = 1;
1857 }
6bb29c8c 1858 if (fWeightsNegative[AliPID::kPion] > fProbThreshold) {
1859 isNegPion = 1;
1860 }
1861
b8f433ab 1862 if (fWeightsPositive[AliPID::kKaon] > fProbThreshold) {
1863 isPosKaon = 1;
1864 }
6bb29c8c 1865
1866 if (fWeightsPositive[AliPID::kPion] > fProbThreshold) {
1867 isPosPion = 1;
1868 }
b8f433ab 1869
1870 break;
6bb29c8c 1871 }
e2aa82b6 1872
775d1f1c 1873
1874 //Momentum-based selection (also applied to filtered weighted method)
1875
1876 if (fBayesianStrategy == kBayesMomentum || fBayesianCondition == kBayesWeight) {
1877 if (isNegKaon && isPosKaon) { // If both are kaons, reject
775d1f1c 1878 isD0 = 1;
6bb29c8c 1879 isD0bar = 1;
1880 } else if (isNegKaon && isPosPion) { //If negative kaon present, D0
1881 isD0 = 1;
1882 } else if (isPosKaon && isNegPion) { //If positive kaon present, D0bar
775d1f1c 1883 isD0bar = 1;
1884 } else { //If neither ID'd as kaon, subject to extra tests
1885 isD0 = 1;
1886 isD0bar = 1;
1887 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
1888 if (aodtrack[0]->Charge() == -1) {
1889 isD0 = 0; //Check charge and reject based on pion hypothesis
1890 }
1891 if (aodtrack[0]->Charge() == 1) {
1892 isD0bar = 0;
1893 }
1894 }
1895 if (aodtrack[1]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[1], "TOF"))) {
1896 if (aodtrack[1]->Charge() == -1) {
1897 isD0 = 0;
1898 }
1899 if (aodtrack[1]->Charge() == 1) {
1900 isD0bar = 0;
1901 }
1902 }
1903 }
e2aa82b6 1904
775d1f1c 1905
e2aa82b6 1906
775d1f1c 1907 if (isD0 && isD0bar) {
1908 returnvalue = 3;
b8f433ab 1909 }
775d1f1c 1910 if (isD0 && !isD0bar) {
1911 returnvalue = 1;
b8f433ab 1912 }
775d1f1c 1913 if (!isD0 && isD0bar) {
1914 returnvalue = 2;
b8f433ab 1915 }
775d1f1c 1916 if (!isD0 && !isD0bar) {
1917 returnvalue = 0;
b8f433ab 1918 }
775d1f1c 1919 }
b8f433ab 1920
775d1f1c 1921 //Simple Bayesian
1922 if (fBayesianStrategy == kBayesSimple) {
1923
731b87ba 1924 if (isPosKaon && isNegKaon) { //If both are ID'd as kaons, accept as possible
1925 returnvalue = 3;
6bb29c8c 1926 } else if (isNegKaon && isPosPion) { //If negative kaon, D0
775d1f1c 1927 returnvalue = 1;
6bb29c8c 1928 } else if (isPosKaon && isNegPion) { //If positive kaon, D0-bar
775d1f1c 1929 returnvalue = 2;
6bb29c8c 1930 } else if (isPosPion && isNegPion) {
775d1f1c 1931 returnvalue = 0; //If neither kaon, reject
6bb29c8c 1932 } else {returnvalue = 0;} //default
1933
775d1f1c 1934 }
1935
1936 return returnvalue;
e2aa82b6 1937
1938
e2aa82b6 1939
e2aa82b6 1940}
1941
1942
1943
e2aa82b6 1944//---------------------------------------------------------------------------
1945void AliRDHFCutsD0toKpi::CalculateBayesianWeights(AliAODRecoDecayHF* d)
1946
1947{
1948 //Function to compute weights for Bayesian method
e2aa82b6 1949
1950 AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
1951 if ((aodtrack[0]->Charge() * aodtrack[1]->Charge()) != -1) {
1952 return; //Reject if charges do not oppose
1953 }
1954 for (Int_t daught = 0; daught < 2; daught++) {
1955 //Loop over prongs
1956
1957 if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
1958 continue;
1959 }
1960
1961
c2ec6218 1962 // identify kaon, define weights
e2aa82b6 1963 if (aodtrack[daught]->Charge() == +1) {
c2ec6218 1964 fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsPositive);
e2aa82b6 1965 }
c2ec6218 1966
e2aa82b6 1967 if (aodtrack[daught]->Charge() == -1) {
c2ec6218 1968 fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsNegative);
e2aa82b6 1969 }
e2aa82b6 1970 }
1971}
1972
1973
1974