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