From 9af36e3ac08b68035ad812be55346b5fc5da9566 Mon Sep 17 00:00:00 2001 From: dainese Date: Fri, 20 Apr 2007 07:46:41 +0000 Subject: [PATCH] move D0toKpi to PWG3 (Andrea) --- PWG3/AliD0toKpi.cxx | 535 +++++++++++++++++++++++++++++ PWG3/AliD0toKpi.h | 195 +++++++++++ PWG3/AliD0toKpiAnalysis.cxx | 657 ++++++++++++++++++++++++++++++++++++ PWG3/AliD0toKpiAnalysis.h | 90 +++++ PWG3/AliD0toKpiPlots.C | 235 +++++++++++++ PWG3/AliD0toKpiTest.C | 6 +- 6 files changed, 1716 insertions(+), 2 deletions(-) create mode 100644 PWG3/AliD0toKpi.cxx create mode 100644 PWG3/AliD0toKpi.h create mode 100644 PWG3/AliD0toKpiAnalysis.cxx create mode 100644 PWG3/AliD0toKpiAnalysis.h create mode 100644 PWG3/AliD0toKpiPlots.C diff --git a/PWG3/AliD0toKpi.cxx b/PWG3/AliD0toKpi.cxx new file mode 100644 index 00000000000..31a4c94d156 --- /dev/null +++ b/PWG3/AliD0toKpi.cxx @@ -0,0 +1,535 @@ +/************************************************************************** + * 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 +#include +#include +#include +#include + +#include "AliD0toKpi.h" + +ClassImp(AliD0toKpi) + +//---------------------------------------------------------------------------- +AliD0toKpi::AliD0toKpi() { + // Default constructor + + fSignal = kFALSE; + + fEvent = 0; + + fTrkNum[0] = 0; + fTrkNum[1] = 0; + + fV1x = 0.; + fV1y = 0.; + fV1z = 0.; + fV2x = 0.; + fV2y = 0.; + fV2z = 0.; + fDCA = 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.; + + fWgtAD0=fWgtAD0bar=fWgtBD0=fWgtBD0bar=fWgtCD0=fWgtCD0bar=fWgtDD0=fWgtDD0bar=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]) { + // Constructor + + fSignal = kFALSE; + + fEvent = ev; + fTrkNum[0] = trkNum[0]; + fTrkNum[1] = trkNum[1]; + + fV1x = v1[0]; + fV1y = v1[1]; + fV1z = v1[2]; + fV2x = v2[0]; + fV2y = v2[1]; + fV2z = v2[2]; + fDCA = dca; + + 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.; + + fWgtAD0=fWgtAD0bar=fWgtBD0=fWgtBD0bar=fWgtCD0=fWgtCD0bar=fWgtDD0=fWgtDD0bar=0; +} +//---------------------------------------------------------------------------- +AliD0toKpi::~AliD0toKpi() {} +//____________________________________________________________________________ +AliD0toKpi::AliD0toKpi( const AliD0toKpi& d0toKpi):TObject(d0toKpi) { + // dummy copy constructor +} +//---------------------------------------------------------------------------- +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< 1.) { cerr<<"AliD0toKpi::CosThetaStar: > 1 "< 1.) { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<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; i1.) value=1.; + + return value; +} +//---------------------------------------------------------------------------- + + + + + + + + + diff --git a/PWG3/AliD0toKpi.h b/PWG3/AliD0toKpi.h new file mode 100644 index 00000000000..716dbb9d4dd --- /dev/null +++ b/PWG3/AliD0toKpi.h @@ -0,0 +1,195 @@ +#ifndef AliD0toKpi_H +#define AliD0toKpi_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//------------------------------------------------------------------------- +// Class AliD0toKpi +// Reconstructed D0 -> K^- pi^+ class +// +// Note: the two decay tracks are labelled: 0 (positive track) +// 1 (negative track) +// +// Origin: A. Dainese andrea.dainese@lnl.infn.it +//------------------------------------------------------------------------- + +#include +#include + +//---------------------------------------------------------------------------- +// Some constants (masses + parameterized TOF PID) +// +// particle masses +const Double_t kMD0 = 1.8645; // D0 mass +const Double_t kMK = 0.49368; // K+ mass +const Double_t kMPi = 0.13957; // pi+ mass + +// --- TOF tagging probabilities --- +// central HIJING +// B = 0.4 T +// tracking errors in TPC included +// With TRD +// +// *** Pb-Pb dNch/dy=6000 *** +// +// PIONS +const Int_t kPiBinsPbPb = 10; +const Double_t kPiBinWidthPbPb = 0.250; +const Double_t kPiTagPiPbPb[kPiBinsPbPb] = {0.211421,0.652184,0.624421,0.614727,0.610777,0.628015,0.631520,0.630324,0.637551,0.575235}; +const Double_t kPiTagNidPbPb[kPiBinsPbPb] = {0.788579,0.347816,0.375579,0.385273,0.389223,0.371985,0.368480,0.369676,0.362449,0.424765}; +// KAONS +const Int_t kKBinsPbPb = 10; +const Double_t kKBinWidthPbPb = 0.250; +const Double_t kKTagKPbPb[kKBinsPbPb] = {0.000000,0.101255,0.397662,0.467586,0.517008,0.555023,0.584185,0.519029,0.464117,0.247308}; +const Double_t kKTagPiPbPb[kKBinsPbPb] = {0.102049,0.289930,0.101930,0.057771,0.040286,0.028567,0.053108,0.094369,0.066302,0.247308}; +const Double_t kKTagNidPbPb[kKBinsPbPb] = {0.897951,0.608815,0.500408,0.474643,0.442705,0.416410,0.362707,0.386603,0.469580,0.505383}; +// PROTONS +const Int_t kPBinsPbPb = 9; +const Double_t kPBinWidthPbPb = 0.500; +const Double_t kPTagPPbPb[kPBinsPbPb] = {0.017940,0.350681,0.535286,0.583264,0.562935,0.560524,0.545992,0.598060,0.351245}; +const Double_t kPTagPiPbPb[kPBinsPbPb] = {0.195955,0.094949,0.039962,0.026039,0.007556,0.016986,0.030333,0.000000,0.000000}; +const Double_t kPTagNidPbPb[kPBinsPbPb] = {0.786105,0.554370,0.424751,0.390697,0.429508,0.422491,0.423675,0.401940,0.648755}; +// +// pp PYTHIA +// +// *** cuts for pp *** +// +// PIONS +const Int_t kPiBinsPP = 10; +const Double_t kPiBinWidthPP = 0.250; +const Double_t kPiTagPiPP[kPiBinsPP] = {0.194528,0.447097,0.603364,0.646413,0.647125,0.669157,0.688139,0.682564,0.689910,0.665710}; +const Double_t kPiTagNidPP[kPiBinsPP] = {0.805472,0.552903,0.396636,0.353587,0.352875,0.330843,0.311861,0.317436,0.310090,0.334290}; +// KAONS +const Int_t kKBinsPP = 10; +const Double_t kKBinWidthPP = 0.250; +const Double_t kKTagKPP[kKBinsPP] = {0.000000,0.173393,0.439690,0.519423,0.587025,0.605372,0.586021,0.650139,0.444444,0.299363}; +const Double_t kKTagPiPP[kKBinsPP] = {0.000000,0.001495,0.000000,-0.000000,-0.000000,0.000000,0.032258,0.060572,0.101449,0.242038}; +const Double_t kKTagNidPP[kKBinsPP] = {1.000000,0.825112,0.560310,0.480577,0.412975,0.394628,0.381720,0.289289,0.454106,0.458599}; +// PROTONS +const Int_t kPBinsPP = 9; +const Double_t kPBinWidthPP = 0.500; +const Double_t kPTagPPP[kPBinsPP] = {0.029404,0.438640,0.613710,0.665152,0.634961,0.657711,0.703704,0.685714,0.235294}; +const Double_t kPTagPiPP[kPBinsPP] = {0.000000,0.000000,0.000000,-0.000000,0.000000,0.000000,-0.000000,0.014286,-0.000000}; +const Double_t kPTagNidPP[kPBinsPP] = {0.970596,0.561360,0.386290,0.334848,0.365039,0.342289,0.296296,0.300000,0.764706}; + + + + +//----------------------------------------------------------------------------- +class AliD0toKpi : public TObject { + public: + // + 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]); + virtual ~AliD0toKpi(); + AliD0toKpi(const AliD0toKpi& d0toKpi); + + Double_t Alpha() const { return (Ql(0)-Ql(1))/(Ql(0)+Ql(1)); } + void ApplyPID(TString pidScheme="TOFparamPbPb"); + Double_t ChildrenRelAngle() const; + void ComputeWgts(); + void CorrectWgt4BR(Double_t factor); + Double_t CosPointing() const; + Double_t CosPointingXY() const; + void CosThetaStar(Double_t &ctsD0,Double_t &ctsD0bar) const; + Double_t Ct() const {return Length()*kMD0/P();} + Double_t Energy() const { return TMath::Sqrt(P()*P()+kMD0*kMD0); } + Double_t Eta() const; + Double_t EtaChild(Int_t child) const; + Int_t EventNo() const {return TMath::Abs(fEvent);} + Double_t GetDCA() const { return 10000.*fDCA; } + Int_t GetTrkNum(Int_t child) const { return fTrkNum[child]; } + Double_t Getd0Child(Int_t child) const { return fd0[child]; } + Int_t GetPdgChild(Int_t child) const { return fPdg[child]; } + Int_t GetPdgMum(Int_t child) const {return fMum[child]; } + void GetWgts(Double_t &WgtD0,Double_t &WgtD0bar,TString sample) const; + void GetPrimaryVtx(Double_t vtx[3]) const + { vtx[0]=fV1x; vtx[1]=fV1y; vtx[2]=fV1z; } + void GetSecondaryVtx(Double_t vtx[3]) const + { vtx[0]=fV2x; vtx[1]=fV2y; vtx[2]=fV2z; } + + Double_t ImpPar() const; + void InvMass(Double_t &mD0,Double_t &mD0bar) const; + Bool_t IsSignal() const { if(fSignal) return kTRUE; return kFALSE; } + Double_t Length() const + { return TMath::Sqrt((fV1x-fV2x)*(fV1x-fV2x) + +(fV1y-fV2y)*(fV1y-fV2y)+(fV1z-fV2z)*(fV1z-fV2z)); } + Double_t P() const { return TMath::Sqrt(Pt()*Pt()+Pz()*Pz()); } + Double_t PChild(Int_t child) const { return TMath::Sqrt(fPx[child]*fPx[child]+fPy[child]*fPy[child]+fPz[child]*fPz[child]); } + Double_t ProdImpParams() const { return fd0[0]*fd0[1]; } + Double_t Pt() const { return TMath::Sqrt(Px()*Px()+Py()*Py()); } + Double_t PtChild(Int_t child) const { return TMath::Sqrt(fPx[child]*fPx[child]+fPy[child]*fPy[child]); } + Double_t Px() const { return (fPx[0]+fPx[1]); } + Double_t Py() const { return (fPy[0]+fPy[1]); } + Double_t Pz() const { return (fPz[0]+fPz[1]); } + Double_t Ql(Int_t child) const; + Double_t Qt() const; + Double_t Rapidity() const { return 0.5*TMath::Log((Energy()+Pz())/(Energy()-Pz()+1.e-13)); } + Bool_t Select(const Double_t* cuts,Int_t& okD0,Int_t& okD0bar) const; + void SetPrimaryVtx(Double_t vtx[3]) + { fV1x=vtx[0]; fV1y=vtx[1]; fV1z=vtx[2]; } + void SetSignal() { fSignal = kTRUE; } + void SetTOFmasses(Double_t mass[2]) + { fTOFmass[0]=mass[0]; fTOFmass[1]=mass[1]; } + void SetPIDresponse(Double_t resp0[5],Double_t resp1[5]); + void SetPdgCodes(Int_t pdg[2]) {fPdg[0]=pdg[0];fPdg[1]=pdg[1]; } + void SetMumPdgCodes(Int_t mum[2]) {fMum[0]=mum[0];fMum[1]=mum[1]; } + Double_t LinearInterpolation(Double_t p,Int_t nBins,Double_t Bin, + const Double_t *values) const; + // + private: + // + Bool_t fSignal; // TRUE if signal, FALSE if background (for simulation) + Int_t fEvent; // number of the event this D0 comes from + // -1 if the D0 comes from ev. mixing + + Int_t fTrkNum[2]; // numbers of the two decay tracks + + Double_t fV1x; // X-position of the primary vertex of the event + Double_t fV1y; // Y-position of the primary vertex of the event + Double_t fV1z; // Z-position of the primary vertex of the event + Double_t fV2x; // X-position of the reconstructed secondary vertex + Double_t fV2y; // Y-position of the reconstructed secondary vertex + Double_t fV2z; // Z-position of the reconstructed secondary vertex + Double_t fDCA; // DCA of the two tracks + + Double_t fPx[2]; // X,Y,Z + Double_t fPy[2]; // momenta of the two tracks + Double_t fPz[2]; // at the reconstructed vertex + + Double_t fd0[2]; // impact parameters in the bending plane + + Int_t fPdg[2]; // PDG codes of the two tracks (for sim.) + Int_t fMum[2]; // PDG codes of the mothers (for sim.) + + Double_t fTagPi[2]; // probability to be tagged as pion + Double_t fTagKa[2]; // probability to be tagged as kaon + Double_t fTagPr[2]; // probability to be tagged as proton + Double_t fTagNid[2]; // probability to be tagged as "non-identified" + + Double_t fPIDrespEl[2]; // det. response to be electron + Double_t fPIDrespMu[2]; // det. response to be muon + Double_t fPIDrespPi[2]; // det. response to be pion + Double_t fPIDrespKa[2]; // det. response to be kaon + Double_t fPIDrespPr[2]; // det. response to be proton + Double_t fTOFmass[2]; // mass estimated by the TOF (-1000. if track not reached TOF) + + Double_t fWgtAD0,fWgtAD0bar; // weights for the 3 samples + Double_t fWgtBD0,fWgtBD0bar; // weights for the 3 samples + Double_t fWgtCD0,fWgtCD0bar; // A: (K,Pi)+(K,?) B: (?,Pi) C: (?,?) + Double_t fWgtDD0,fWgtDD0bar; // D: all other pairs + + ClassDef(AliD0toKpi,1) // Reconstructed D0 candidate class +}; + +#endif + + + + + + + + diff --git a/PWG3/AliD0toKpiAnalysis.cxx b/PWG3/AliD0toKpiAnalysis.cxx new file mode 100644 index 00000000000..e9294d2dd6c --- /dev/null +++ b/PWG3/AliD0toKpiAnalysis.cxx @@ -0,0 +1,657 @@ +/************************************************************************** + * 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 reconstruction and analysis class +// Note: the two decay tracks are labelled: 0 (positive track) +// 1 (negative track) +// An example of usage can be found in the macro AliD0toKpiTest.C +// Origin: A. Dainese andrea.dainese@lnl.infn.it +//---------------------------------------------------------------------------- +#include +#include +#include +#include +#include +#include +#include "AliESD.h" +#include "AliMC.h" +#include "AliRun.h" +#include "AliRunLoader.h" +#include "AliVertexerTracks.h" +#include "AliESDVertex.h" +#include "AliESDv0.h" +#include "AliD0toKpi.h" +#include "AliD0toKpiAnalysis.h" +#include "AliLog.h" + +typedef struct { + Int_t lab; + Int_t pdg; + Int_t mumlab; + Int_t mumpdg; + Int_t mumprongs; + Float_t Vx,Vy,Vz; + Float_t Px,Py,Pz; +} REFTRACK; + +ClassImp(AliD0toKpiAnalysis) + +//---------------------------------------------------------------------------- +AliD0toKpiAnalysis::AliD0toKpiAnalysis() { + // Default constructor + + SetPtCut(); + Setd0Cut(); + SetMassCut(); + SetD0Cuts(); + SetVertex1(); + SetPID(); + fVertexOnTheFly = kFALSE; + fSim = kFALSE; + fOnlySignal = kFALSE; +} +//---------------------------------------------------------------------------- +AliD0toKpiAnalysis::~AliD0toKpiAnalysis() {} +//---------------------------------------------------------------------------- +void AliD0toKpiAnalysis::ApplySelection(const Char_t *inName,const Char_t *outName) const { + // select candidates that pass fD0Cuts and write them to a new file + + TFile *inFile = TFile::Open(inName); + + TTree *treeD0in=(TTree*)inFile->Get("TreeD0"); + AliD0toKpiAnalysis *inAnalysis = (AliD0toKpiAnalysis*)inFile->Get("D0toKpiAnalysis"); + printf("+++\n+++ I N P U T S T A T U S:\n+++\n"); + inAnalysis->PrintStatus(); + + + AliD0toKpi *d = 0; + treeD0in->SetBranchAddress("D0toKpi",&d); + Int_t entries = (Int_t)treeD0in->GetEntries(); + + printf("+++\n+++ Number of D0 in input tree: %d\n+++\n",entries); + + TTree *treeD0out = new TTree("TreeD0","Tree with selected D0 candidates"); + treeD0out->Branch("D0toKpi","AliD0toKpi",&d,200000,0); + + + Int_t okD0=0,okD0bar=0; + Int_t nSel = 0; + + for(Int_t i=0; iGetEvent(i); + + if(fSim && fOnlySignal && !d->IsSignal()) continue; + + // check if candidate passes selection (as D0 or D0bar) + if(d->Select(fD0Cuts,okD0,okD0bar)) { + nSel++; + treeD0out->Fill(); + } + + } + + AliD0toKpiAnalysis *outAnalysis = (AliD0toKpiAnalysis*)inAnalysis->Clone("D0toKpiAnalysis"); + outAnalysis->SetD0Cuts(fD0Cuts); + printf("------------------------------------------\n"); + printf("+++\n+++ O U T P U T S T A T U S:\n+++\n"); + outAnalysis->PrintStatus(); + + printf("+++\n+++ Number of D0 in output tree: %d\n+++\n",nSel); + + TFile* outFile = new TFile(outName,"recreate"); + treeD0out->Write(); + outAnalysis->Write(); + outFile->Close(); + + return; +} +//---------------------------------------------------------------------------- +Double_t AliD0toKpiAnalysis::CalculateTOFmass(Double_t mom,Double_t length, + Double_t time) const { + // calculated the mass from momentum, track length from vertex to TOF + // and time measured by the TOF + if(length==0.) return -1000.; + Double_t a = time*time/length/length; + if(a > 1.) { + a = TMath::Sqrt(a-1.); + } else { + a = -TMath::Sqrt(1.-a); + } + + return mom*a; +} +//---------------------------------------------------------------------------- +void AliD0toKpiAnalysis::FindCandidates(Int_t evFirst,Int_t evLast, + const Char_t *outName) { + // Find D0 candidates and calculate parameters + + + TString esdName="AliESDs.root"; + if(gSystem->AccessPathName(esdName.Data(),kFileExists)) { + printf("AliD0toKpiAnalysis::FindCandidatesESD(): No ESDs file found!\n"); + return; + } + + TString outName1=outName; + TString outName2="nTotEvents.dat"; + + Int_t nTotEv=0,nD0rec=0,nD0rec1ev=0; + Double_t dca; + Double_t v2[3],mom[6],d0[2]; + Int_t iTrkP,iTrkN,trkEntries; + Int_t nTrksP=0,nTrksN=0; + Int_t trkNum[2]; + Double_t tofmass[2]; + Int_t okD0=0,okD0bar=0; + AliESDtrack *postrack = 0; + AliESDtrack *negtrack = 0; + + // create the AliVertexerTracks object + // (it will be used only if fVertexOnTheFly=kTrue) + AliVertexerTracks *vertexer1 = new AliVertexerTracks; + if(fVertexOnTheFly) { + // open the mean vertex + TFile *invtx = new TFile("AliESDVertexMean.root"); + AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean"); + invtx->Close(); + vertexer1->SetVtxStart(initVertex); + delete invtx; + } + Int_t skipped[2]; + Bool_t goodVtx1; + + // create tree for reconstructed D0s + AliD0toKpi *ioD0toKpi=0; + TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates"); + treeD0->Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0); + + // open file with tracks + TFile *esdFile = TFile::Open(esdName.Data()); + AliESD* event = new AliESD; + TTree* tree = (TTree*) esdFile->Get("esdTree"); + if(!tree) { + Error("FindCandidatesESD", "no ESD tree found"); + return; + } + tree->SetBranchAddress("ESD",&event); + + // loop on events in file + for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) { + if(iEvent > evLast) break; + tree->GetEvent(iEvent); + Int_t ev = (Int_t)event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number. + printf("--- Finding D0 -> Kpi in event %d\n",ev); + // count the total number of events + nTotEv++; + + trkEntries = (Int_t)event->GetNumberOfTracks(); + printf(" Number of tracks: %d\n",trkEntries); + if(trkEntries<2) continue; + + // retrieve primary vertex from file + if(!event->GetPrimaryVertex()) { + printf(" No vertex\n"); + continue; + } + event->GetPrimaryVertex()->GetXYZ(fV1); + + // call function which applies sigle-track selection and + // separetes positives and negatives + TObjArray trksP(trkEntries/2); + Int_t *trkEntryP = new Int_t[trkEntries]; + TObjArray trksN(trkEntries/2); + Int_t *trkEntryN = new Int_t[trkEntries]; + TTree *trkTree = new TTree(); + SelectTracks(event,trksP,trkEntryP,nTrksP, + trksN,trkEntryN,nTrksN); + + printf(" pos. tracks: %d neg .tracks: %d\n",nTrksP,nTrksN); + + + nD0rec1ev = 0; + + // LOOP ON POSITIVE TRACKS + for(iTrkP=0; iTrkPGetMagneticField(); + AliESDtrack nt(*negtrack), pt(*postrack); + dca = nt.PropagateToDCA(&pt,b); + + // define the AliESDv0 object + AliESDv0 vertex2(nt,trkNum[0],pt,trkNum[1]); + + // get position of the secondary vertex + vertex2.GetXYZ(v2[0],v2[1],v2[2]); + vertex2.GetPPxPyPz(mom[0],mom[1],mom[2]); + vertex2.GetNPxPyPz(mom[3],mom[4],mom[5]); + // impact parameters of the tracks w.r.t. the primary vertex + d0[0] = 10000.*pt.GetD(fV1[0],fV1[1],b); + d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b); + goodVtx1 = kTRUE; + + // no vertexing if DeltaMass > fMassCut + if(fVertexOnTheFly) { + goodVtx1 = kFALSE; + if(SelectInvMass(mom)) { + // primary vertex from *other* tracks in the event + skipped[0] = trkEntryP[iTrkP]; + skipped[1] = trkEntryN[iTrkN]; + vertexer1->SetSkipTracks(2,skipped); + AliESDVertex *vertex1onfly = + (AliESDVertex*)vertexer1->FindPrimaryVertex(event); + if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE; + vertex1onfly->GetXYZ(fV1); + //vertex1onfly->PrintStatus(); + delete vertex1onfly; + } + } + + + // create the object AliD0toKpi + AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0); + // select D0s + if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) { + // get PID info from ESD + AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]); + Double_t esdpid0[5]; + t0->GetESDpid(esdpid0); + if(t0->GetStatus()&AliESDtrack::kTOFpid) { + tofmass[0] = CalculateTOFmass(t0->GetP(), + t0->GetIntegratedLength(), + t0->GetTOFsignal()); + } else { + tofmass[0] = -1000.; + } + AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]); + Double_t esdpid1[5]; + t1->GetESDpid(esdpid1); + if(t1->GetStatus()&AliESDtrack::kTOFpid) { + tofmass[1] = CalculateTOFmass(t1->GetP(), + t1->GetIntegratedLength(), + t1->GetTOFsignal()); + } else { + tofmass[1] = -1000.; + } + + theD0.SetPIDresponse(esdpid0,esdpid1); + theD0.SetTOFmasses(tofmass); + + // fill the tree + ioD0toKpi=&theD0; + treeD0->Fill(); + + nD0rec++; nD0rec1ev++; + ioD0toKpi=0; + } + + negtrack = 0; + } // loop on negative tracks + postrack = 0; + } // loop on positive tracks + + delete [] trkEntryP; + delete [] trkEntryN; + delete trkTree; + + printf(" Number of D0 candidates: %d\n",nD0rec1ev); + } // loop on events in file + + + printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv); + printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec); + + delete vertexer1; + + esdFile->Close(); + + // create a copy of this class to be written to output file + AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis"); + + // add PDG codes to decay tracks in found candidates (in simulation mode) + // and store tree in the output file + if(!fSim) { + TFile *outroot = new TFile(outName1.Data(),"recreate"); + treeD0->Write(); + copy->Write(); + outroot->Close(); + delete outroot; + } else { + printf(" Now adding information from simulation (PDG codes) ...\n"); + TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates"); + SimulationInfo(treeD0,treeD0sim); + delete treeD0; + TFile *outroot = new TFile(outName1.Data(),"recreate"); + treeD0sim->Write(); + copy->Write(); + outroot->Close(); + delete outroot; + } + + // write to a file the total number of events + FILE *outdat = fopen(outName2.Data(),"w"); + fprintf(outdat,"%d\n",nTotEv); + fclose(outdat); + + return; +} +//----------------------------------------------------------------------------- +void AliD0toKpiAnalysis::PrintStatus() const { + // Print parameters being used + + printf("Preselections:\n"); + printf(" fPtCut = %f GeV\n",fPtCut); + printf(" fd0Cut = %f micron\n",fd0Cut); + printf(" fMassCut = %f GeV\n",fMassCut); + if(fVertexOnTheFly) printf("Primary vertex on the fly\n"); + if(fSim) { + printf("Simulation mode\n"); + if(fOnlySignal) printf(" Only signal goes to file\n"); + } + printf("Cuts on candidates:\n"); + printf(" |M-MD0| [GeV] < %f\n",fD0Cuts[0]); + printf(" dca [micron] < %f\n",fD0Cuts[1]); + printf(" cosThetaStar < %f\n",fD0Cuts[2]); + printf(" pTK [GeV] > %f\n",fD0Cuts[3]); + printf(" pTpi [GeV] > %f\n",fD0Cuts[4]); + printf(" |d0K| [micron] < %f\n",fD0Cuts[5]); + printf(" |d0pi| [micron] < %f\n",fD0Cuts[6]); + printf(" d0d0 [micron^2] < %f\n",fD0Cuts[7]); + printf(" cosThetaPoint > %f\n",fD0Cuts[8]); + + return; +} +//----------------------------------------------------------------------------- +Bool_t AliD0toKpiAnalysis::SelectInvMass(const Double_t p[6]) const { + // Apply preselection in the invariant mass of the pair + + Double_t mD0 = 1.8645; + Double_t mPi = 0.13957; + Double_t mKa = 0.49368; + + Double_t energy[2]; + Double_t mom2[2],momTot2; + + mom2[0] = p[0]*p[0] + p[1]*p[1] + p[2]*p[2]; + mom2[1] = p[3]*p[3] + p[4]*p[4] + p[5]*p[5]; + + momTot2 = (p[0]+p[3])*(p[0]+p[3])+ + (p[1]+p[4])*(p[1]+p[4])+ + (p[2]+p[5])*(p[2]+p[5]); + + // D0 -> K- pi+ + energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]); + energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]); + + Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2); + + // D0bar -> K+ pi- + energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]); + energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]); + + Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2); + + if(TMath::Abs(minvD0-mD0) < fMassCut) return kTRUE; + if(TMath::Abs(minvD0bar-mD0) < fMassCut) return kTRUE; + return kFALSE; +} +//----------------------------------------------------------------------------- +void AliD0toKpiAnalysis::SelectTracks(AliESD *event, + TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP, + TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const { + // Create two TObjArrays with positive and negative tracks and + // apply single-track preselection + + nTrksP=0,nTrksN=0; + + Int_t entr = event->GetNumberOfTracks(); + + // transfer ITS tracks from ESD to arrays and to a tree + for(Int_t i=0; iGetTrack(i); + UInt_t status = esdtrack->GetStatus(); + + if(!(status&AliESDtrack::kITSin)) continue; + + // single track selection + if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue; + + if(esdtrack->GetSign()<0) { // negative track + trksN.AddLast(esdtrack); + trkEntryN[nTrksN] = i; + nTrksN++; + } else { // positive track + trksP.AddLast(esdtrack); + trkEntryP[nTrksP] = i; + nTrksP++; + } + + } // loop on ESD tracks + + return; +} +//----------------------------------------------------------------------------- +void AliD0toKpiAnalysis::SetD0Cuts(Double_t cut0,Double_t cut1, + Double_t cut2,Double_t cut3,Double_t cut4, + Double_t cut5,Double_t cut6, + Double_t cut7,Double_t cut8) { + // Set the cuts for D0 selection + fD0Cuts[0] = cut0; + fD0Cuts[1] = cut1; + fD0Cuts[2] = cut2; + fD0Cuts[3] = cut3; + fD0Cuts[4] = cut4; + fD0Cuts[5] = cut5; + fD0Cuts[6] = cut6; + fD0Cuts[7] = cut7; + fD0Cuts[8] = cut8; + + return; +} +//----------------------------------------------------------------------------- +void AliD0toKpiAnalysis::SetD0Cuts(const Double_t cuts[9]) { + // Set the cuts for D0 selection + + for(Int_t i=0; i<9; i++) fD0Cuts[i] = cuts[i]; + + return; +} +//----------------------------------------------------------------------------- +Bool_t +AliD0toKpiAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const { + // Check if track passes some kinematical cuts + // Magnetic field "b" (kG) + + if(TMath::Abs(1./trk.GetParameter()[4]) < fPtCut) + return kFALSE; + if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1],b)) < fd0Cut) + return kFALSE; + + return kTRUE; +} +//---------------------------------------------------------------------------- +void AliD0toKpiAnalysis::MakeTracksRefFile(AliRun *gAlice, + Int_t evFirst,Int_t evLast) const { + // Create a file with simulation info for the reconstructed tracks + + TFile *outFile = TFile::Open("D0TracksRefFile.root","recreate"); + TFile *esdFile = TFile::Open("AliESDs.root"); + + AliMC *mc = gAlice->GetMCApp(); + + Int_t label; + TParticle *part; + TParticle *mumpart; + REFTRACK reftrk; + + AliESD* event = new AliESD; + TTree* tree = (TTree*) esdFile->Get("esdTree"); + tree->SetBranchAddress("ESD",&event); + // loop on events in file + for(Int_t iEvent=evFirst; iEventGetEntries(); iEvent++) { + if(iEvent>evLast) break; + tree->GetEvent(iEvent); + Int_t ev = (Int_t)event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number. + + gAlice->GetEvent(ev); + + Int_t nentr=(Int_t)event->GetNumberOfTracks(); + + // Tree for true track parameters + char ttname[100]; + sprintf(ttname,"Tree_Ref_%d",ev); + TTree *reftree = new TTree(ttname,"Tree with true track params"); + reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz"); + + for(Int_t i=0; iGetTrack(i); + label = TMath::Abs(esdtrack->GetLabel()); + + part = (TParticle*)mc->Particle(label); + reftrk.lab = label; + reftrk.pdg = part->GetPdgCode(); + reftrk.mumlab = part->GetFirstMother(); + if(part->GetFirstMother()>=0) { + mumpart = (TParticle*)gAlice->GetMCApp()->Particle(part->GetFirstMother()); + reftrk.mumpdg = mumpart->GetPdgCode(); + reftrk.mumprongs = mumpart->GetNDaughters(); + } else { + reftrk.mumpdg=-1; + reftrk.mumprongs=-1; + } + reftrk.Vx = part->Vx(); + reftrk.Vy = part->Vy(); + reftrk.Vz = part->Vz(); + reftrk.Px = part->Px(); + reftrk.Py = part->Py(); + reftrk.Pz = part->Pz(); + + reftree->Fill(); + + } // loop on tracks + + outFile->cd(); + reftree->Write(); + + delete reftree; + } // loop on events + + esdFile->Close(); + outFile->Close(); + + return; +} +//----------------------------------------------------------------------------- +void AliD0toKpiAnalysis::SimulationInfo(TTree *treeD0in,TTree *treeD0out) const { + // add pdg codes to candidate decay tracks (for sim) + + TString refFileName("D0TracksRefFile.root"); + if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) { + printf("AliD0toKpiAnalysis::SimulationInfo: no reference file found!\n"); + return; + } + TFile *refFile = TFile::Open(refFileName.Data()); + + Char_t refTreeName[100]; + Int_t event; + Int_t pdg[2],mumpdg[2],mumlab[2]; + REFTRACK reftrk; + + // read-in reference tree for event 0 (the only event for Pb-Pb) + sprintf(refTreeName,"Tree_Ref_%d",0); + TTree *refTree0 = (TTree*)refFile->Get(refTreeName); + refTree0->SetBranchAddress("rectracks",&reftrk); + + AliD0toKpi *theD0 = 0; + treeD0in->SetBranchAddress("D0toKpi",&theD0); + treeD0out->Branch("D0toKpi","AliD0toKpi",&theD0,200000,0); + + Int_t entries = (Int_t)treeD0in->GetEntries(); + + for(Int_t i=0; iGetEvent(i); + event = theD0->EventNo(); + + if(event==0) { // always true for Pb-Pb (avoid to read-in tree every time) + refTree0->GetEvent(theD0->GetTrkNum(0)); + pdg[0] = reftrk.pdg; + mumpdg[0] = reftrk.mumpdg; + mumlab[0] = reftrk.mumlab; + refTree0->GetEvent(theD0->GetTrkNum(1)); + pdg[1] = reftrk.pdg; + mumpdg[1] = reftrk.mumpdg; + mumlab[1] = reftrk.mumlab; + } else { + sprintf(refTreeName,"Tree_Ref_%d",event); + TTree *refTree = (TTree*)refFile->Get(refTreeName); + refTree->SetBranchAddress("rectracks",&reftrk); + refTree->GetEvent(theD0->GetTrkNum(0)); + pdg[0] = reftrk.pdg; + mumpdg[0] = reftrk.mumpdg; + mumlab[0] = reftrk.mumlab; + refTree->GetEvent(theD0->GetTrkNum(1)); + pdg[1] = reftrk.pdg; + mumpdg[1] = reftrk.mumpdg; + mumlab[1] = reftrk.mumlab; + delete refTree; + } + + theD0->SetPdgCodes(pdg); + theD0->SetMumPdgCodes(mumpdg); + + if(TMath::Abs(mumpdg[0])==421 && + TMath::Abs(mumpdg[1])==421 && + mumlab[0]==mumlab[1] && + reftrk.mumprongs==2 && + ((TMath::Abs(pdg[0])==211 && TMath::Abs(pdg[1])==321) || + (TMath::Abs(pdg[0])==321 && TMath::Abs(pdg[1])==211)) + ) theD0->SetSignal(); + + if(!fOnlySignal || theD0->IsSignal()) treeD0out->Fill(); + + } + + delete refTree0; + + refFile->Close(); + + return; +} + + + + + + diff --git a/PWG3/AliD0toKpiAnalysis.h b/PWG3/AliD0toKpiAnalysis.h new file mode 100644 index 00000000000..bc6ebcff56b --- /dev/null +++ b/PWG3/AliD0toKpiAnalysis.h @@ -0,0 +1,90 @@ +#ifndef AliD0toKpiAnalysis_H +#define AliD0toKpiAnalysis_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//------------------------------------------------------------------------- +// Class AliD0toKpiAnalysis +// Reconstruction and analysis D0 -> K^- pi^+ +// +// Origin: A. Dainese andrea.dainese@lnl.infn.it +//------------------------------------------------------------------------- + +#include +#include +#include "AliESD.h" +#include "AliRun.h" + +//----------------------------------------------------------------------------- +class AliD0toKpiAnalysis : public TNamed { + public: + // + AliD0toKpiAnalysis(); + virtual ~AliD0toKpiAnalysis(); + + void ApplySelection(const Char_t *inName="AliD0toKpi.root", + const Char_t *outName="AliD0toKpi_sele.root") const; + void FindCandidates(Int_t evFirst=0,Int_t evLast=0, + const Char_t *outName="AliD0toKpi.root"); + void MakeTracksRefFile(AliRun *gAlice,Int_t evFirst=0,Int_t evLast=0) const; + void PrintStatus() const; + void SetVertexOnTheFly() { fVertexOnTheFly=kTRUE; } + void SetSimulation() { fSim=kTRUE; } + void SetOnlySignal() { fOnlySignal=kTRUE; } + void SetPtCut(Double_t pt=0.) { fPtCut=pt; } + void Setd0Cut(Double_t d0=0.) { fd0Cut=d0; } + void SetMassCut(Double_t deltaM=1000.) { fMassCut=deltaM; } + void SetD0Cuts(Double_t cut0=1000.,Double_t cut1=100000., + Double_t cut2=1.1,Double_t cut3=0.,Double_t cut4=0., + Double_t cut5=100000.,Double_t cut6=100000., + Double_t cut7=100000000.,Double_t cut8=-1.1); + void SetD0Cuts(const Double_t cuts[9]); + void SetPID(const Char_t * pid="TOFparam_PbPb") { fPID=pid; } + // + private: + // + Bool_t fVertexOnTheFly; // flag for primary vertex reco on the fly + Bool_t fSim; // flag for the analysis of simulated events + Bool_t fOnlySignal; // write to file only signal candidates (for sim) + TString fPID; // PID scheme + + Double_t fV1[3]; // primary vertex position (in cm) + Double_t fPtCut; // minimum track pt (in GeV/c) + Double_t fd0Cut; // minimum track |rphi impact parameter| (in micron) + Double_t fMassCut; // maximum of |InvMass-MD0| (in GeV) + Double_t fD0Cuts[9]; // cuts on D0 candidates (see SetD0Cuts()) + // (to be passed to function AliD0toKpi::Select()) + // 0 = inv. mass half width [GeV] + // 1 = dca [micron] + // 2 = cosThetaStar + // 3 = pTK [GeV/c] + // 4 = pTPi [GeV/c] + // 5 = d0K [micron] upper limit! + // 6 = d0Pi [micron] upper limit! + // 7 = d0d0 [micron^2] + // 8 = cosThetaPoint + + // + Double_t CalculateTOFmass(Double_t mom,Double_t length,Double_t time) const; + Bool_t SelectInvMass(const Double_t p[6]) const; + void SelectTracks(AliESD *event, + TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP, + TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const; + void SetVertex1(Double_t x=0.,Double_t y=0.,Double_t z=0.) + { fV1[0]=x;fV1[1]=y;fV1[2]=z; } + void SimulationInfo(TTree *treeD0in,TTree *treeD0out) const; + Bool_t SingleTrkCuts(const AliESDtrack& trk, Double_t b) const; + // + ClassDef(AliD0toKpiAnalysis,3) // Reconstruction of D0 candidates class +}; + + +#endif + + + + + + + + diff --git a/PWG3/AliD0toKpiPlots.C b/PWG3/AliD0toKpiPlots.C new file mode 100644 index 00000000000..2096a59c1ef --- /dev/null +++ b/PWG3/AliD0toKpiPlots.C @@ -0,0 +1,235 @@ +void AliD0toKpiPlots(const Char_t *inName="AliD0toKpi.root", + const Char_t *outName="D0histograms.root") { + //-------------------------------------------------------------------------- + // This macro histograms many variables of D0->Kpi candidates + // + // Andrea Dainese, andrea.dainese@lnl.infn.it + //-------------------------------------------------------------------------- + + gSystem->Load("libAOD.so"); + gSystem->Load("libPWG3base.so"); + + // set of cuts + Double_t D0Cuts[9] = {0.1, // mass [GeV] + 1000000., // dca [micron] + 1.1, // cosThetaStar + 0., // pT K [GeV/c] + 0., // pT Pi [GeV/c] + 100000., // d0K upper [micron] + 100000., // d0Pi upper [micron] + 10000000000., // d0d0 [micron^2] + -1.1}; // cosThetaPointing + + // number of events (for normalization) + Bool_t normalize = kFALSE; + Double_t events = 1.; + + + + // define histograms + TH1F *hptK = new TH1F("hptK","\"K\" p_{t} distribution",50,0,10); + hptK->SetXTitle("p_{t} [GeV]"); + + TH1F *hptPi = new TH1F("hptPi","\"#pi\" p_{t} distribution",50,0,10); + hptPi->SetXTitle("p_{t} [GeV]"); + + TH1F *hDCA = new TH1F("hDCA","DCA",50,0,1000); + hDCA->SetXTitle("dca [#mu m]"); + + TH1F *hptD0 = new TH1F("hptD0","D^{0} p_{t} distribution",40,0,40); + hptD0->SetXTitle("p_{t} [GeV]"); + + TH1F *hyD0 = new TH1F("hyD0","D^{0} rapidity distribution",50,-2,2); + hyD0->SetXTitle("y"); + + TH1F *hCPtaD0 = new TH1F("hCPtaD0","cosine of pointing angle distribution",100,-1,1); + hCPtaD0->SetXTitle("cos #theta_{point}"); + + TH1F *hCPtaXY = new TH1F("hCPtaXY","cosine of pointing angle in (x,y) plane",100,-1,1); + hCPtaXY->SetXTitle("cos #theta_{point}"); + + TH1F *hCts = new TH1F("hCts","cosine of decay angle",50,-1.2,1.2); + hCts->SetXTitle("cos #theta^{*}"); + + TH2F *hCtsVsPtK = new TH2F("hCtsVsPtK","cosine of decay angle VS \"K\" p_{t}",50,0,5,50,-1,1); + hCtsVsPtK->SetYTitle("cos #theta^{*}"); + hCtsVsPtK->SetXTitle("p_{t} [GeV]"); + + TH1F *hd0d0 = new TH1F("hd0d0","Product of the impact parameters",100,-100000,100000); + hd0d0->SetXTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]"); + + TH1F *hd0K = new TH1F("hd0K","Impact parameter of \"K\"",100,-5000,5000); + hd0K->SetXTitle("d_{0}^{K} [#mu m]"); + + TH1F *hd0Pi = new TH1F("hd0Pi","Impact parameter of \"#pi\"",100,-5000,5000); + hd0Pi->SetXTitle("d_{0}^{#pi} [#mu m]"); + + TH2F *hCPtaVsd0d0 = new TH2F("hCPtaVsd0d0","cos #theta_{point} vs d_{0}^{K} #times d_{0}^{#pi}",100,-100000,100000,100,-1,1); + hCPtaVsd0d0->SetXTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]"); + hCPtaVsd0d0->SetYTitle("cos #theta_{point}"); + + TH2F *hCPtaVsd0d0zoom = new TH2F("hCPtaVsd0d0zoom","cos #theta_{point} vs d_{0}^{K} #times d_{0}^{#pi}",100,-100000,0,100,.9,1); + hCPtaVsd0d0zoom->SetXTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]"); + hCPtaVsd0d0zoom->SetYTitle("cos #theta_{point}"); + + TH2F *hd0d0VSptD0 = new TH2F("hd0d0VSptD0","d_{0}^{K} #times d_{0}^{#pi} VS D^{0} p_{t}",50,0,25,100,-120000,120000); + hd0d0VSptD0->SetYTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]"); + hd0d0VSptD0->SetXTitle("D^{0} p_{t} [GeV]"); + + TH1F *hMass = new TH1F("hMass","Invariant mass distribution",50,1.765,1.965); + hMass->SetXTitle("M[K,#pi] [GeV]"); + + TH2F *hArm = new TH2F("hArm","Armenteros plot",50,-2,2,50,0,1); + hArm->SetXTitle("#alpha"); + hArm->SetYTitle("q_{t}"); + + // open input file and get tree + TFile *inFile = TFile::Open(inName); + + TTree *treeD0 = (TTree*)inFile->Get("TreeD0"); + AliD0toKpi *D = 0; + treeD0->SetBranchAddress("D0toKpi",&D); + Int_t entries = (Int_t)treeD0->GetEntries(); + + printf("+++\n+++ Number of D0 in tree: %d\n+++\n",entries); + + Double_t MD0,MD0bar,ctsD0,ctsD0bar,ctsPiD0,ctsPiD0bar; + Double_t WgtD0,WgtD0bar; + Double_t sampleABC=0.; + Int_t okD0=0,okD0bar=0; + Int_t nSel = 0; + Int_t ptbin; + + // loop on D0 + for(Int_t i=0; iGetEvent(i); + //--- select the PID strategy & compute weights + // D->ApplyPID("TOFparam_PbPb"); + // D->ComputeWgts(); + // get weights for the three samples A+B+C + // D->GetWgts(WgtD0,WgtD0bar,"ABC"); + WgtD0 = 1.; WgtD0bar = 1.; + + // normalize to 1 event + if(normalize) { WgtD0 /= events; WgtD0bar /= events; } + + // check if candidate passes selection (as D0 or D0bar) + D->Select(D0Cuts,okD0,okD0bar); + + // set weights to 0 if the candidate doesn't pass selection + if(!okD0) WgtD0=0.; + if(!okD0bar) WgtD0bar=0.; + if(okD0 || okD0bar) nSel++; + + // count selected candidates + sampleABC += WgtD0 + WgtD0bar; + + // inv mass and cosThetaStar + D->InvMass(MD0,MD0bar); + D->CosThetaStar(ctsD0,ctsD0bar); + + // fill histograms + hptK->Fill(D->PtChild(1),WgtD0); + hptK->Fill(D->PtChild(0),WgtD0bar); + hptPi->Fill(D->PtChild(0),WgtD0); + hptPi->Fill(D->PtChild(1),WgtD0bar); + hd0K->Fill(D->Getd0Child(1),WgtD0); + hd0K->Fill(D->Getd0Child(0),WgtD0bar); + hd0Pi->Fill(D->Getd0Child(0),WgtD0); + hd0Pi->Fill(D->Getd0Child(1),WgtD0bar); + hMass->Fill(MD0,WgtD0); + hMass->Fill(MD0bar,WgtD0bar); + hCts->Fill(ctsD0,WgtD0); + hCts->Fill(ctsD0bar,WgtD0bar); + hCtsVsPtK->Fill(D->PtChild(1),ctsD0,WgtD0); + hCtsVsPtK->Fill(D->PtChild(0),ctsD0bar,WgtD0bar); + hDCA->Fill(D->GetDCA(),WgtD0+WgtD0bar); + hptD0->Fill(D->Pt(),WgtD0+WgtD0bar); + hyD0->Fill(D->Rapidity(),WgtD0+WgtD0bar); + hd0d0->Fill(D->ProdImpParams(),WgtD0+WgtD0bar); + hCPtaD0->Fill(D->CosPointing(),WgtD0+WgtD0bar); + hCPtaXY->Fill(D->CosPointingXY(),WgtD0+WgtD0bar); + hCPtaVsd0d0->Fill(D->ProdImpParams(),D->CosPointing(),WgtD0+WgtD0bar); + hd0d0VSptD0->Fill(D->Pt(),D->ProdImpParams(),WgtD0+WgtD0bar); + hCPtaVsd0d0zoom->Fill(D->ProdImpParams(),D->CosPointing(),WgtD0+WgtD0bar); + hArm->Fill(D->Alpha(),D->Qt(),WgtD0+WgtD0bar); + + + } // end loop on D0 candidates + + inFile->Close(); + + printf("\n\n --- Total number of candidates passing selection: %d\n\n --- Sum of weights sample A+B+C: %f\n\n",nSel,sampleABC); + + // draw histograms + TCanvas *c1 = new TCanvas("c1","pt K & pi",0,0,700,700); + c1->SetLogy(); + hptK->Draw(); + hptPi->Draw("same"); + + TCanvas *c2 = new TCanvas("c2","pt D0",0,0,700,700); + c2->SetLogy(); + hptD0->Draw(); + + TCanvas *c3 = new TCanvas("c3","rapidity D0",0,0,700,700); + hyD0->Draw(); + + TCanvas *c4 = new TCanvas("c4","pointing angle",0,0,700,700); + hCPtaD0->Draw(); + + TCanvas *c5 = new TCanvas("c5","d0 x d0",0,0,700,700); + c5->SetLogy(); + hd0d0->Draw(); + + TCanvas *c6 = new TCanvas("c6","pointing angle VS d0d0",0,0,700,700); + c6->SetLogz(); + hCPtaVsd0d0->Draw("box"); + + TCanvas *c7 = new TCanvas("c7","mass",0,0,700,700); + hMass->Draw(); + + TCanvas *c8 = new TCanvas("c8","armenteros",0,0,700,700); + hArm->Draw("box"); + + TCanvas *c9 = new TCanvas("c9","decay angle",0,0,700,700); + hCts->Draw(); + + TCanvas *c10 = new TCanvas("c10","dca",0,0,700,700); + c10->SetLogy(); + hDCA->Draw(); + + TCanvas *c11 = new TCanvas("c11","d0 K & pi",0,0,700,700); + c11->SetLogy(); + hd0K->Draw(); + hd0Pi->Draw("same"); + + // write all histograms to file + TFile *outFile = new TFile(outName,"recreate"); + hMass->Write(); + hDCA->Write(); + hCts->Write(); + hCtsVsPtK->Write(); + hArm->Write(); + hCPtaVsd0d0->Write(); + hd0d0VSptD0->Write(); + hCPtaVsd0d0zoom->Write(); + hd0d0->Write(); + hCPtaD0->Write(); + hCPtaXY->Write(); + hptK->Write(); + hptPi->Write(); + hptD0->Write(); + hyD0->Write(); + hd0K->Write(); + hd0Pi->Write(); + outFile->Close(); + + return; +} + + + + diff --git a/PWG3/AliD0toKpiTest.C b/PWG3/AliD0toKpiTest.C index 83a93f3c261..5ddb852aeb8 100644 --- a/PWG3/AliD0toKpiTest.C +++ b/PWG3/AliD0toKpiTest.C @@ -6,7 +6,8 @@ void AliD0toKpiReco() { - gSystem->Load("libANALYSIS.so"); + gSystem->Load("libAOD.so"); + gSystem->Load("libPWG3base.so"); //============== R E C O N S T R U C T I O N ============================== @@ -67,7 +68,8 @@ void AliD0toKpiReco() { //========================================================================== void AliD0toKpiSele() { - gSystem->Load("libANALYSIS.so"); + gSystem->Load("libAOD.so"); + gSystem->Load("libPWG3base.so"); //======================== S E L E C T I O N ============================ -- 2.39.3