ANALYSIS/AliD0toKpi.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //----------------------------------------------------------------------------
  17 //               Implementation of the D0toKpi class
  18 //
  19 // Note: the two decay tracks are labelled: 0 (positive track)
  20 //                                          1 (negative track)
  21 //
  22 //            Origin: A. Dainese    andrea.dainese@pd.infn.it
  23 //----------------------------------------------------------------------------
  24 #include <Riostream.h>
  25 #include <TH1.h>
  26 #include <TH2.h>
  27 #include <TCanvas.h>
  28 #include <TPaveLabel.h>
  29 #include <TVector3.h>
  30
  31 #include "AliD0toKpi.h"
  32
  33 ClassImp(AliD0toKpi)
  34
  35 //----------------------------------------------------------------------------
  36 AliD0toKpi::AliD0toKpi() {
  37   // Default constructor
  38
  39   fSignal = kFALSE;
  40
  41   fEvent = 0;
  42
  43   fTrkNum[0] = 0;
  44   fTrkNum[1] = 0;
  45
  46   fV1x = 0.;
  47   fV1y = 0.;
  48   fV1z = 0.;
  49   fV2x = 0.;
  50   fV2y = 0.;
  51   fV2z = 0.;
  52   fDCA = 0.;
  53
  54   fPx[0] = 0.;
  55   fPy[0] = 0.;
  56   fPz[0] = 0.;
  57   fPx[1] = 0.;
  58   fPy[1] = 0.;
  59   fPz[1] = 0.;
  60
  61   fd0[0] = 0.;
  62   fd0[1] = 0.;
  63
  64   fPdg[0] = 0;
  65   fPdg[1] = 0;
  66   fMum[0] = 0;
  67   fMum[1] = 0;
  68
  69   fTagPi[0] = 0.;
  70   fTagPi[1] = 0.;
  71   fTagKa[0] = 0.;
  72   fTagKa[1] = 0.;
  73   fTagNid[0] = 0.;
  74   fTagNid[1] = 0.;
  75
  76   fWgtAD0=fWgtAD0bar=fWgtBD0=fWgtBD0bar=fWgtCD0=fWgtCD0bar=fWgtDD0=fWgtDD0bar=0;
  77
  78 }
  79 //----------------------------------------------------------------------------
  80 AliD0toKpi::AliD0toKpi(Int_t ev,Int_t trkNum[2],
  81                        Double_t v1[3],Double_t v2[3],
  82                        Double_t dca,
  83                        Double_t mom[6],Double_t d0[2]) {
  84   // Constructor
  85
  86   fSignal = kFALSE;
  87
  88   fEvent = ev;
  89   fTrkNum[0] = trkNum[0];
  90   fTrkNum[1] = trkNum[1];
  91
  92   fV1x = v1[0];
  93   fV1y = v1[1];
  94   fV1z = v1[2];
  95   fV2x = v2[0];
  96   fV2y = v2[1];
  97   fV2z = v2[2];
  98   fDCA = dca;
  99
 100   fPx[0] = mom[0];
 101   fPy[0] = mom[1];
 102   fPz[0] = mom[2];
 103   fPx[1] = mom[3];
 104   fPy[1] = mom[4];
 105   fPz[1] = mom[5];
 106
 107   fd0[0] = d0[0];
 108   fd0[1] = d0[1];
 109
 110   fPdg[0] = 0;
 111   fPdg[1] = 0;
 112   fMum[0] = 0;
 113   fMum[1] = 0;
 114
 115   fTagPi[0]  = 0.;
 116   fTagPi[1]  = 0.;
 117   fTagKa[0]  = 0.;
 118   fTagKa[1]  = 0.;
 119   fTagNid[0] = 0.;
 120   fTagNid[1] = 0.;
 121
 122   fWgtAD0=fWgtAD0bar=fWgtBD0=fWgtBD0bar=fWgtCD0=fWgtCD0bar=fWgtDD0=fWgtDD0bar=0;
 123 }
 124 //----------------------------------------------------------------------------
 125 AliD0toKpi::~AliD0toKpi() {}
 126 //____________________________________________________________________________
 127 AliD0toKpi::AliD0toKpi( const AliD0toKpi& d0toKpi):TObject(d0toKpi) {
 128   // dummy copy constructor
 129 }
 130 //----------------------------------------------------------------------------
 131 void AliD0toKpi::ApplyPID(TString pidScheme) {
 132
 133   const char *tofparampbpb = strstr(pidScheme.Data(),"TOFparam_PbPb");
 134   const char *tofparampp   = strstr(pidScheme.Data(),"TOFparam_pp");
 135
 136   if((tofparampbpb || tofparampp) && fPdg[0]==0) {
 137     printf("AliD0toKpi::ApplyPID :\n Warning: TOF parameterized PID can be used only for simulation!\n");
 138     return;
 139   }
 140
 141   if(tofparampbpb) {
 142     // tagging of the positive track
 143     if(TMath::Abs(fPdg[0])==211 || TMath::Abs(fPdg[0])==13
 144        || TMath::Abs(fPdg[0])==11) { // pion,muon,electron
 145       fTagPi[0]  = LinearInterpolation(PChild(0),kPiBins_PbPb,kPiBinWidth_PbPb,kPiTagPi_PbPb);
 146       fTagNid[0] = 1.-fTagPi[0];
 147       fTagKa[0]   = 0.;
 148       fTagPr[0]   = 0.;
 149     }
 150     if(TMath::Abs(fPdg[0])==321) { // kaon
 151       fTagKa[0]   = LinearInterpolation(PChild(0),kKBins_PbPb,kKBinWidth_PbPb,kKTagK_PbPb);
 152       fTagNid[0] = LinearInterpolation(PChild(0),kKBins_PbPb,kKBinWidth_PbPb,kKTagNid_PbPb);
 153       if((fTagNid[0]+fTagKa[0])>1.) fTagNid[0] = 1.-fTagKa[0];
 154       fTagPi[0] = 1.-fTagNid[0]-fTagKa[0];
 155       fTagPr[0] = 0.;
 156     }
 157     if(TMath::Abs(fPdg[0])==2212) { // proton
 158       fTagPr[0]  = LinearInterpolation(PChild(0),kPBins_PbPb,kPBinWidth_PbPb,kPTagP_PbPb);
 159       fTagNid[0] = LinearInterpolation(PChild(0),kPBins_PbPb,kPBinWidth_PbPb,kPTagNid_PbPb);
 160       if((fTagNid[0]+fTagPr[0])>1.) fTagNid[0] = 1.-fTagPr[0];
 161       fTagPi[0] = 1.-fTagNid[0]-fTagPr[0];
 162       fTagKa[0]   = 0.;
 163     }
 164     // tagging of the negative track
 165     if(TMath::Abs(fPdg[1])==211 || TMath::Abs(fPdg[1])==13
 166        || TMath::Abs(fPdg[1])==11) { // pion,muon,electron
 167       fTagPi[1]  = LinearInterpolation(PChild(1),kPiBins_PbPb,kPiBinWidth_PbPb,kPiTagPi_PbPb);
 168       fTagNid[1] = 1.-fTagPi[1];
 169       fTagKa[1]   = 0.;
 170       fTagPr[1]   = 0.;
 171     }
 172     if(TMath::Abs(fPdg[1])==321) { // kaon
 173       fTagKa[1]   = LinearInterpolation(PChild(1),kKBins_PbPb,kKBinWidth_PbPb,kKTagK_PbPb);
 174       fTagNid[1] = LinearInterpolation(PChild(1),kKBins_PbPb,kKBinWidth_PbPb,kKTagNid_PbPb);
 175       if((fTagNid[1]+fTagKa[1])>1.) fTagNid[1] = 1.-fTagKa[1];
 176       fTagPi[1] = 1.-fTagNid[1]-fTagKa[1];
 177       fTagPr[1] = 0.;
 178     }
 179     if(TMath::Abs(fPdg[1])==2212) { // proton
 180       fTagPr[1]  = LinearInterpolation(PChild(1),kPBins_PbPb,kPBinWidth_PbPb,kPTagP_PbPb);
 181       fTagNid[1] = LinearInterpolation(PChild(1),kPBins_PbPb,kPBinWidth_PbPb,kPTagNid_PbPb);
 182       if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1];
 183       fTagPi[1] = 1.-fTagNid[1]-fTagPr[1];
 184       fTagKa[1]   = 0.;
 185     }
 186   }
 187
 188
 189   if(tofparampp) {
 190     // tagging of the positive track
 191     if(TMath::Abs(fPdg[0])==211 || TMath::Abs(fPdg[0])==13
 192        || TMath::Abs(fPdg[0])==11) { // pion,muon,electron
 193       fTagPi[0]  = LinearInterpolation(PChild(0),kPiBins_pp,kPiBinWidth_pp,kPiTagPi_pp);
 194       fTagNid[0] = 1.-fTagPi[0];
 195       fTagKa[0]   = 0.;
 196       fTagPr[0]   = 0.;
 197     }
 198     if(TMath::Abs(fPdg[0])==321) { // kaon
 199       fTagKa[0]   = LinearInterpolation(PChild(0),kKBins_pp,kKBinWidth_pp,kKTagK_pp);
 200       fTagNid[0] = LinearInterpolation(PChild(0),kKBins_pp,kKBinWidth_pp,kKTagNid_pp);
 201       if((fTagNid[0]+fTagKa[0])>1.) fTagNid[0] = 1.-fTagKa[0];
 202       fTagPi[0] = 1.-fTagNid[0]-fTagKa[0];
 203       fTagPr[0] = 0.;
 204     }
 205     if(TMath::Abs(fPdg[0])==2212) { // proton
 206       fTagPr[0]  = LinearInterpolation(PChild(0),kPBins_pp,kPBinWidth_pp,kPTagP_pp);
 207       fTagNid[0] = LinearInterpolation(PChild(0),kPBins_pp,kPBinWidth_pp,kPTagNid_pp);
 208       if((fTagNid[0]+fTagPr[0])>1.) fTagNid[0] = 1.-fTagPr[0];
 209       fTagPi[0] = 1.-fTagNid[0]-fTagPr[0];
 210       fTagKa[0]   = 0.;
 211     }
 212     // tagging of the negative track
 213     if(TMath::Abs(fPdg[1])==211 || TMath::Abs(fPdg[1])==13
 214        || TMath::Abs(fPdg[1])==11) { // pion,muon,electron
 215       fTagPi[1]  = LinearInterpolation(PChild(1),kPiBins_pp,kPiBinWidth_pp,kPiTagPi_pp);
 216       fTagNid[1] = 1.-fTagPi[1];
 217       fTagKa[1]   = 0.;
 218       fTagPr[1]   = 0.;
 219     }
 220     if(TMath::Abs(fPdg[1])==321) { // kaon
 221       fTagKa[1]   = LinearInterpolation(PChild(1),kKBins_pp,kKBinWidth_pp,kKTagK_pp);
 222       fTagNid[1] = LinearInterpolation(PChild(1),kKBins_pp,kKBinWidth_pp,kKTagNid_pp);
 223       if((fTagNid[1]+fTagKa[1])>1.) fTagNid[1] = 1.-fTagKa[1];
 224       fTagPi[1] = 1.-fTagNid[1]-fTagKa[1];
 225       fTagPr[1] = 0.;
 226     }
 227     if(TMath::Abs(fPdg[1])==2212) { // proton
 228       fTagPr[1]  = LinearInterpolation(PChild(1),kPBins_pp,kPBinWidth_pp,kPTagP_pp);
 229       fTagNid[1] = LinearInterpolation(PChild(1),kPBins_pp,kPBinWidth_pp,kPTagNid_pp);
 230       if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1];
 231       fTagPi[1] = 1.-fTagNid[1]-fTagPr[1];
 232       fTagKa[1]   = 0.;
 233     }
 234   }
 235
 236   return;
 237 }
 238 //----------------------------------------------------------------------------
 239 Double_t AliD0toKpi::ChildrenRelAngle() const {
 240   // relative angle between K and pi
 241
 242   TVector3 mom0(fPx[0],fPy[0],fPz[0]);
 243   TVector3 mom1(fPx[1],fPy[1],fPz[1]);
 244
 245   Double_t angle = mom0.Angle(mom1);
 246
 247   return angle;
 248 }
 249 //----------------------------------------------------------------------------
 250 void AliD0toKpi::ComputeWgts() {
 251   // calculate the weights for PID
 252
 253
 254   // assignement of the weights from PID
 255   fWgtAD0    = fTagKa[1]*(fTagPi[0]+fTagNid[0]);
 256   fWgtAD0bar = fTagKa[0]*(fTagPi[1]+fTagNid[1]);
 257   fWgtBD0    = fTagPi[0]*fTagNid[1];
 258   fWgtBD0bar = fTagPi[1]*fTagNid[0];
 259   fWgtCD0    = fTagNid[0]*fTagNid[1];
 260   fWgtCD0bar = fTagNid[0]*fTagNid[1];
 261   fWgtDD0    = 1.-fWgtAD0-fWgtBD0-fWgtCD0;
 262   fWgtDD0bar = 1.-fWgtAD0bar-fWgtBD0bar-fWgtCD0bar;
 263
 264   /*
 265   cerr<<fWgtAD0<<"  "<<fWgtAD0bar<<endl;
 266   cerr<<fWgtBD0<<"  "<<fWgtBD0bar<<endl;
 267   cerr<<fWgtCD0<<"  "<<fWgtCD0bar<<endl;
 268
 269   if(fWgtAD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
 270   if(fWgtAD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
 271   if(fWgtBD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
 272   if(fWgtBD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
 273   if(fWgtCD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
 274   if(fWgtCD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
 275   */
 276
 277   return;
 278 }
 279 //----------------------------------------------------------------------------
 280 void AliD0toKpi::CorrectWgt4BR(Double_t factor) {
 281   // correct weights of background from charm
 282
 283   fWgtAD0    *= factor;
 284   fWgtAD0bar *= factor;
 285   fWgtBD0    *= factor;
 286   fWgtBD0bar *= factor;
 287   fWgtCD0    *= factor;
 288   fWgtCD0bar *= factor;
 289   fWgtDD0    *= factor;
 290   fWgtDD0bar *= factor;
 291
 292   return;
 293 }
 294 //----------------------------------------------------------------------------
 295 Double_t AliD0toKpi::CosPointing() const {
 296   // cosine of pointing angle in space
 297
 298   TVector3 mom(Px(),Py(),Pz());
 299   TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
 300
 301   Double_t pta = mom.Angle(flight);
 302
 303   return TMath::Cos(pta);
 304 }
 305 //----------------------------------------------------------------------------
 306 Double_t AliD0toKpi::CosPointingXY() const {
 307   // cosine of pointing angle in transverse plane
 308
 309   TVector3 momXY(Px(),Py(),0.);
 310   TVector3 flightXY(fV2x-fV1x,fV2y-fV1y,0.);
 311
 312   Double_t ptaXY = momXY.Angle(flightXY);
 313
 314   return TMath::Cos(ptaXY);
 315 }
 316 //----------------------------------------------------------------------------
 317 void AliD0toKpi::CosThetaStar(Double_t &ctsD0,Double_t &ctsD0bar) const {
 318   // cosine of decay angle in the D0 rest frame
 319
 320   Double_t pStar = TMath::Sqrt(TMath::Power(kMD0*kMD0-kMK*kMK-kMPi*kMPi,2.)-4.*kMK*kMK*kMPi*kMPi)/(2.*kMD0);
 321
 322   Double_t beta = P()/Energy();
 323   Double_t gamma = Energy()/kMD0;
 324
 325   ctsD0 = (Ql(1)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
 326   // if(ctsD0 > 1.)  { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0<<"!\n"; }
 327   // if(ctsD0 < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0<<"!\n"; }
 328
 329   ctsD0bar = (Ql(0)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
 330   // if(ctsD0bar > 1.)  { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0bar<<"!\n"; }
 331   // if(ctsD0bar < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0bar<<"!\n";}
 332
 333   return;
 334 }
 335 //----------------------------------------------------------------------------
 336 Double_t AliD0toKpi::Eta() const {
 337   // pseudorapidity of the D0
 338
 339   Double_t theta = TMath::Pi()/2.-TMath::ATan2(Pz(),Pt());
 340   Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
 341   return eta;
 342 }
 343 //----------------------------------------------------------------------------
 344 Double_t AliD0toKpi::EtaChild(Int_t child) const {
 345   // pseudorapidity of the decay tracks
 346
 347   Double_t theta = TMath::Pi()/2.-TMath::ATan2(fPz[child],PtChild(child));
 348   Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
 349   return eta;
 350 }
 351 //----------------------------------------------------------------------------
 352 void AliD0toKpi::GetWgts(Double_t &WgtD0,Double_t &WgtD0bar,TString sample)
 353   const {
 354   // returns the weights for pid
 355
 356   const char *sampleA = strstr(sample.Data(),"A");
 357   const char *sampleB = strstr(sample.Data(),"B");
 358   const char *sampleC = strstr(sample.Data(),"C");
 359   const char *sampleD = strstr(sample.Data(),"D");
 360   const char *sampleABCD = strstr(sample.Data(),"ABCD");
 361   const char *sampleABC = strstr(sample.Data(),"ABC");
 362   const char *sampleBC = strstr(sample.Data(),"BC");
 363
 364   if(sampleA) { WgtD0 = fWgtAD0;  WgtD0bar = fWgtAD0bar; }
 365   if(sampleB) { WgtD0 = fWgtBD0;  WgtD0bar = fWgtBD0bar; }
 366   if(sampleC) { WgtD0 = fWgtCD0;  WgtD0bar = fWgtCD0bar; }
 367   if(sampleD) { WgtD0 = fWgtDD0;  WgtD0bar = fWgtDD0bar; }
 368   if(sampleABCD) {
 369     WgtD0    = fWgtAD0+fWgtBD0+fWgtCD0+fWgtDD0;
 370     WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar+fWgtDD0bar;
 371   }
 372   if(sampleABC) {
 373     WgtD0    = fWgtAD0+fWgtBD0+fWgtCD0;
 374     WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar;
 375   }
 376   if(sampleBC) {
 377     WgtD0    = fWgtBD0+fWgtCD0;
 378     WgtD0bar = fWgtBD0bar+fWgtCD0bar;
 379   }
 380
 381
 382   if(fSignal) {
 383     if(fMum[0]==421)  WgtD0bar = 0.;
 384     if(fMum[0]==-421) WgtD0 = 0.;
 385   }
 386
 387   return;
 388 }
 389 //----------------------------------------------------------------------------
 390 void AliD0toKpi::InvMass(Double_t &mD0,Double_t &mD0bar) const {
 391   // invariant mass as D0 and as D0bar
 392
 393   Double_t energy[2];
 394
 395   // D0 -> K- Pi+
 396   energy[1] = TMath::Sqrt(kMK*kMK+PChild(1)*PChild(1));
 397   energy[0] = TMath::Sqrt(kMPi*kMPi+PChild(0)*PChild(0));
 398
 399   mD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
 400
 401
 402   // D0bar -> K+ Pi-
 403   energy[0] = TMath::Sqrt(kMK*kMK+PChild(0)*PChild(0));
 404   energy[1] = TMath::Sqrt(kMPi*kMPi+PChild(1)*PChild(1));
 405
 406   mD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
 407
 408   return;
 409
 410 }
 411 //----------------------------------------------------------------------------
 412 Double_t AliD0toKpi::Ql(Int_t child) const {
 413   // longitudinal momentum of decay tracks w.r.t. to D0 momentum
 414
 415   Double_t qL;
 416   TVector3 mom(fPx[child],fPy[child],fPz[child]);
 417   TVector3 momD(Px(),Py(),Pz());
 418
 419   qL = mom.Dot(momD)/momD.Mag();
 420
 421   return qL ;
 422 }
 423 //----------------------------------------------------------------------------
 424 Double_t AliD0toKpi::Qt() const {
 425   // transverse momentum of decay tracks w.r.t. to D0 momentum
 426
 427   TVector3 mom0(fPx[0],fPy[0],fPz[0]);
 428   TVector3 momD(Px(),Py(),Pz());
 429
 430   return mom0.Perp(momD);
 431 }
 432 //----------------------------------------------------------------------------
 433 Bool_t AliD0toKpi::Select(const Double_t* cuts,Int_t& okD0,Int_t& okD0bar)
 434   const {
 435 //
 436 // This function compares the D0 with a set of cuts:
 437 //
 438 // cuts[0] = inv. mass half width [GeV]
 439 // cuts[1] = dca [micron]
 440 // cuts[2] = cosThetaStar
 441 // cuts[3] = pTK [GeV/c]
 442 // cuts[4] = pTPi [GeV/c]
 443 // cuts[5] = d0K [micron]   upper limit!
 444 // cuts[6] = d0Pi [micron]  upper limit!
 445 // cuts[7] = d0d0 [micron^2]
 446 // cuts[8] = cosThetaPoint
 447 //
 448 // If the D0/D0bar doesn't pass the cuts it sets the weights to 0
 449 // If neither D0 nor D0bar pass the cuts return kFALSE
 450 //
 451   Double_t mD0,mD0bar,ctsD0,ctsD0bar;
 452   okD0=1; okD0bar=1;
 453
 454   if(PtChild(1) < cuts[3] || PtChild(0) < cuts[4]) okD0 = 0;
 455   if(PtChild(0) < cuts[3] || PtChild(1) < cuts[4]) okD0bar = 0;
 456   if(!okD0 && !okD0bar) return kFALSE;
 457
 458   if(TMath::Abs(Getd0Child(1)) > cuts[5] ||
 459      TMath::Abs(Getd0Child(0)) > cuts[6]) okD0 = 0;
 460   if(TMath::Abs(Getd0Child(0)) > cuts[6] ||
 461      TMath::Abs(Getd0Child(1)) > cuts[5]) okD0bar = 0;
 462   if(!okD0 && !okD0bar) return kFALSE;
 463
 464   if(GetDCA() > cuts[1]) { okD0 = okD0bar = 0; return kFALSE; }
 465
 466   InvMass(mD0,mD0bar);
 467   if(TMath::Abs(mD0-kMD0)    > cuts[0]) okD0 = 0;
 468   if(TMath::Abs(mD0bar-kMD0) > cuts[0]) okD0bar = 0;
 469   if(!okD0 && !okD0bar) return kFALSE;
 470
 471   CosThetaStar(ctsD0,ctsD0bar);
 472   if(TMath::Abs(ctsD0)    > cuts[2]) okD0 = 0;
 473   if(TMath::Abs(ctsD0bar) > cuts[2]) okD0bar = 0;
 474   if(!okD0 && !okD0bar) return kFALSE;
 475
 476   if(ProdImpParams() > cuts[7]) { okD0 = okD0bar = 0; return kFALSE; }
 477
 478   if(CosPointing()   < cuts[8]) { okD0 = okD0bar = 0; return kFALSE; }
 479
 480   return kTRUE;
 481 }
 482 //-----------------------------------------------------------------------------
 483 void AliD0toKpi::SetPIDresponse(Double_t resp0[5],Double_t resp1[5]) {
 484   // Set combined PID detector response probabilities
 485
 486   fPIDrespEl[0] = resp0[0];
 487   fPIDrespEl[1] = resp1[0];
 488   fPIDrespMu[0] = resp0[1];
 489   fPIDrespMu[1] = resp1[1];
 490   fPIDrespPi[0] = resp0[2];
 491   fPIDrespPi[1] = resp1[2];
 492   fPIDrespKa[0] = resp0[3];
 493   fPIDrespKa[1] = resp1[3];
 494   fPIDrespPr[0] = resp0[4];
 495   fPIDrespPr[1] = resp1[4];
 496
 497   return;
 498 }
 499 //-----------------------------------------------------------------------------
 500 void AliD0toKpi::DrawPIDinTOF(TString pidScheme) const {
 501   // Draw parameterized PID probabilities in TOF
 502
 503   const char *tofparampbpb = strstr(pidScheme.Data(),"TOFparam_PbPb");
 504   const char *tofparampp = strstr(pidScheme.Data(),"TOFparam_pp");
 505
 506   TH2F* framePi = new TH2F("framePi","Tag probabilities for PIONS",2,0,2.5,2,0,1);
 507   framePi->SetXTitle("p [GeV/c]");
 508   framePi->SetStats(0);
 509   TH2F* frameK = new TH2F("frameK","Tag probabilities for KAONS",2,0,2.5,2,0,1);
 510   frameK->SetXTitle("p [GeV/c]");
 511   frameK->SetStats(0);
 512   TH2F* frameP = new TH2F("frameP","Tag probabilities for PROTONS",2,0,4.5,2,0,1);
 513   frameP->SetXTitle("p [GeV/c]");
 514   frameP->SetStats(0);
 515
 516   TH1F* hPiPi = new TH1F("hPiPi","Tag probabilities for PIONS",kPiBins_PbPb,0,2.5);
 517   TH1F* hPiNid = new TH1F("hPiNid","Tag probabilities for PIONS",kPiBins_PbPb,0,2.5);
 518
 519   TH1F* hKK = new TH1F("hKK","Tag probabilities for KAONS",kKBins_PbPb,0,2.5);
 520   TH1F* hKNid = new TH1F("hKNid","Tag probabilities for KAONS",kKBins_PbPb,0,2.5);
 521   TH1F* hKPi = new TH1F("hKPi","Tag probabilities for KAONS",kKBins_PbPb,0,2.5);
 522
 523   TH1F* hPP = new TH1F("hPP","Tag probabilities for PROTONS",kPBins_PbPb,0,4.5);
 524   TH1F* hPNid = new TH1F("hPNid","Tag probabilities for PROTONS",kPBins_PbPb,0,4.5);
 525   TH1F* hPPi = new TH1F("hPPi","Tag probabilities for PROTONS",kPBins_PbPb,0,4.5);
 526
 527
 528   if(tofparampbpb) {
 529
 530     for(Int_t i=1; i<=kPiBins_PbPb; i++) {
 531       hPiPi->SetBinContent(i,kPiTagPi_PbPb[i-1]);
 532       hPiNid->SetBinContent(i,kPiTagPi_PbPb[i-1]+kPiTagNid_PbPb[i-1]);
 533
 534       hKK->SetBinContent(i,kKTagK_PbPb[i-1]);
 535       hKPi->SetBinContent(i,kKTagK_PbPb[i-1]+kKTagPi_PbPb[i-1]);
 536       hKNid->SetBinContent(i,kKTagK_PbPb[i-1]+kKTagPi_PbPb[i-1]+kKTagNid_PbPb[i-1]);
 537     }
 538     for(Int_t i=1; i<=kPBins_PbPb; i++) {
 539       hPP->SetBinContent(i,kPTagP_PbPb[i-1]);
 540       hPPi->SetBinContent(i,kPTagP_PbPb[i-1]+kPTagPi_PbPb[i-1]);
 541       hPNid->SetBinContent(i,kPTagP_PbPb[i-1]+kPTagPi_PbPb[i-1]+kPTagNid_PbPb[i-1]);
 542     }
 543
 544   } else if(tofparampp) {
 545
 546     for(Int_t i=1; i<=kPiBins_pp; i++) {
 547       hPiPi->SetBinContent(i,kPiTagPi_pp[i-1]);
 548       hPiNid->SetBinContent(i,kPiTagPi_pp[i-1]+kPiTagNid_pp[i-1]);
 549
 550       hKK->SetBinContent(i,kKTagK_pp[i-1]);
 551       hKPi->SetBinContent(i,kKTagK_pp[i-1]+kKTagPi_pp[i-1]);
 552       hKNid->SetBinContent(i,kKTagK_pp[i-1]+kKTagPi_pp[i-1]+kKTagNid_pp[i-1]);
 553     }
 554     for(Int_t i=1; i<=kPBins_pp; i++) {
 555       hPP->SetBinContent(i,kPTagP_pp[i-1]);
 556       hPPi->SetBinContent(i,kPTagP_pp[i-1]+kPTagPi_pp[i-1]);
 557       hPNid->SetBinContent(i,kPTagP_pp[i-1]+kPTagPi_pp[i-1]+kPTagNid_pp[i-1]);
 558     }
 559
 560   }
 561
 562
 563   TCanvas* c = new TCanvas("c","Parameterized PID in TOF",0,0,1000,400);
 564   c->Divide(3,1);
 565   c->cd(1);
 566   framePi->Draw();
 567   hPiNid->SetFillColor(18); hPiNid->Draw("same");
 568   hPiPi->SetFillColor(4); hPiPi->Draw("same");
 569   TPaveLabel* pav1 = new TPaveLabel(1,.2,1.4,.3,"#pi");
 570   pav1->SetBorderSize(0);
 571   pav1->Draw("same");
 572   TPaveLabel* pav2 = new TPaveLabel(1,.8,1.8,.9,"non-id");
 573   pav2->SetBorderSize(0);
 574   pav2->Draw("same");
 575
 576   c->cd(2);
 577   frameK->Draw();
 578   hKNid->SetFillColor(18); hKNid->Draw("same");
 579   hKPi->SetFillColor(4); hKPi->Draw("same");
 580   hKK->SetFillColor(7); hKK->Draw("same");
 581   TPaveLabel* pav3 = new TPaveLabel(1,.2,1.5,.3,"K");
 582   pav3->SetBorderSize(0);
 583   pav3->Draw("same");
 584   TPaveLabel* pav4 = new TPaveLabel(1,.8,1.8,.9,"non-id");
 585   pav4->SetBorderSize(0);
 586   pav4->Draw("same");
 587   TPaveLabel* pav5 = new TPaveLabel(.4,.5,.8,.6,"#pi");
 588   pav5->SetBorderSize(0);
 589   pav5->Draw("same");
 590
 591   c->cd(3);
 592   frameP->Draw();
 593   hPNid->SetFillColor(18); hPNid->Draw("same");
 594   hPPi->SetFillColor(4); hPPi->Draw("same");
 595   hPP->SetFillColor(3); hPP->Draw("same");
 596   TPaveLabel* pav6 = new TPaveLabel(1,.2,1.5,.3,"p");
 597   pav6->SetBorderSize(0);
 598   pav6->Draw("same");
 599   TPaveLabel* pav7 = new TPaveLabel(1,.8,2.6,.9,"non-id");
 600   pav7->SetBorderSize(0);
 601   pav7->Draw("same");
 602   TPaveLabel* pav8 = new TPaveLabel(.2,.5,1,.6,"#pi");
 603   pav8->SetBorderSize(0);
 604   pav8->Draw("same");
 605
 606
 607   return;
 608 }
 609 //----------------------------------------------------------------------------
 610 Double_t AliD0toKpi::LinearInterpolation(Double_t p,Int_t nBins,Double_t Bin,
 611                                          const Double_t *values) const {
 612   // a linear interpolation method
 613
 614   Double_t value=0;
 615   Double_t slope;
 616
 617   if(p<0.5*Bin) {
 618     value = values[0];
 619   } else if(p>=(nBins-0.5)*Bin) {
 620     slope = (2*values[nBins-1]-values[nBins-2]-values[nBins-3])/Bin/2;
 621     value = values[nBins-2]+slope*(p-Bin*(nBins-1.5));
 622   } else {
 623     for(Int_t i=0; i<nBins; i++) {
 624       if(p<(i+0.5)*Bin) {
 625         slope = (values[i]-values[i-1])/Bin;
 626         value = values[i-1]+slope*(p-Bin*(i-0.5));
 627         break;
 628       }
 629     }
 630   }
 631
 632   if(value<0.) value=0.;
 633   if(value>1.) value=1.;
 634
 635   return value;
 636 }
 637 //----------------------------------------------------------------------------
 638
 639
 640
 641
 642
 643
 644 /*
 645 //____________________________________________________________________________
 646 void AliD0toKpi::SetPtWgts4pp() {
 647   // Correct pt distribution in order to reproduce MNR pt slope
 648   // (for pp generated with PYTHIA min. bias)
 649
 650   if(TMath::Abs(fMum[0]) != 421 && TMath::Abs(fMum[1]) != 421 &&
 651      TMath::Abs(fMum[0]) != 411 && TMath::Abs(fMum[1]) != 411) return;
 652
 653   Double_t ptWgt = 1.;
 654   ptWgt = 2.05-0.47*Pt()+0.02*Pt()*Pt();
 655   if(Pt() >= 5.) ptWgt = 0.56*TMath::Exp(-0.12*Pt());
 656
 657   fWgtAD0    *= ptWgt;
 658   fWgtAD0bar *= ptWgt;
 659   fWgtBD0    *= ptWgt;
 660   fWgtBD0bar *= ptWgt;
 661   fWgtCD0    *= ptWgt;
 662   fWgtCD0bar *= ptWgt;
 663   fWgtDD0    *= ptWgt;
 664   fWgtDD0bar *= ptWgt;
 665
 666   return;
 667 }
 668 //____________________________________________________________________________
 669 */
 670
 671
 672
 673
 674
 675
 676
 677
 678