[u/mrichter/AliRoot.git] / PWG3 / AliD0toKpi.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//----------------------------------------------------------------------------
//               Implementation of the D0toKpi class
//                  for pp and PbPb interactions
// Note: the two decay tracks are labelled: 0 (positive track)
//                                          1 (negative track)
//            Origin: A. Dainese    andrea.dainese@lnl.infn.it            
//----------------------------------------------------------------------------

#include <TH1.h>
#include <TH2.h>
#include <TCanvas.h>
#include <TPaveLabel.h>
#include <TVector3.h>

#include "AliD0toKpi.h"

ClassImp(AliD0toKpi)

//----------------------------------------------------------------------------
AliD0toKpi::AliD0toKpi():
fSignal(kFALSE),
fEvent(0),
fV1x(0.),
fV1y(0.),
fV1z(0.),
fV2x(0.),
fV2y(0.),
fV2z(0.),
fDCA(0.),
fWgtAD0(0.),
fWgtAD0bar(0.),
fWgtBD0(0.),
fWgtBD0bar(0.),
fWgtCD0(0.),
fWgtCD0bar(0.),
fWgtDD0(0.),
fWgtDD0bar(0.)
{
  // Default constructor

  fTrkNum[0] = 0;
  fTrkNum[1] = 0;

  fPx[0] = 0.;
  fPy[0] = 0.;
  fPz[0] = 0.;
  fPx[1] = 0.;
  fPy[1] = 0.;
  fPz[1] = 0.;

  fd0[0] = 0.;
  fd0[1] = 0.;

  fPdg[0] = 0;
  fPdg[1] = 0;
  fMum[0] = 0;
  fMum[1] = 0;

  fTagPi[0] = 0.;
  fTagPi[1] = 0.;
  fTagKa[0] = 0.;
  fTagKa[1] = 0.;
  fTagNid[0] = 0.;
  fTagNid[1] = 0.;

}
//----------------------------------------------------------------------------
AliD0toKpi::AliD0toKpi(Int_t ev,Int_t trkNum[2],
		       Double_t v1[3],Double_t v2[3], 
		       Double_t dca,
		       Double_t mom[6],Double_t d0[2]): 
fSignal(kFALSE),
fEvent(ev),
fV1x(v1[0]),
fV1y(v1[1]),
fV1z(v1[2]),
fV2x(v2[0]),
fV2y(v2[1]),
fV2z(v2[2]),
fDCA(dca),
fWgtAD0(0.),
fWgtAD0bar(0.),
fWgtBD0(0.),
fWgtBD0bar(0.),
fWgtCD0(0.),
fWgtCD0bar(0.),
fWgtDD0(0.),
fWgtDD0bar(0.)
{
  // Constructor

  fTrkNum[0] = trkNum[0];
  fTrkNum[1] = trkNum[1];

  fPx[0] = mom[0];
  fPy[0] = mom[1];
  fPz[0] = mom[2];
  fPx[1] = mom[3];
  fPy[1] = mom[4];
  fPz[1] = mom[5];

  fd0[0] = d0[0];
  fd0[1] = d0[1];

  fPdg[0] = 0;
  fPdg[1] = 0;
  fMum[0] = 0;
  fMum[1] = 0;

  fTagPi[0]  = 0.;
  fTagPi[1]  = 0.;
  fTagKa[0]  = 0.;
  fTagKa[1]  = 0.;
  fTagNid[0] = 0.;
  fTagNid[1] = 0.;
}
//----------------------------------------------------------------------------
AliD0toKpi::~AliD0toKpi() {}
//----------------------------------------------------------------------------
void AliD0toKpi::ApplyPID(TString pidScheme) {
  // Applies particle identification

  if((!pidScheme.CompareTo("TOFparamPbPb") || !pidScheme.CompareTo("TOFparamPP")) && fPdg[0]==0) {
    printf("AliD0toKpi::ApplyPID :\n Warning: TOF parameterized PID can be used only for simulation!\n"); 
    return;
  }

  if(!pidScheme.CompareTo("TOFparamPbPb")) {
    // tagging of the positive track
    if(TMath::Abs(fPdg[0])==211 || TMath::Abs(fPdg[0])==13 
       || TMath::Abs(fPdg[0])==11) { // pion,muon,electron
      fTagPi[0]  = LinearInterpolation(PChild(0),kPiBinsPbPb,kPiBinWidthPbPb,kPiTagPiPbPb);
      fTagNid[0] = 1.-fTagPi[0];
      fTagKa[0]   = 0.;
      fTagPr[0]   = 0.;
    } 
    if(TMath::Abs(fPdg[0])==321) { // kaon
      fTagKa[0]   = LinearInterpolation(PChild(0),kKBinsPbPb,kKBinWidthPbPb,kKTagKPbPb);
      fTagNid[0] = LinearInterpolation(PChild(0),kKBinsPbPb,kKBinWidthPbPb,kKTagNidPbPb);
      if((fTagNid[0]+fTagKa[0])>1.) fTagNid[0] = 1.-fTagKa[0];
      fTagPi[0] = 1.-fTagNid[0]-fTagKa[0];
      fTagPr[0] = 0.;
    } 
    if(TMath::Abs(fPdg[0])==2212) { // proton
      fTagPr[0]  = LinearInterpolation(PChild(0),kPBinsPbPb,kPBinWidthPbPb,kPTagPPbPb);
      fTagNid[0] = LinearInterpolation(PChild(0),kPBinsPbPb,kPBinWidthPbPb,kPTagNidPbPb);
      if((fTagNid[0]+fTagPr[0])>1.) fTagNid[0] = 1.-fTagPr[0];
      fTagPi[0] = 1.-fTagNid[0]-fTagPr[0];
      fTagKa[0]   = 0.;
    } 
    // tagging of the negative track
    if(TMath::Abs(fPdg[1])==211 || TMath::Abs(fPdg[1])==13 
       || TMath::Abs(fPdg[1])==11) { // pion,muon,electron
      fTagPi[1]  = LinearInterpolation(PChild(1),kPiBinsPbPb,kPiBinWidthPbPb,kPiTagPiPbPb);
      fTagNid[1] = 1.-fTagPi[1];
      fTagKa[1]   = 0.;
      fTagPr[1]   = 0.;
    } 
    if(TMath::Abs(fPdg[1])==321) { // kaon
      fTagKa[1]   = LinearInterpolation(PChild(1),kKBinsPbPb,kKBinWidthPbPb,kKTagKPbPb);
      fTagNid[1] = LinearInterpolation(PChild(1),kKBinsPbPb,kKBinWidthPbPb,kKTagNidPbPb);
      if((fTagNid[1]+fTagKa[1])>1.) fTagNid[1] = 1.-fTagKa[1];
      fTagPi[1] = 1.-fTagNid[1]-fTagKa[1];
      fTagPr[1] = 0.;
    } 
    if(TMath::Abs(fPdg[1])==2212) { // proton
      fTagPr[1]  = LinearInterpolation(PChild(1),kPBinsPbPb,kPBinWidthPbPb,kPTagPPbPb);
      fTagNid[1] = LinearInterpolation(PChild(1),kPBinsPbPb,kPBinWidthPbPb,kPTagNidPbPb);
      if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1];
      fTagPi[1] = 1.-fTagNid[1]-fTagPr[1];
      fTagKa[1]   = 0.;
    } 
  }


  if(!pidScheme.CompareTo("TOFparamPP")) {
    // tagging of the positive track
    if(TMath::Abs(fPdg[0])==211 || TMath::Abs(fPdg[0])==13 
       || TMath::Abs(fPdg[0])==11) { // pion,muon,electron
      fTagPi[0]  = LinearInterpolation(PChild(0),kPiBinsPP,kPiBinWidthPP,kPiTagPiPP);
      fTagNid[0] = 1.-fTagPi[0];
      fTagKa[0]   = 0.;
      fTagPr[0]   = 0.;
    } 
    if(TMath::Abs(fPdg[0])==321) { // kaon
      fTagKa[0]   = LinearInterpolation(PChild(0),kKBinsPP,kKBinWidthPP,kKTagKPP);
      fTagNid[0] = LinearInterpolation(PChild(0),kKBinsPP,kKBinWidthPP,kKTagNidPP);
      if((fTagNid[0]+fTagKa[0])>1.) fTagNid[0] = 1.-fTagKa[0];
      fTagPi[0] = 1.-fTagNid[0]-fTagKa[0];
      fTagPr[0] = 0.;
    } 
    if(TMath::Abs(fPdg[0])==2212) { // proton
      fTagPr[0]  = LinearInterpolation(PChild(0),kPBinsPP,kPBinWidthPP,kPTagPPP);
      fTagNid[0] = LinearInterpolation(PChild(0),kPBinsPP,kPBinWidthPP,kPTagNidPP);
      if((fTagNid[0]+fTagPr[0])>1.) fTagNid[0] = 1.-fTagPr[0];
      fTagPi[0] = 1.-fTagNid[0]-fTagPr[0];
      fTagKa[0]   = 0.;
    } 
    // tagging of the negative track
    if(TMath::Abs(fPdg[1])==211 || TMath::Abs(fPdg[1])==13 
       || TMath::Abs(fPdg[1])==11) { // pion,muon,electron
      fTagPi[1]  = LinearInterpolation(PChild(1),kPiBinsPP,kPiBinWidthPP,kPiTagPiPP);
      fTagNid[1] = 1.-fTagPi[1];
      fTagKa[1]   = 0.;
      fTagPr[1]   = 0.;
    } 
    if(TMath::Abs(fPdg[1])==321) { // kaon
      fTagKa[1]   = LinearInterpolation(PChild(1),kKBinsPP,kKBinWidthPP,kKTagKPP);
      fTagNid[1] = LinearInterpolation(PChild(1),kKBinsPP,kKBinWidthPP,kKTagNidPP);
      if((fTagNid[1]+fTagKa[1])>1.) fTagNid[1] = 1.-fTagKa[1];
      fTagPi[1] = 1.-fTagNid[1]-fTagKa[1];
      fTagPr[1] = 0.;
    } 
    if(TMath::Abs(fPdg[1])==2212) { // proton
      fTagPr[1]  = LinearInterpolation(PChild(1),kPBinsPP,kPBinWidthPP,kPTagPPP);
      fTagNid[1] = LinearInterpolation(PChild(1),kPBinsPP,kPBinWidthPP,kPTagNidPP);
      if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1];
      fTagPi[1] = 1.-fTagNid[1]-fTagPr[1];
      fTagKa[1]   = 0.;
    } 
  }

  return;
}
//----------------------------------------------------------------------------
Double_t AliD0toKpi::ChildrenRelAngle() const {
  // relative angle between K and pi

  TVector3 mom0(fPx[0],fPy[0],fPz[0]);
  TVector3 mom1(fPx[1],fPy[1],fPz[1]);

  Double_t angle = mom0.Angle(mom1);

  return angle; 
}
//----------------------------------------------------------------------------
void AliD0toKpi::ComputeWgts() {
  // calculate the weights for PID


  // assignement of the weights from PID
  fWgtAD0    = fTagKa[1]*(fTagPi[0]+fTagNid[0]);
  fWgtAD0bar = fTagKa[0]*(fTagPi[1]+fTagNid[1]);
  fWgtBD0    = fTagPi[0]*fTagNid[1];
  fWgtBD0bar = fTagPi[1]*fTagNid[0];
  fWgtCD0    = fTagNid[0]*fTagNid[1];
  fWgtCD0bar = fTagNid[0]*fTagNid[1];
  fWgtDD0    = 1.-fWgtAD0-fWgtBD0-fWgtCD0;
  fWgtDD0bar = 1.-fWgtAD0bar-fWgtBD0bar-fWgtCD0bar;

  /*  
  cerr<<fWgtAD0<<"  "<<fWgtAD0bar<<endl;
  cerr<<fWgtBD0<<"  "<<fWgtBD0bar<<endl;
  cerr<<fWgtCD0<<"  "<<fWgtCD0bar<<endl;
  cerr<<fWgtDD0<<"  "<<fWgtDD0bar<<endl;
  */
  /*
  if(fWgtAD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
  if(fWgtAD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
  if(fWgtBD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
  if(fWgtBD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
  if(fWgtCD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
  if(fWgtCD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
  */

  return;
}
//----------------------------------------------------------------------------
void AliD0toKpi::CorrectWgt4BR(Double_t factor) {
  // correct weights of background from charm 

  fWgtAD0    *= factor;
  fWgtAD0bar *= factor;
  fWgtBD0    *= factor;
  fWgtBD0bar *= factor;
  fWgtCD0    *= factor;
  fWgtCD0bar *= factor;
  fWgtDD0    *= factor;
  fWgtDD0bar *= factor;

  return;
}
//----------------------------------------------------------------------------
Double_t AliD0toKpi::CosPointing() const {
  // cosine of pointing angle in space

  TVector3 mom(Px(),Py(),Pz());
  TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);

  Double_t pta = mom.Angle(flight);

  return TMath::Cos(pta); 
}
//----------------------------------------------------------------------------
Double_t AliD0toKpi::CosPointingXY() const {
  // cosine of pointing angle in transverse plane

  TVector3 momXY(Px(),Py(),0.);
  TVector3 flightXY(fV2x-fV1x,fV2y-fV1y,0.);

  Double_t ptaXY = momXY.Angle(flightXY);

  return TMath::Cos(ptaXY); 
}
//----------------------------------------------------------------------------
void AliD0toKpi::CosThetaStar(Double_t &ctsD0,Double_t &ctsD0bar) const {
  // cosine of decay angle in the D0 rest frame

  Double_t pStar = TMath::Sqrt(TMath::Power(kMD0*kMD0-kMK*kMK-kMPi*kMPi,2.)-4.*kMK*kMK*kMPi*kMPi)/(2.*kMD0);

  Double_t beta = P()/Energy();
  Double_t gamma = Energy()/kMD0;

  ctsD0 = (Ql(1)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
  // if(ctsD0 > 1.)  { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0<<"!\n"; }
  // if(ctsD0 < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0<<"!\n"; }

  ctsD0bar = (Ql(0)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
  // if(ctsD0bar > 1.)  { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0bar<<"!\n"; }
  // if(ctsD0bar < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0bar<<"!\n";} 

  return;
}
//----------------------------------------------------------------------------
Double_t AliD0toKpi::Eta() const {
  // pseudorapidity of the D0

  Double_t theta = TMath::Pi()/2.-TMath::ATan2(Pz(),Pt());
  Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
  return eta;
}
//----------------------------------------------------------------------------
Double_t AliD0toKpi::EtaChild(Int_t child) const {
  // pseudorapidity of the decay tracks

  Double_t theta = TMath::Pi()/2.-TMath::ATan2(fPz[child],PtChild(child));
  Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
  return eta;
}
//----------------------------------------------------------------------------
void AliD0toKpi::GetWgts(Double_t &WgtD0,Double_t &WgtD0bar,TString sample) 
  const {
  // returns the weights for pid

  if(!sample.CompareTo("A")) { WgtD0 = fWgtAD0;  WgtD0bar = fWgtAD0bar; }
  if(!sample.CompareTo("B")) { WgtD0 = fWgtBD0;  WgtD0bar = fWgtBD0bar; }
  if(!sample.CompareTo("C")) { WgtD0 = fWgtCD0;  WgtD0bar = fWgtCD0bar; }
  if(!sample.CompareTo("D")) { WgtD0 = fWgtDD0;  WgtD0bar = fWgtDD0bar; }
  if(!sample.CompareTo("ABCD")) { 
    WgtD0    = fWgtAD0+fWgtBD0+fWgtCD0+fWgtDD0; 
    WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar+fWgtDD0bar; 
  }
  if(!sample.CompareTo("ABC")) { 
    WgtD0    = fWgtAD0+fWgtBD0+fWgtCD0; 
    WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar; 
  }
  if(!sample.CompareTo("BC")) { 
    WgtD0    = fWgtBD0+fWgtCD0; 
    WgtD0bar = fWgtBD0bar+fWgtCD0bar; 
  }

  return;
}
//----------------------------------------------------------------------------
Double_t AliD0toKpi::ImpPar() const {
  // D0 impact parameter in the bending plane
  
    Double_t k = -(fV2x-fV1x)*Px()-(fV2y-fV1y)*Py();
    k /= Pt()*Pt();
    Double_t dx = fV2x-fV1x+k*Px();
    Double_t dy = fV2y-fV1y+k*Py();
    Double_t absDD = TMath::Sqrt(dx*dx+dy*dy);
    TVector3 mom(Px(),Py(),Pz());
    TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
    TVector3 cross = mom.Cross(flight);
    return (cross.Z()>0. ? absDD : -absDD);
}
//----------------------------------------------------------------------------
void AliD0toKpi::InvMass(Double_t &mD0,Double_t &mD0bar) const {
  // invariant mass as D0 and as D0bar

  Double_t energy[2];

  // D0 -> K- Pi+
  energy[1] = TMath::Sqrt(kMK*kMK+PChild(1)*PChild(1));
  energy[0] = TMath::Sqrt(kMPi*kMPi+PChild(0)*PChild(0));

  mD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
    

  // D0bar -> K+ Pi-
  energy[0] = TMath::Sqrt(kMK*kMK+PChild(0)*PChild(0));
  energy[1] = TMath::Sqrt(kMPi*kMPi+PChild(1)*PChild(1));

  mD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
    
  return;

}
//----------------------------------------------------------------------------
Double_t AliD0toKpi::Ql(Int_t child) const {
  // longitudinal momentum of decay tracks w.r.t. to D0 momentum

  Double_t qL;
  TVector3 mom(fPx[child],fPy[child],fPz[child]);
  TVector3 momD(Px(),Py(),Pz());

  qL = mom.Dot(momD)/momD.Mag();

  return qL ;
}
//----------------------------------------------------------------------------
Double_t AliD0toKpi::Qt() const {
  // transverse momentum of decay tracks w.r.t. to D0 momentum  

  TVector3 mom0(fPx[0],fPy[0],fPz[0]);
  TVector3 momD(Px(),Py(),Pz());

  return mom0.Perp(momD);
}
//----------------------------------------------------------------------------
Bool_t AliD0toKpi::Select(const Double_t* cuts,Int_t& okD0,Int_t& okD0bar) 
  const {
//
// This function compares the D0 with a set of cuts:
//
// cuts[0] = inv. mass half width [GeV]   
// cuts[1] = dca [micron]
// cuts[2] = cosThetaStar 
// cuts[3] = pTK [GeV/c]
// cuts[4] = pTPi [GeV/c]
// cuts[5] = d0K [micron]   upper limit!
// cuts[6] = d0Pi [micron]  upper limit!
// cuts[7] = d0d0 [micron^2]
// cuts[8] = cosThetaPoint
//
// If the D0/D0bar doesn't pass the cuts it sets the weights to 0
// If neither D0 nor D0bar pass the cuts return kFALSE
//
  Double_t mD0,mD0bar,ctsD0,ctsD0bar;
  okD0=1; okD0bar=1;

  if(PtChild(1) < cuts[3] || PtChild(0) < cuts[4]) okD0 = 0;
  if(PtChild(0) < cuts[3] || PtChild(1) < cuts[4]) okD0bar = 0;
  if(!okD0 && !okD0bar) return kFALSE;

  if(TMath::Abs(Getd0Child(1)) > cuts[5] || 
     TMath::Abs(Getd0Child(0)) > cuts[6]) okD0 = 0;
  if(TMath::Abs(Getd0Child(0)) > cuts[6] ||
     TMath::Abs(Getd0Child(1)) > cuts[5]) okD0bar = 0;
  if(!okD0 && !okD0bar) return kFALSE;

  if(GetDCA() > cuts[1]) { okD0 = okD0bar = 0; return kFALSE; }

  InvMass(mD0,mD0bar);
  if(TMath::Abs(mD0-kMD0)    > cuts[0]) okD0 = 0;
  if(TMath::Abs(mD0bar-kMD0) > cuts[0]) okD0bar = 0;
  if(!okD0 && !okD0bar) return kFALSE;

  CosThetaStar(ctsD0,ctsD0bar);
  if(TMath::Abs(ctsD0)    > cuts[2]) okD0 = 0;
  if(TMath::Abs(ctsD0bar) > cuts[2]) okD0bar = 0;
  if(!okD0 && !okD0bar) return kFALSE;

  if(ProdImpParams() > cuts[7]) { okD0 = okD0bar = 0; return kFALSE; }

  if(CosPointing()   < cuts[8]) { okD0 = okD0bar = 0; return kFALSE; }

  return kTRUE;
}
//-----------------------------------------------------------------------------
void AliD0toKpi::SetPIDresponse(Double_t resp0[5],Double_t resp1[5]) {
  // Set combined PID detector response probabilities

  fPIDrespEl[0] = resp0[0];
  fPIDrespEl[1] = resp1[0];
  fPIDrespMu[0] = resp0[1];
  fPIDrespMu[1] = resp1[1];
  fPIDrespPi[0] = resp0[2];
  fPIDrespPi[1] = resp1[2];
  fPIDrespKa[0] = resp0[3];
  fPIDrespKa[1] = resp1[3];
  fPIDrespPr[0] = resp0[4];
  fPIDrespPr[1] = resp1[4];

  return;
} 
//----------------------------------------------------------------------------
Double_t AliD0toKpi::LinearInterpolation(Double_t p,Int_t nBins,Double_t Bin,
					 const Double_t *values) const {
  // a linear interpolation method

  Double_t value=0; 
  Double_t slope;

  if(p<0.5*Bin) {
    value = values[0];
  } else if(p>=(nBins-0.5)*Bin) {
    slope = (2*values[nBins-1]-values[nBins-2]-values[nBins-3])/Bin/2;
    value = values[nBins-2]+slope*(p-Bin*(nBins-1.5));
  } else {
    for(Int_t i=0; i<nBins; i++) {
      if(p<(i+0.5)*Bin) {
	slope = (values[i]-values[i-1])/Bin;
	value = values[i-1]+slope*(p-Bin*(i-0.5));
	break;
      }
    }
  }

  if(value<0.) value=0.;
  if(value>1.) value=1.;

  return value;
}
//----------------------------------------------------------------------------
Commit	Line	Data
9af36e3a	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	// for pp and PbPb interactions
	19	// Note: the two decay tracks are labelled: 0 (positive track)
	20	// 1 (negative track)
	21	// Origin: A. Dainese andrea.dainese@lnl.infn.it
	22	//----------------------------------------------------------------------------
	23
	24	#include <TH1.h>
	25	#include <TH2.h>
	26	#include <TCanvas.h>
	27	#include <TPaveLabel.h>
	28	#include <TVector3.h>
	29
	30	#include "AliD0toKpi.h"
	31
	32	ClassImp(AliD0toKpi)
	33
	34	//----------------------------------------------------------------------------
b30302f6	35	AliD0toKpi::AliD0toKpi():
	36	fSignal(kFALSE),
	37	fEvent(0),
	38	fV1x(0.),
	39	fV1y(0.),
	40	fV1z(0.),
	41	fV2x(0.),
	42	fV2y(0.),
	43	fV2z(0.),
	44	fDCA(0.),
	45	fWgtAD0(0.),
	46	fWgtAD0bar(0.),
	47	fWgtBD0(0.),
	48	fWgtBD0bar(0.),
	49	fWgtCD0(0.),
	50	fWgtCD0bar(0.),
	51	fWgtDD0(0.),
	52	fWgtDD0bar(0.)
	53	{
9af36e3a	54	// Default constructor
9af36e3a	55
	56	fTrkNum[0] = 0;
	57	fTrkNum[1] = 0;
	58
9af36e3a	59	fPx[0] = 0.;
	60	fPy[0] = 0.;
	61	fPz[0] = 0.;
	62	fPx[1] = 0.;
	63	fPy[1] = 0.;
	64	fPz[1] = 0.;
	65
	66	fd0[0] = 0.;
	67	fd0[1] = 0.;
	68
	69	fPdg[0] = 0;
	70	fPdg[1] = 0;
	71	fMum[0] = 0;
	72	fMum[1] = 0;
	73
	74	fTagPi[0] = 0.;
	75	fTagPi[1] = 0.;
	76	fTagKa[0] = 0.;
	77	fTagKa[1] = 0.;
	78	fTagNid[0] = 0.;
	79	fTagNid[1] = 0.;
	80
9af36e3a	81	}
	82	//----------------------------------------------------------------------------
	83	AliD0toKpi::AliD0toKpi(Int_t ev,Int_t trkNum[2],
	84	Double_t v1[3],Double_t v2[3],
	85	Double_t dca,
b30302f6	86	Double_t mom[6],Double_t d0[2]):
	87	fSignal(kFALSE),
	88	fEvent(ev),
	89	fV1x(v1[0]),
	90	fV1y(v1[1]),
	91	fV1z(v1[2]),
	92	fV2x(v2[0]),
	93	fV2y(v2[1]),
	94	fV2z(v2[2]),
	95	fDCA(dca),
	96	fWgtAD0(0.),
	97	fWgtAD0bar(0.),
	98	fWgtBD0(0.),
	99	fWgtBD0bar(0.),
	100	fWgtCD0(0.),
	101	fWgtCD0bar(0.),
	102	fWgtDD0(0.),
	103	fWgtDD0bar(0.)
	104	{
9af36e3a	105	// Constructor
9af36e3a	106
9af36e3a	107	fTrkNum[0] = trkNum[0];
	108	fTrkNum[1] = trkNum[1];
	109
9af36e3a	110	fPx[0] = mom[0];
	111	fPy[0] = mom[1];
	112	fPz[0] = mom[2];
	113	fPx[1] = mom[3];
	114	fPy[1] = mom[4];
	115	fPz[1] = mom[5];
	116
	117	fd0[0] = d0[0];
	118	fd0[1] = d0[1];
	119
	120	fPdg[0] = 0;
	121	fPdg[1] = 0;
	122	fMum[0] = 0;
	123	fMum[1] = 0;
	124
	125	fTagPi[0] = 0.;
	126	fTagPi[1] = 0.;
	127	fTagKa[0] = 0.;
	128	fTagKa[1] = 0.;
	129	fTagNid[0] = 0.;
	130	fTagNid[1] = 0.;
9af36e3a	131	}
	132	//----------------------------------------------------------------------------
	133	AliD0toKpi::~AliD0toKpi() {}
9af36e3a	134	//----------------------------------------------------------------------------
	135	void AliD0toKpi::ApplyPID(TString pidScheme) {
	136	// Applies particle identification
	137
	138	if((!pidScheme.CompareTo("TOFparamPbPb") \|\| !pidScheme.CompareTo("TOFparamPP")) && fPdg[0]==0) {
	139	printf("AliD0toKpi::ApplyPID :\n Warning: TOF parameterized PID can be used only for simulation!\n");
	140	return;
	141	}
	142
	143	if(!pidScheme.CompareTo("TOFparamPbPb")) {
	144	// tagging of the positive track
	145	if(TMath::Abs(fPdg[0])==211 \|\| TMath::Abs(fPdg[0])==13
	146	\|\| TMath::Abs(fPdg[0])==11) { // pion,muon,electron
	147	fTagPi[0] = LinearInterpolation(PChild(0),kPiBinsPbPb,kPiBinWidthPbPb,kPiTagPiPbPb);
	148	fTagNid[0] = 1.-fTagPi[0];
	149	fTagKa[0] = 0.;
	150	fTagPr[0] = 0.;
	151	}
	152	if(TMath::Abs(fPdg[0])==321) { // kaon
	153	fTagKa[0] = LinearInterpolation(PChild(0),kKBinsPbPb,kKBinWidthPbPb,kKTagKPbPb);
	154	fTagNid[0] = LinearInterpolation(PChild(0),kKBinsPbPb,kKBinWidthPbPb,kKTagNidPbPb);
	155	if((fTagNid[0]+fTagKa[0])>1.) fTagNid[0] = 1.-fTagKa[0];
	156	fTagPi[0] = 1.-fTagNid[0]-fTagKa[0];
	157	fTagPr[0] = 0.;
	158	}
	159	if(TMath::Abs(fPdg[0])==2212) { // proton
	160	fTagPr[0] = LinearInterpolation(PChild(0),kPBinsPbPb,kPBinWidthPbPb,kPTagPPbPb);
	161	fTagNid[0] = LinearInterpolation(PChild(0),kPBinsPbPb,kPBinWidthPbPb,kPTagNidPbPb);
	162	if((fTagNid[0]+fTagPr[0])>1.) fTagNid[0] = 1.-fTagPr[0];
	163	fTagPi[0] = 1.-fTagNid[0]-fTagPr[0];
	164	fTagKa[0] = 0.;
	165	}
	166	// tagging of the negative track
	167	if(TMath::Abs(fPdg[1])==211 \|\| TMath::Abs(fPdg[1])==13
	168	\|\| TMath::Abs(fPdg[1])==11) { // pion,muon,electron
	169	fTagPi[1] = LinearInterpolation(PChild(1),kPiBinsPbPb,kPiBinWidthPbPb,kPiTagPiPbPb);
	170	fTagNid[1] = 1.-fTagPi[1];
	171	fTagKa[1] = 0.;
	172	fTagPr[1] = 0.;
	173	}
	174	if(TMath::Abs(fPdg[1])==321) { // kaon
	175	fTagKa[1] = LinearInterpolation(PChild(1),kKBinsPbPb,kKBinWidthPbPb,kKTagKPbPb);
	176	fTagNid[1] = LinearInterpolation(PChild(1),kKBinsPbPb,kKBinWidthPbPb,kKTagNidPbPb);
	177	if((fTagNid[1]+fTagKa[1])>1.) fTagNid[1] = 1.-fTagKa[1];
	178	fTagPi[1] = 1.-fTagNid[1]-fTagKa[1];
	179	fTagPr[1] = 0.;
	180	}
	181	if(TMath::Abs(fPdg[1])==2212) { // proton
	182	fTagPr[1] = LinearInterpolation(PChild(1),kPBinsPbPb,kPBinWidthPbPb,kPTagPPbPb);
	183	fTagNid[1] = LinearInterpolation(PChild(1),kPBinsPbPb,kPBinWidthPbPb,kPTagNidPbPb);
	184	if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1];
	185	fTagPi[1] = 1.-fTagNid[1]-fTagPr[1];
	186	fTagKa[1] = 0.;
	187	}
	188	}
	189
	190
	191	if(!pidScheme.CompareTo("TOFparamPP")) {
	192	// tagging of the positive track
	193	if(TMath::Abs(fPdg[0])==211 \|\| TMath::Abs(fPdg[0])==13
	194	\|\| TMath::Abs(fPdg[0])==11) { // pion,muon,electron
	195	fTagPi[0] = LinearInterpolation(PChild(0),kPiBinsPP,kPiBinWidthPP,kPiTagPiPP);
	196	fTagNid[0] = 1.-fTagPi[0];
	197	fTagKa[0] = 0.;
198	fTagPr[0] = 0.;
199	}
200	if(TMath::Abs(fPdg[0])==321) { // kaon
201	fTagKa[0] = LinearInterpolation(PChild(0),kKBinsPP,kKBinWidthPP,kKTagKPP);
202	fTagNid[0] = LinearInterpolation(PChild(0),kKBinsPP,kKBinWidthPP,kKTagNidPP);
203	if((fTagNid[0]+fTagKa[0])>1.) fTagNid[0] = 1.-fTagKa[0];
204	fTagPi[0] = 1.-fTagNid[0]-fTagKa[0];
205	fTagPr[0] = 0.;
206	}
207	if(TMath::Abs(fPdg[0])==2212) { // proton
208	fTagPr[0] = LinearInterpolation(PChild(0),kPBinsPP,kPBinWidthPP,kPTagPPP);
209	fTagNid[0] = LinearInterpolation(PChild(0),kPBinsPP,kPBinWidthPP,kPTagNidPP);
210	if((fTagNid[0]+fTagPr[0])>1.) fTagNid[0] = 1.-fTagPr[0];
211	fTagPi[0] = 1.-fTagNid[0]-fTagPr[0];
212	fTagKa[0] = 0.;
213	}
214	// tagging of the negative track
215	if(TMath::Abs(fPdg[1])==211 \|\| TMath::Abs(fPdg[1])==13
216	\|\| TMath::Abs(fPdg[1])==11) { // pion,muon,electron
217	fTagPi[1] = LinearInterpolation(PChild(1),kPiBinsPP,kPiBinWidthPP,kPiTagPiPP);
218	fTagNid[1] = 1.-fTagPi[1];
219	fTagKa[1] = 0.;
220	fTagPr[1] = 0.;
221	}
222	if(TMath::Abs(fPdg[1])==321) { // kaon
223	fTagKa[1] = LinearInterpolation(PChild(1),kKBinsPP,kKBinWidthPP,kKTagKPP);
224	fTagNid[1] = LinearInterpolation(PChild(1),kKBinsPP,kKBinWidthPP,kKTagNidPP);
225	if((fTagNid[1]+fTagKa[1])>1.) fTagNid[1] = 1.-fTagKa[1];
226	fTagPi[1] = 1.-fTagNid[1]-fTagKa[1];
227	fTagPr[1] = 0.;
228	}
229	if(TMath::Abs(fPdg[1])==2212) { // proton
230	fTagPr[1] = LinearInterpolation(PChild(1),kPBinsPP,kPBinWidthPP,kPTagPPP);
231	fTagNid[1] = LinearInterpolation(PChild(1),kPBinsPP,kPBinWidthPP,kPTagNidPP);
232	if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1];
233	fTagPi[1] = 1.-fTagNid[1]-fTagPr[1];
234	fTagKa[1] = 0.;
235	}
236	}
237
238	return;
239	}
240	//----------------------------------------------------------------------------
241	Double_t AliD0toKpi::ChildrenRelAngle() const {
242	// relative angle between K and pi
243
244	TVector3 mom0(fPx[0],fPy[0],fPz[0]);
245	TVector3 mom1(fPx[1],fPy[1],fPz[1]);
246
247	Double_t angle = mom0.Angle(mom1);
248
249	return angle;
250	}
251	//----------------------------------------------------------------------------
252	void AliD0toKpi::ComputeWgts() {
253	// calculate the weights for PID
254
255
256	// assignement of the weights from PID
257	fWgtAD0 = fTagKa[1]*(fTagPi[0]+fTagNid[0]);
258	fWgtAD0bar = fTagKa[0]*(fTagPi[1]+fTagNid[1]);
259	fWgtBD0 = fTagPi[0]*fTagNid[1];
260	fWgtBD0bar = fTagPi[1]*fTagNid[0];
261	fWgtCD0 = fTagNid[0]*fTagNid[1];
262	fWgtCD0bar = fTagNid[0]*fTagNid[1];
263	fWgtDD0 = 1.-fWgtAD0-fWgtBD0-fWgtCD0;
264	fWgtDD0bar = 1.-fWgtAD0bar-fWgtBD0bar-fWgtCD0bar;
265
266	/*
267	cerr<<fWgtAD0<<" "<<fWgtAD0bar<<endl;
268	cerr<<fWgtBD0<<" "<<fWgtBD0bar<<endl;
269	cerr<<fWgtCD0<<" "<<fWgtCD0bar<<endl;
270	cerr<<fWgtDD0<<" "<<fWgtDD0bar<<endl;
271	*/
272	/*
273	if(fWgtAD0<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
274	if(fWgtAD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
275	if(fWgtBD0<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
276	if(fWgtBD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
277	if(fWgtCD0<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
278	if(fWgtCD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
279	*/
280
281	return;
282	}
283	//----------------------------------------------------------------------------
284	void AliD0toKpi::CorrectWgt4BR(Double_t factor) {
285	// correct weights of background from charm
286
287	fWgtAD0 *= factor;
288	fWgtAD0bar *= factor;
289	fWgtBD0 *= factor;
290	fWgtBD0bar *= factor;
291	fWgtCD0 *= factor;
292	fWgtCD0bar *= factor;
293	fWgtDD0 *= factor;
294	fWgtDD0bar *= factor;
295
296	return;
297	}
298	//----------------------------------------------------------------------------
299	Double_t AliD0toKpi::CosPointing() const {
300	// cosine of pointing angle in space
301
302	TVector3 mom(Px(),Py(),Pz());
303	TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
304
305	Double_t pta = mom.Angle(flight);
306
307	return TMath::Cos(pta);
308	}
309	//----------------------------------------------------------------------------
310	Double_t AliD0toKpi::CosPointingXY() const {
311	// cosine of pointing angle in transverse plane
312
313	TVector3 momXY(Px(),Py(),0.);
314	TVector3 flightXY(fV2x-fV1x,fV2y-fV1y,0.);
315
316	Double_t ptaXY = momXY.Angle(flightXY);
317
318	return TMath::Cos(ptaXY);
319	}
320	//----------------------------------------------------------------------------
321	void AliD0toKpi::CosThetaStar(Double_t &ctsD0,Double_t &ctsD0bar) const {
322	// cosine of decay angle in the D0 rest frame
323
324	Double_t pStar = TMath::Sqrt(TMath::Power(kMD0kMD0-kMKkMK-kMPikMPi,2.)-4.kMKkMKkMPikMPi)/(2.kMD0);
325
326	Double_t beta = P()/Energy();
327	Double_t gamma = Energy()/kMD0;
328
329	ctsD0 = (Ql(1)/gamma-betaTMath::Sqrt(pStarpStar+kMK*kMK))/pStar;
330	// if(ctsD0 > 1.) { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0<<"!\n"; }
331	// if(ctsD0 < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0<<"!\n"; }
332
333	ctsD0bar = (Ql(0)/gamma-betaTMath::Sqrt(pStarpStar+kMK*kMK))/pStar;
334	// if(ctsD0bar > 1.) { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0bar<<"!\n"; }
335	// if(ctsD0bar < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0bar<<"!\n";}
336
337	return;
338	}
339	//----------------------------------------------------------------------------
340	Double_t AliD0toKpi::Eta() const {
341	// pseudorapidity of the D0
342
343	Double_t theta = TMath::Pi()/2.-TMath::ATan2(Pz(),Pt());
344	Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
345	return eta;
346	}
347	//----------------------------------------------------------------------------
348	Double_t AliD0toKpi::EtaChild(Int_t child) const {
349	// pseudorapidity of the decay tracks
350
351	Double_t theta = TMath::Pi()/2.-TMath::ATan2(fPz[child],PtChild(child));
352	Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
353	return eta;
354	}
355	//----------------------------------------------------------------------------
356	void AliD0toKpi::GetWgts(Double_t &WgtD0,Double_t &WgtD0bar,TString sample)
357	const {
358	// returns the weights for pid
359
360	if(!sample.CompareTo("A")) { WgtD0 = fWgtAD0; WgtD0bar = fWgtAD0bar; }
361	if(!sample.CompareTo("B")) { WgtD0 = fWgtBD0; WgtD0bar = fWgtBD0bar; }
362	if(!sample.CompareTo("C")) { WgtD0 = fWgtCD0; WgtD0bar = fWgtCD0bar; }
363	if(!sample.CompareTo("D")) { WgtD0 = fWgtDD0; WgtD0bar = fWgtDD0bar; }
364	if(!sample.CompareTo("ABCD")) {
365	WgtD0 = fWgtAD0+fWgtBD0+fWgtCD0+fWgtDD0;
366	WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar+fWgtDD0bar;
367	}
368	if(!sample.CompareTo("ABC")) {
369	WgtD0 = fWgtAD0+fWgtBD0+fWgtCD0;
370	WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar;
371	}
372	if(!sample.CompareTo("BC")) {
373	WgtD0 = fWgtBD0+fWgtCD0;
374	WgtD0bar = fWgtBD0bar+fWgtCD0bar;
375	}
376
377	return;
378	}
379	//----------------------------------------------------------------------------
380	Double_t AliD0toKpi::ImpPar() const {
381	// D0 impact parameter in the bending plane
382
383	Double_t k = -(fV2x-fV1x)Px()-(fV2y-fV1y)Py();
384	k /= Pt()*Pt();
385	Double_t dx = fV2x-fV1x+k*Px();
386	Double_t dy = fV2y-fV1y+k*Py();
387	Double_t absDD = TMath::Sqrt(dxdx+dydy);
388	TVector3 mom(Px(),Py(),Pz());
389	TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
390	TVector3 cross = mom.Cross(flight);
391	return (cross.Z()>0. ? absDD : -absDD);
392	}
393	//----------------------------------------------------------------------------
394	void AliD0toKpi::InvMass(Double_t &mD0,Double_t &mD0bar) const {
395	// invariant mass as D0 and as D0bar
396
397	Double_t energy[2];
398
399	// D0 -> K- Pi+
400	energy[1] = TMath::Sqrt(kMKkMK+PChild(1)PChild(1));
401	energy[0] = TMath::Sqrt(kMPikMPi+PChild(0)PChild(0));
402
403	mD0 = TMath::Sqrt((energy[0]+energy[1])(energy[0]+energy[1])-P()P());
404
405
406	// D0bar -> K+ Pi-
407	energy[0] = TMath::Sqrt(kMKkMK+PChild(0)PChild(0));
408	energy[1] = TMath::Sqrt(kMPikMPi+PChild(1)PChild(1));
409
410	mD0bar = TMath::Sqrt((energy[0]+energy[1])(energy[0]+energy[1])-P()P());
411
412	return;
413
414	}
415	//----------------------------------------------------------------------------
416	Double_t AliD0toKpi::Ql(Int_t child) const {
417	// longitudinal momentum of decay tracks w.r.t. to D0 momentum
418
419	Double_t qL;
420	TVector3 mom(fPx[child],fPy[child],fPz[child]);
421	TVector3 momD(Px(),Py(),Pz());
422
423	qL = mom.Dot(momD)/momD.Mag();
424
425	return qL ;
426	}
427	//----------------------------------------------------------------------------
428	Double_t AliD0toKpi::Qt() const {
429	// transverse momentum of decay tracks w.r.t. to D0 momentum
430
431	TVector3 mom0(fPx[0],fPy[0],fPz[0]);
432	TVector3 momD(Px(),Py(),Pz());
433
434	return mom0.Perp(momD);
435	}
436	//----------------------------------------------------------------------------
437	Bool_t AliD0toKpi::Select(const Double_t* cuts,Int_t& okD0,Int_t& okD0bar)
438	const {
439	//
440	// This function compares the D0 with a set of cuts:
441	//
442	// cuts[0] = inv. mass half width [GeV]
443	// cuts[1] = dca [micron]
444	// cuts[2] = cosThetaStar
445	// cuts[3] = pTK [GeV/c]
446	// cuts[4] = pTPi [GeV/c]
447	// cuts[5] = d0K [micron] upper limit!
448	// cuts[6] = d0Pi [micron] upper limit!
449	// cuts[7] = d0d0 [micron^2]
450	// cuts[8] = cosThetaPoint
451	//
452	// If the D0/D0bar doesn't pass the cuts it sets the weights to 0
453	// If neither D0 nor D0bar pass the cuts return kFALSE
454	//
455	Double_t mD0,mD0bar,ctsD0,ctsD0bar;
456	okD0=1; okD0bar=1;
457
458	if(PtChild(1) < cuts[3] \|\| PtChild(0) < cuts[4]) okD0 = 0;
459	if(PtChild(0) < cuts[3] \|\| PtChild(1) < cuts[4]) okD0bar = 0;
460	if(!okD0 && !okD0bar) return kFALSE;
461
462	if(TMath::Abs(Getd0Child(1)) > cuts[5] \|\|
463	TMath::Abs(Getd0Child(0)) > cuts[6]) okD0 = 0;
464	if(TMath::Abs(Getd0Child(0)) > cuts[6] \|\|
465	TMath::Abs(Getd0Child(1)) > cuts[5]) okD0bar = 0;
466	if(!okD0 && !okD0bar) return kFALSE;
467
468	if(GetDCA() > cuts[1]) { okD0 = okD0bar = 0; return kFALSE; }
469
470	InvMass(mD0,mD0bar);
471	if(TMath::Abs(mD0-kMD0) > cuts[0]) okD0 = 0;
472	if(TMath::Abs(mD0bar-kMD0) > cuts[0]) okD0bar = 0;
473	if(!okD0 && !okD0bar) return kFALSE;
474
475	CosThetaStar(ctsD0,ctsD0bar);
476	if(TMath::Abs(ctsD0) > cuts[2]) okD0 = 0;
477	if(TMath::Abs(ctsD0bar) > cuts[2]) okD0bar = 0;
478	if(!okD0 && !okD0bar) return kFALSE;
479
480	if(ProdImpParams() > cuts[7]) { okD0 = okD0bar = 0; return kFALSE; }
481
482	if(CosPointing() < cuts[8]) { okD0 = okD0bar = 0; return kFALSE; }
483
484	return kTRUE;
485	}
486	//-----------------------------------------------------------------------------
487	void AliD0toKpi::SetPIDresponse(Double_t resp0[5],Double_t resp1[5]) {
488	// Set combined PID detector response probabilities
489
490	fPIDrespEl[0] = resp0[0];
491	fPIDrespEl[1] = resp1[0];
492	fPIDrespMu[0] = resp0[1];
493	fPIDrespMu[1] = resp1[1];
494	fPIDrespPi[0] = resp0[2];
495	fPIDrespPi[1] = resp1[2];
496	fPIDrespKa[0] = resp0[3];
497	fPIDrespKa[1] = resp1[3];
498	fPIDrespPr[0] = resp0[4];
499	fPIDrespPr[1] = resp1[4];
500
501	return;
502	}
503	//----------------------------------------------------------------------------
504	Double_t AliD0toKpi::LinearInterpolation(Double_t p,Int_t nBins,Double_t Bin,
505	const Double_t *values) const {
506	// a linear interpolation method
507
508	Double_t value=0;
509	Double_t slope;
510
511	if(p<0.5*Bin) {
512	value = values[0];
513	} else if(p>=(nBins-0.5)*Bin) {
514	slope = (2*values[nBins-1]-values[nBins-2]-values[nBins-3])/Bin/2;
515	value = values[nBins-2]+slope(p-Bin(nBins-1.5));
516	} else {
517	for(Int_t i=0; i<nBins; i++) {
518	if(p<(i+0.5)*Bin) {
519	slope = (values[i]-values[i-1])/Bin;
520	value = values[i-1]+slope(p-Bin(i-0.5));
521	break;
522	}
523	}
524	}
525
526	if(value<0.) value=0.;
527	if(value>1.) value=1.;
528
529	return value;
530	}
531	//----------------------------------------------------------------------------
532
533
534
535
536
537
538
539
540