From: hristov Date: Mon, 16 Feb 2004 15:52:15 +0000 (+0000) Subject: New analysis code (A.Dainese) X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=3a9a34872cda8f00f7ad553db52916d4f36c6dcc;hp=4133c0aaa75f1d4539fd64d981748ef73cbf8871;ds=sidebyside New analysis code (A.Dainese) --- diff --git a/ANALYSIS/ANALYSISLinkDef.h b/ANALYSIS/ANALYSISLinkDef.h new file mode 100644 index 00000000000..b01effe4aae --- /dev/null +++ b/ANALYSIS/ANALYSISLinkDef.h @@ -0,0 +1,11 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class AliD0toKpi+; +#pragma link C++ class AliD0toKpiAnalysis+; + + +#endif diff --git a/ANALYSIS/AliD0toKpi.cxx b/ANALYSIS/AliD0toKpi.cxx new file mode 100644 index 00000000000..4ddcb216096 --- /dev/null +++ b/ANALYSIS/AliD0toKpi.cxx @@ -0,0 +1,677 @@ +/************************************************************************** + * 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 +// +// Note: the two decay tracks are labelled: 0 (positive track) +// 1 (negative track) +// +// Origin: A. Dainese andrea.dainese@pd.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) { + + const char *tofparampbpb = strstr(pidScheme.Data(),"TOFparam_PbPb"); + const char *tofparampp = strstr(pidScheme.Data(),"TOFparam_pp"); + + if((tofparampbpb || tofparampp) && fPdg[0]==0) { + printf("AliD0toKpi::ApplyPID :\n Warning: TOF parameterized PID can be used only for simulation!\n"); + return; + } + + if(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),kPiBins_PbPb,kPiBinWidth_PbPb,kPiTagPi_PbPb); + fTagNid[0] = 1.-fTagPi[0]; + fTagKa[0] = 0.; + fTagPr[0] = 0.; + } + if(TMath::Abs(fPdg[0])==321) { // kaon + fTagKa[0] = LinearInterpolation(PChild(0),kKBins_PbPb,kKBinWidth_PbPb,kKTagK_PbPb); + fTagNid[0] = LinearInterpolation(PChild(0),kKBins_PbPb,kKBinWidth_PbPb,kKTagNid_PbPb); + 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),kPBins_PbPb,kPBinWidth_PbPb,kPTagP_PbPb); + fTagNid[0] = LinearInterpolation(PChild(0),kPBins_PbPb,kPBinWidth_PbPb,kPTagNid_PbPb); + 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),kPiBins_PbPb,kPiBinWidth_PbPb,kPiTagPi_PbPb); + fTagNid[1] = 1.-fTagPi[1]; + fTagKa[1] = 0.; + fTagPr[1] = 0.; + } + if(TMath::Abs(fPdg[1])==321) { // kaon + fTagKa[1] = LinearInterpolation(PChild(1),kKBins_PbPb,kKBinWidth_PbPb,kKTagK_PbPb); + fTagNid[1] = LinearInterpolation(PChild(1),kKBins_PbPb,kKBinWidth_PbPb,kKTagNid_PbPb); + 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),kPBins_PbPb,kPBinWidth_PbPb,kPTagP_PbPb); + fTagNid[1] = LinearInterpolation(PChild(1),kPBins_PbPb,kPBinWidth_PbPb,kPTagNid_PbPb); + if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1]; + fTagPi[1] = 1.-fTagNid[1]-fTagPr[1]; + fTagKa[1] = 0.; + } + } + + + if(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),kPiBins_pp,kPiBinWidth_pp,kPiTagPi_pp); + fTagNid[0] = 1.-fTagPi[0]; + fTagKa[0] = 0.; + fTagPr[0] = 0.; + } + if(TMath::Abs(fPdg[0])==321) { // kaon + fTagKa[0] = LinearInterpolation(PChild(0),kKBins_pp,kKBinWidth_pp,kKTagK_pp); + fTagNid[0] = LinearInterpolation(PChild(0),kKBins_pp,kKBinWidth_pp,kKTagNid_pp); + 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),kPBins_pp,kPBinWidth_pp,kPTagP_pp); + fTagNid[0] = LinearInterpolation(PChild(0),kPBins_pp,kPBinWidth_pp,kPTagNid_pp); + 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),kPiBins_pp,kPiBinWidth_pp,kPiTagPi_pp); + fTagNid[1] = 1.-fTagPi[1]; + fTagKa[1] = 0.; + fTagPr[1] = 0.; + } + if(TMath::Abs(fPdg[1])==321) { // kaon + fTagKa[1] = LinearInterpolation(PChild(1),kKBins_pp,kKBinWidth_pp,kKTagK_pp); + fTagNid[1] = LinearInterpolation(PChild(1),kKBins_pp,kKBinWidth_pp,kKTagNid_pp); + 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),kPBins_pp,kPBinWidth_pp,kPTagP_pp); + fTagNid[1] = LinearInterpolation(PChild(1),kPBins_pp,kPBinWidth_pp,kPTagNid_pp); + 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 "< 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; +} +//----------------------------------------------------------------------------- +void AliD0toKpi::DrawPIDinTOF(TString pidScheme) const { + // Draw parameterized PID probabilities in TOF + + const char *tofparampbpb = strstr(pidScheme.Data(),"TOFparam_PbPb"); + const char *tofparampp = strstr(pidScheme.Data(),"TOFparam_pp"); + + TH2F* framePi = new TH2F("framePi","Tag probabilities for PIONS",2,0,2.5,2,0,1); + framePi->SetXTitle("p [GeV/c]"); + framePi->SetStats(0); + TH2F* frameK = new TH2F("frameK","Tag probabilities for KAONS",2,0,2.5,2,0,1); + frameK->SetXTitle("p [GeV/c]"); + frameK->SetStats(0); + TH2F* frameP = new TH2F("frameP","Tag probabilities for PROTONS",2,0,4.5,2,0,1); + frameP->SetXTitle("p [GeV/c]"); + frameP->SetStats(0); + + TH1F* hPiPi = new TH1F("hPiPi","Tag probabilities for PIONS",kPiBins_PbPb,0,2.5); + TH1F* hPiNid = new TH1F("hPiNid","Tag probabilities for PIONS",kPiBins_PbPb,0,2.5); + + TH1F* hKK = new TH1F("hKK","Tag probabilities for KAONS",kKBins_PbPb,0,2.5); + TH1F* hKNid = new TH1F("hKNid","Tag probabilities for KAONS",kKBins_PbPb,0,2.5); + TH1F* hKPi = new TH1F("hKPi","Tag probabilities for KAONS",kKBins_PbPb,0,2.5); + + TH1F* hPP = new TH1F("hPP","Tag probabilities for PROTONS",kPBins_PbPb,0,4.5); + TH1F* hPNid = new TH1F("hPNid","Tag probabilities for PROTONS",kPBins_PbPb,0,4.5); + TH1F* hPPi = new TH1F("hPPi","Tag probabilities for PROTONS",kPBins_PbPb,0,4.5); + + + if(tofparampbpb) { + + for(Int_t i=1; i<=kPiBins_PbPb; i++) { + hPiPi->SetBinContent(i,kPiTagPi_PbPb[i-1]); + hPiNid->SetBinContent(i,kPiTagPi_PbPb[i-1]+kPiTagNid_PbPb[i-1]); + + hKK->SetBinContent(i,kKTagK_PbPb[i-1]); + hKPi->SetBinContent(i,kKTagK_PbPb[i-1]+kKTagPi_PbPb[i-1]); + hKNid->SetBinContent(i,kKTagK_PbPb[i-1]+kKTagPi_PbPb[i-1]+kKTagNid_PbPb[i-1]); + } + for(Int_t i=1; i<=kPBins_PbPb; i++) { + hPP->SetBinContent(i,kPTagP_PbPb[i-1]); + hPPi->SetBinContent(i,kPTagP_PbPb[i-1]+kPTagPi_PbPb[i-1]); + hPNid->SetBinContent(i,kPTagP_PbPb[i-1]+kPTagPi_PbPb[i-1]+kPTagNid_PbPb[i-1]); + } + + } else if(tofparampp) { + + for(Int_t i=1; i<=kPiBins_pp; i++) { + hPiPi->SetBinContent(i,kPiTagPi_pp[i-1]); + hPiNid->SetBinContent(i,kPiTagPi_pp[i-1]+kPiTagNid_pp[i-1]); + + hKK->SetBinContent(i,kKTagK_pp[i-1]); + hKPi->SetBinContent(i,kKTagK_pp[i-1]+kKTagPi_pp[i-1]); + hKNid->SetBinContent(i,kKTagK_pp[i-1]+kKTagPi_pp[i-1]+kKTagNid_pp[i-1]); + } + for(Int_t i=1; i<=kPBins_pp; i++) { + hPP->SetBinContent(i,kPTagP_pp[i-1]); + hPPi->SetBinContent(i,kPTagP_pp[i-1]+kPTagPi_pp[i-1]); + hPNid->SetBinContent(i,kPTagP_pp[i-1]+kPTagPi_pp[i-1]+kPTagNid_pp[i-1]); + } + + } + + + TCanvas* c = new TCanvas("c","Parameterized PID in TOF",0,0,1000,400); + c->Divide(3,1); + c->cd(1); + framePi->Draw(); + hPiNid->SetFillColor(18); hPiNid->Draw("same"); + hPiPi->SetFillColor(4); hPiPi->Draw("same"); + TPaveLabel* pav1 = new TPaveLabel(1,.2,1.4,.3,"#pi"); + pav1->SetBorderSize(0); + pav1->Draw("same"); + TPaveLabel* pav2 = new TPaveLabel(1,.8,1.8,.9,"non-id"); + pav2->SetBorderSize(0); + pav2->Draw("same"); + + c->cd(2); + frameK->Draw(); + hKNid->SetFillColor(18); hKNid->Draw("same"); + hKPi->SetFillColor(4); hKPi->Draw("same"); + hKK->SetFillColor(7); hKK->Draw("same"); + TPaveLabel* pav3 = new TPaveLabel(1,.2,1.5,.3,"K"); + pav3->SetBorderSize(0); + pav3->Draw("same"); + TPaveLabel* pav4 = new TPaveLabel(1,.8,1.8,.9,"non-id"); + pav4->SetBorderSize(0); + pav4->Draw("same"); + TPaveLabel* pav5 = new TPaveLabel(.4,.5,.8,.6,"#pi"); + pav5->SetBorderSize(0); + pav5->Draw("same"); + + c->cd(3); + frameP->Draw(); + hPNid->SetFillColor(18); hPNid->Draw("same"); + hPPi->SetFillColor(4); hPPi->Draw("same"); + hPP->SetFillColor(3); hPP->Draw("same"); + TPaveLabel* pav6 = new TPaveLabel(1,.2,1.5,.3,"p"); + pav6->SetBorderSize(0); + pav6->Draw("same"); + TPaveLabel* pav7 = new TPaveLabel(1,.8,2.6,.9,"non-id"); + pav7->SetBorderSize(0); + pav7->Draw("same"); + TPaveLabel* pav8 = new TPaveLabel(.2,.5,1,.6,"#pi"); + pav8->SetBorderSize(0); + pav8->Draw("same"); + + + 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; +} +//---------------------------------------------------------------------------- + + + + + + +/* +//____________________________________________________________________________ +void AliD0toKpi::SetPtWgts4pp() { + // Correct pt distribution in order to reproduce MNR pt slope + // (for pp generated with PYTHIA min. bias) + + if(TMath::Abs(fMum[0]) != 421 && TMath::Abs(fMum[1]) != 421 && + TMath::Abs(fMum[0]) != 411 && TMath::Abs(fMum[1]) != 411) return; + + Double_t ptWgt = 1.; + ptWgt = 2.05-0.47*Pt()+0.02*Pt()*Pt(); + if(Pt() >= 5.) ptWgt = 0.56*TMath::Exp(-0.12*Pt()); + + fWgtAD0 *= ptWgt; + fWgtAD0bar *= ptWgt; + fWgtBD0 *= ptWgt; + fWgtBD0bar *= ptWgt; + fWgtCD0 *= ptWgt; + fWgtCD0bar *= ptWgt; + fWgtDD0 *= ptWgt; + fWgtDD0bar *= ptWgt; + + return; +} +//____________________________________________________________________________ +*/ + + + + + + + + + diff --git a/ANALYSIS/AliD0toKpi.h b/ANALYSIS/AliD0toKpi.h new file mode 100644 index 00000000000..9448cb59cb1 --- /dev/null +++ b/ANALYSIS/AliD0toKpi.h @@ -0,0 +1,198 @@ +#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 +// +// Origin: A. Dainese andrea.dainese@pd.infn.it +//------------------------------------------------------------------------- + +#include +#include +#include +#include +#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 kPiBins_PbPb = 10; +const Double_t kPiBinWidth_PbPb = 0.250; +const Double_t kPiTagPi_PbPb[kPiBins_PbPb] = {0.211421,0.652184,0.624421,0.614727,0.610777,0.628015,0.631520,0.630324,0.637551,0.575235}; +const Double_t kPiTagNid_PbPb[kPiBins_PbPb] = {0.788579,0.347816,0.375579,0.385273,0.389223,0.371985,0.368480,0.369676,0.362449,0.424765}; +// KAONS +const Int_t kKBins_PbPb = 10; +const Double_t kKBinWidth_PbPb = 0.250; +const Double_t kKTagK_PbPb[kKBins_PbPb] = {0.000000,0.101255,0.397662,0.467586,0.517008,0.555023,0.584185,0.519029,0.464117,0.247308}; +const Double_t kKTagPi_PbPb[kKBins_PbPb] = {0.102049,0.289930,0.101930,0.057771,0.040286,0.028567,0.053108,0.094369,0.066302,0.247308}; +const Double_t kKTagNid_PbPb[kKBins_PbPb] = {0.897951,0.608815,0.500408,0.474643,0.442705,0.416410,0.362707,0.386603,0.469580,0.505383}; +// PROTONS +const Int_t kPBins_PbPb = 9; +const Double_t kPBinWidth_PbPb = 0.500; +const Double_t kPTagP_PbPb[kPBins_PbPb] = {0.017940,0.350681,0.535286,0.583264,0.562935,0.560524,0.545992,0.598060,0.351245}; +const Double_t kPTagPi_PbPb[kPBins_PbPb] = {0.195955,0.094949,0.039962,0.026039,0.007556,0.016986,0.030333,0.000000,0.000000}; +const Double_t kPTagNid_PbPb[kPBins_PbPb] = {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 kPiBins_pp = 10; +const Double_t kPiBinWidth_pp = 0.250; +const Double_t kPiTagPi_pp[kPiBins_pp] = {0.194528,0.447097,0.603364,0.646413,0.647125,0.669157,0.688139,0.682564,0.689910,0.665710}; +const Double_t kPiTagNid_pp[kPiBins_pp] = {0.805472,0.552903,0.396636,0.353587,0.352875,0.330843,0.311861,0.317436,0.310090,0.334290}; +// KAONS +const Int_t kKBins_pp = 10; +const Double_t kKBinWidth_pp = 0.250; +const Double_t kKTagK_pp[kKBins_pp] = {0.000000,0.173393,0.439690,0.519423,0.587025,0.605372,0.586021,0.650139,0.444444,0.299363}; +const Double_t kKTagPi_pp[kKBins_pp] = {0.000000,0.001495,0.000000,-0.000000,-0.000000,0.000000,0.032258,0.060572,0.101449,0.242038}; +const Double_t kKTagNid_pp[kKBins_pp] = {1.000000,0.825112,0.560310,0.480577,0.412975,0.394628,0.381720,0.289289,0.454106,0.458599}; +// PROTONS +const Int_t kPBins_pp = 9; +const Double_t kPBinWidth_pp = 0.500; +const Double_t kPTagP_pp[kPBins_pp] = {0.029404,0.438640,0.613710,0.665152,0.634961,0.657711,0.703704,0.685714,0.235294}; +const Double_t kPTagPi_pp[kPBins_pp] = {0.000000,0.000000,0.000000,-0.000000,0.000000,0.000000,-0.000000,0.014286,-0.000000}; +const Double_t kPTagNid_pp[kPBins_pp] = {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="TOFparam_PbPb"); + Double_t ChildrenRelAngle() const; + void ComputeWgts(); + void CorrectWgt4BR(Double_t factor); + Double_t CosPointing() const; + Double_t CosPointingXY() const; + void CosThetaStar(Double_t&,Double_t&) 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&,Double_t&,TString sample) const; + void GetPrimaryVtx(Double_t vtx[3]) const + { vtx[0]=fV1x; vtx[1]=fV1y; vtx[2]=fV1z; return; } + void GetSecondaryVtx(Double_t vtx[3]) const + { vtx[0]=fV2x; vtx[1]=fV2y; vtx[2]=fV2z; return; } + + void InvMass(Double_t&,Double_t&) 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&,Int_t&) const; + void SetPrimaryVtx(Double_t vtx[3]) + { fV1x=vtx[0]; fV1y=vtx[1]; fV1z=vtx[2]; return; } + void SetSignal() { fSignal = kTRUE; return; } + void SetTOFmasses(Double_t mass[2]) + { fTOFmass[0]=mass[0]; fTOFmass[1]=mass[1]; return; } + void SetPIDresponse(Double_t resp0[5],Double_t resp1[5]); + void SetPdgCodes(Int_t pdg[2]) {fPdg[0]=pdg[0];fPdg[1]=pdg[1];return;} + void SetMumPdgCodes(Int_t mum[2]) {fMum[0]=mum[0];fMum[1]=mum[1];return;} + + void DrawPIDinTOF(TString pidScheme="TOFparam_PbPb") const; + Double_t LinearInterpolation(Double_t p,Int_t nBins,Double_t Bin, + const Double_t *values) const; + // void SetPtWgts4pp(); + // + 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; // + Double_t fV1y; // position of the primary vertex of the event + Double_t fV1z; // + Double_t fV2x; // + Double_t fV2y; // position of the reconstructed secondary vertex + Double_t fV2z; // + Double_t fDCA; // DCA of the two tracks + + Double_t fPx[2]; // + 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; // + 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/ANALYSIS/AliD0toKpiAnalysis.cxx b/ANALYSIS/AliD0toKpiAnalysis.cxx new file mode 100644 index 00000000000..8f1fed8460b --- /dev/null +++ b/ANALYSIS/AliD0toKpiAnalysis.cxx @@ -0,0 +1,1092 @@ +/************************************************************************** + * 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) +// +// Origin: A. Dainese andrea.dainese@pd.infn.it +//---------------------------------------------------------------------------- +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "AliESD.h" +#include "AliStack.h" +#include "AliRunLoader.h" +#include "AliITStrackV2.h" +#include "AliITSVertexerTracks.h" +#include "AliITSVertex.h" +#include "AliV0vertexer.h" +#include "AliV0vertex.h" +#include "AliD0toKpi.h" +#include "AliD0toKpiAnalysis.h" + +typedef struct { + Int_t lab; + Int_t pdg; + Int_t mumlab; + Int_t mumpdg; + Float_t Vx,Vy,Vz; + Float_t Px,Py,Pz; +} REFTRACK; + +ClassImp(AliD0toKpiAnalysis) + +//---------------------------------------------------------------------------- +AliD0toKpiAnalysis::AliD0toKpiAnalysis() { + // Default constructor + + SetBz(); + SetPtCut(); + Setd0Cut(); + SetMassCut(); + SetD0Cuts(); + SetVertex1(); + SetPID(); + fVertexOnTheFly = kFALSE; + fSim = kFALSE; + fOnlySignal = kFALSE; + fDebug = 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 + + if(fBz<-9000.) { + printf("AliD0toKpiAnalysis::FindCandidates(): Set B!\n"); + return; + } + AliKalmanTrack::SetConvConst(100/0.299792458/fBz); + + TString trkName("AliITStracksV2.root"); + if(gSystem->AccessPathName(trkName.Data(),kFileExists)) { + printf("AliD0toKpiAnalysis::FindCandidates(): No tracks file found!\n"); + return; + } + + TString outName1=outName; + TString outName2("nTotEvents.dat"); + + Int_t ev; + Int_t nTotEv=0,nD0rec=0,nD0rec1ev=0; + Double_t dca; + Double_t v2[3],mom[6],d0[2]; + Double_t alphaP,alphaN,ptP,ptN,phiP,phiN; + Int_t iTrkP,iTrkN,trkEntries; + Int_t nTrksP=0,nTrksN=0; + Int_t trkNum[2]; + Int_t okD0=0,okD0bar=0; + Char_t trkTreeName[100],vtxName[100]; + AliITStrackV2 *postrack = 0; + AliITStrackV2 *negtrack = 0; + + // create the AliITSVertexerTracks object + // (it will be used only if fVertexOnTheFly=kTrue) + AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks; + vertexer1->SetMinTracks(2); + vertexer1->SetDebug(0); + Int_t skipped[2]; + Bool_t goodVtx1; + + + // define the cuts for vertexing + Double_t vtxcuts[]={50., // max. allowed chi2 + 0.0, // min. allowed negative daughter's impact param + 0.0, // min. allowed positive daughter's impact param + 1.0, // max. allowed DCA between the daughter tracks + -1.0, // min. allowed cosine of pointing angle + 0.0, // min. radius of the fiducial volume + 2.9};// max. radius of the fiducial volume + + // create the AliV0vertexer object + AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts); + + // 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 *trkFile = TFile::Open(trkName.Data()); + + // loop on events in file + for(ev=evFirst; ev<=evLast; ev++) { + printf(" --- Looking for D0->Kpi in event %d ---\n",ev); + + // retrieve primary vertex from file + sprintf(vtxName,"VertexTracks_%d",ev); + AliITSVertex *vertex1stored = (AliITSVertex*)trkFile->Get(vtxName); + vertex1stored->GetXYZ(fV1); + delete vertex1stored; + + // retrieve tracks from file + sprintf(trkTreeName,"TreeT_ITS_%d",ev); + TTree *trkTree=(TTree*)trkFile->Get(trkTreeName); + if(!trkTree) { + printf("AliD0toKpiAnalysis::FindCandidates():\n tracks tree not found for evet %d\n",ev); + continue; + } + trkEntries = (Int_t)trkTree->GetEntries(); + printf(" Number of tracks: %d\n",trkEntries); + + // count the total number of events + nTotEv++; + + // 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]; + SelectTracks(*trkTree,trksP,trkEntryP,nTrksP,trksN,trkEntryN,nTrksN); + + nD0rec1ev = 0; + + // loop on positive tracks + for(iTrkP=0; iTrkPPropagateToDCA(pnt,ppt); + + // define the AliV0vertex object + AliV0vertex *vertex2 = new AliV0vertex(*pnt,*ppt); + + // get position of the secondary vertex + vertex2->GetXYZ(v2[0],v2[1],v2[2]); + + delete vertex2; + + // momenta of the tracks at the vertex + ptP = 1./TMath::Abs(ppt->Get1Pt()); + alphaP = ppt->GetAlpha(); + phiP = alphaP+TMath::ASin(ppt->GetSnp()); + mom[0] = ptP*TMath::Cos(phiP); + mom[1] = ptP*TMath::Sin(phiP); + mom[2] = ptP*ppt->GetTgl(); + + ptN = 1./TMath::Abs(pnt->Get1Pt()); + alphaN = pnt->GetAlpha(); + phiN = alphaN+TMath::ASin(pnt->GetSnp()); + mom[3] = ptN*TMath::Cos(phiN); + mom[4] = ptN*TMath::Sin(phiN); + mom[5] = ptN*pnt->GetTgl(); + + goodVtx1 = kTRUE; + // no vertexing if DeltaMass > fMassCut + if(fVertexOnTheFly) { + goodVtx1 = kFALSE; + if(SelectInvMass(mom)) { + // primary vertex from *other* tracks in event + vertexer1->SetVtxStart(fV1[0],fV1[1]); + skipped[0] = trkEntryP[iTrkP]; + skipped[1] = trkEntryN[iTrkN]; + vertexer1->SetSkipTracks(2,skipped); + AliITSVertex *vertex1onfly = + (AliITSVertex*)vertexer1->VertexOnTheFly(*trkTree); + if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE; + vertex1onfly->GetXYZ(fV1); + //vertex1onfly->PrintStatus(); + delete vertex1onfly; + } + } + + // impact parameters of the tracks w.r.t. the primary vertex + d0[0] = 10000.*ppt->GetD(fV1[0],fV1[1]); + d0[1] = -10000.*pnt->GetD(fV1[0],fV1[1]); + + // create the object AliD0toKpi + AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0); + + // select D0s + if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) { + + // fill the tree + ioD0toKpi=&theD0; + treeD0->Fill(); + + nD0rec++; nD0rec1ev++; + ioD0toKpi=0; + } + + negtrack = 0; + } // loop on negative tracks + postrack = 0; + } // loop on positive tracks + + trksP.Delete(); + trksN.Delete(); + 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; + delete vertexer2; + + trkFile->Close(); + + // create a copy of this class to be written to output file + AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis"); + copy->PrintStatus(); + + // 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"); + MakeTracksRefFile(evFirst,evLast); + 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::FindCandidatesESD(Int_t evFirst,Int_t evLast, + const Char_t *outName) { + // Find D0 candidates and calculate parameters + + if(fBz<-9000.) { + printf("AliD0toKpiAnalysis::FindCandidatesESD(): Set B!\n"); + return; + } + AliKalmanTrack::SetConvConst(100/0.299792458/fBz); + + 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"); + + Double_t covV1[6]; + Int_t nTotEv=0,nD0rec=0,nD0rec1ev=0; + Double_t dca; + Double_t v2[3],mom[6],d0[2]; + Double_t alphaP,alphaN,ptP,ptN,phiP,phiN; + 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; + AliITStrackV2 *postrack = 0; + AliITStrackV2 *negtrack = 0; + + // create the AliITSVertexerTracks object + // (it will be used only if fVertexOnTheFly=kTrue) + AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks; + vertexer1->SetMinTracks(2); + vertexer1->SetDebug(0); + Int_t skipped[2]; + Bool_t goodVtx1; + + + // define the cuts for vertexing + Double_t vtxcuts[]={50., // max. allowed chi2 + 0.0, // min. allowed negative daughter's impact param + 0.0, // min. allowed positive daughter's impact param + 1.0, // max. allowed DCA between the daughter tracks + -1.0, // min. allowed cosine of pointing angle + 0.0, // min. radius of the fiducial volume + 2.9};// max. radius of the fiducial volume + + // create the AliV0vertexer object + AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts); + + // 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()); + + TKey *key=0; + TIter next(esdFile->GetListOfKeys()); + // loop on events in file + while ((key=(TKey*)next())!=0) { + AliESD *event=(AliESD*)key->ReadObj(); + Int_t ev = (Int_t)event->GetEventNumber(); + if(evevLast) { delete event; continue; } + printf("--- Finding D0 -> Kpi in event %d\n",ev); + + // retrieve primary vertex from file + //sprintf(vtxName,"Vertex_%d",ev); + //AliITSVertex *vertex1stored = (AliITSVertex*)trkFile->Get(vtxName); + //vertex1stored->GetXYZ(fV1); + //delete vertex1stored; + event->GetVertex(fV1,covV1); + + trkEntries = event->GetNumberOfTracks(); + printf(" Number of tracks: %d\n",trkEntries); + + // count the total number of events + nTotEv++; + + // 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(); + if(fVertexOnTheFly) { + SelectTracksESDvtx(*event,trkTree,trksP,trkEntryP,nTrksP, + trksN,trkEntryN,nTrksN); + } else { + SelectTracksESD(*event,trksP,trkEntryP,nTrksP, + trksN,trkEntryN,nTrksN); + } + + if(fDebug) printf(" pos. tracks: %d neg .tracks: %d\n",nTrksP,nTrksN); + + nD0rec1ev = 0; + + // loop on positive tracks + for(iTrkP=0; iTrkPPropagateToDCA(pnt,ppt); + + // define the AliV0vertex object + AliV0vertex *vertex2 = new AliV0vertex(*pnt,*ppt); + + // get position of the secondary vertex + vertex2->GetXYZ(v2[0],v2[1],v2[2]); + + delete vertex2; + + // momenta of the tracks at the vertex + ptP = 1./TMath::Abs(ppt->Get1Pt()); + alphaP = ppt->GetAlpha(); + phiP = alphaP+TMath::ASin(ppt->GetSnp()); + mom[0] = ptP*TMath::Cos(phiP); + mom[1] = ptP*TMath::Sin(phiP); + mom[2] = ptP*ppt->GetTgl(); + + ptN = 1./TMath::Abs(pnt->Get1Pt()); + alphaN = pnt->GetAlpha(); + phiN = alphaN+TMath::ASin(pnt->GetSnp()); + mom[3] = ptN*TMath::Cos(phiN); + mom[4] = ptN*TMath::Sin(phiN); + mom[5] = ptN*pnt->GetTgl(); + + goodVtx1 = kTRUE; + // no vertexing if DeltaMass > fMassCut + if(fVertexOnTheFly) { + goodVtx1 = kFALSE; + if(SelectInvMass(mom)) { + // primary vertex from *other* tracks in event + vertexer1->SetVtxStart(fV1[0],fV1[1]); + skipped[0] = trkEntryP[iTrkP]; + skipped[1] = trkEntryN[iTrkN]; + vertexer1->SetSkipTracks(2,skipped); + AliITSVertex *vertex1onfly = + (AliITSVertex*)vertexer1->VertexOnTheFly(*trkTree); + if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE; + vertex1onfly->GetXYZ(fV1); + //vertex1onfly->PrintStatus(); + delete vertex1onfly; + } + } + + // impact parameters of the tracks w.r.t. the primary vertex + d0[0] = 10000.*ppt->GetD(fV1[0],fV1[1]); + d0[1] = -10000.*pnt->GetD(fV1[0],fV1[1]); + + // 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 + + trksP.Delete(); + trksN.Delete(); + delete [] trkEntryP; + delete [] trkEntryN; + delete trkTree; + delete event; + + + 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; + delete vertexer2; + + 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"); + MakeTracksRefFileESD(); + 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(" fBz = %f T\n",fBz); + 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(TTree &trkTree, + 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 = (Int_t)trkTree.GetEntries(); + + // trasfer tracks from tree to arrays + for(Int_t i=0; iGet1Pt()>0.) { // negative track + trksN.AddLast(track); + trkEntryN[nTrksN] = i; + nTrksN++; + } else { // positive track + trksP.AddLast(track); + trkEntryP[nTrksP] = i; + nTrksP++; + } + + } + + return; +} +//----------------------------------------------------------------------------- +void AliD0toKpiAnalysis::SelectTracksESD(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; iGetStatus(); + + if(!(status&AliESDtrack::kITSrefit)) continue; + + AliITStrackV2 *itstrack = new AliITStrackV2(*esdtrack); + + // single track selection + if(!SingleTrkCuts(*itstrack)) + { delete itstrack; continue; } + + if(itstrack->Get1Pt()>0.) { // negative track + trksN.AddLast(itstrack); + trkEntryN[nTrksN] = i; + nTrksN++; + } else { // positive track + trksP.AddLast(itstrack); + trkEntryP[nTrksP] = i; + nTrksP++; + } + + } // loop on ESD tracks + + return; +} +//----------------------------------------------------------------------------- +void AliD0toKpiAnalysis::SelectTracksESDvtx(AliESD &event,TTree *trkTree, + 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(); + + AliITStrackV2 *itstrackfortree = 0; + trkTree->Branch("tracks","AliITStrackV2",&itstrackfortree,entr,0); + + + // transfer ITS tracks from ESD to arrays and to a tree + for(Int_t i=0; iGetStatus(); + + if(!(status&AliESDtrack::kITSrefit)) continue; + + AliITStrackV2 *itstrack = new AliITStrackV2(*esdtrack); + + // store track in the tree to be used for primary vertex finding + itstrackfortree = new AliITStrackV2(*esdtrack); + trkTree->Fill(); + itstrackfortree = 0; + + // single track selection + if(!SingleTrkCuts(*itstrack)) + { delete itstrack; continue; } + + if(itstrack->Get1Pt()>0.) { // negative track + trksN.AddLast(itstrack); + trkEntryN[nTrksN] = i; + nTrksN++; + } else { // positive track + trksP.AddLast(itstrack); + trkEntryP[nTrksP] = i; + nTrksP++; + } + + } // loop on esd tracks + + delete itstrackfortree; + + 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 AliITStrackV2& trk) const { + // Check if track passes some kinematical cuts + + if(TMath::Abs(1./trk.Get1Pt()) < fPtCut) + return kFALSE; + if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1])) < fd0Cut) + return kFALSE; + + return kTRUE; +} +//---------------------------------------------------------------------------- +void AliD0toKpiAnalysis::MakeTracksRefFile(Int_t evFirst,Int_t evLast) + const { + // Create a file with simulation info for the reconstructed tracks + + TFile *out = TFile::Open("ITStracksRefFile.root","recreate"); + TFile *trk = TFile::Open("AliITStracksV2.root"); + AliRunLoader *rl = AliRunLoader::Open("galice.root"); + + // load kinematics + rl->LoadKinematics(); + + Int_t label; + TParticle *part; + TParticle *mumpart; + REFTRACK reftrk; + + for(Int_t ev=evFirst; ev<=evLast; ev++){ + rl->GetEvent(ev); + AliStack *stack = rl->Stack(); + + trk->cd(); + + // Tree with tracks + char tname[100]; + sprintf(tname,"TreeT_ITS_%d",ev); + + TTree *tracktree=(TTree*)trk->Get(tname); + if(!tracktree) continue; + AliITStrackV2 *track = new AliITStrackV2; + tracktree->SetBranchAddress("tracks",&track); + Int_t nentr=(Int_t)tracktree->GetEntries(); + + // 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; iGetEvent(i); + label = TMath::Abs(track->GetLabel()); + + part = (TParticle*)stack->Particle(label); + reftrk.lab = label; + reftrk.pdg = part->GetPdgCode(); + reftrk.mumlab = part->GetFirstMother(); + if(part->GetFirstMother()>=0) { + mumpart = (TParticle*)stack->Particle(part->GetFirstMother()); + reftrk.mumpdg = mumpart->GetPdgCode(); + } else { + reftrk.mumpdg=-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 + + out->cd(); + reftree->Write(); + + delete track; + delete reftree; + delete tracktree; + delete stack; + } // loop on events + + trk->Close(); + out->Close(); + + return; +} +//---------------------------------------------------------------------------- +void AliD0toKpiAnalysis::MakeTracksRefFileESD() const { + // Create a file with simulation info for the reconstructed tracks + + TFile *outFile = TFile::Open("ITStracksRefFile.root","recreate"); + TFile *esdFile = TFile::Open("AliESDs.root"); + AliRunLoader *rl = AliRunLoader::Open("galice.root"); + + // load kinematics + rl->LoadKinematics(); + + Int_t label; + TParticle *part; + TParticle *mumpart; + REFTRACK reftrk; + + TKey *key=0; + TIter next(esdFile->GetListOfKeys()); + // loop on events in file + while ((key=(TKey*)next())!=0) { + AliESD *event=(AliESD*)key->ReadObj(); + Int_t ev = (Int_t)event->GetEventNumber(); + + rl->GetEvent(ev); + AliStack *stack = rl->Stack(); + + 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*)stack->Particle(label); + reftrk.lab = label; + reftrk.pdg = part->GetPdgCode(); + reftrk.mumlab = part->GetFirstMother(); + if(part->GetFirstMother()>=0) { + mumpart = (TParticle*)stack->Particle(part->GetFirstMother()); + reftrk.mumpdg = mumpart->GetPdgCode(); + } else { + reftrk.mumpdg=-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; + delete event; + delete stack; + } // 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("ITStracksRefFile.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]) theD0->SetSignal(); + + if(!fOnlySignal || theD0->IsSignal()) treeD0out->Fill(); + + } + + delete refTree0; + + refFile->Close(); + + return; +} + + + + + + diff --git a/ANALYSIS/AliD0toKpiAnalysis.h b/ANALYSIS/AliD0toKpiAnalysis.h new file mode 100644 index 00000000000..cf87106e30e --- /dev/null +++ b/ANALYSIS/AliD0toKpiAnalysis.h @@ -0,0 +1,109 @@ +#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@pd.infn.it +//------------------------------------------------------------------------- + +#include +#include +#include +#include "AliESD.h" +#include "AliITStrackV2.h" +#include "AliITSVertexerTracks.h" +#include "AliITSVertex.h" +#include "AliV0vertexer.h" +#include "AliV0vertex.h" +#include "AliD0toKpi.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 FindCandidatesESD(Int_t evFirst=0,Int_t evLast=0, + const Char_t *outName="AliD0toKpi.root"); + void PrintStatus() const; + void SetBz(Double_t bz=-9999.) { fBz=bz; return; } + void SetVertexOnTheFly() { fVertexOnTheFly=kTRUE; return; } + void SetSimulation() { fSim=kTRUE; return; } + void SetOnlySignal() { fOnlySignal=kTRUE; return; } + void SetPtCut(Double_t pt=0.) { fPtCut=pt; return; } + void Setd0Cut(Double_t d0=0.) { fd0Cut=d0; return; } + void SetMassCut(Double_t deltaM=1000.) { fMassCut=deltaM; return; } + 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(TString pid="TOFparam_PbPb") { fPID=pid; return; } + void SetDebug() { fDebug=kTRUE; return; } + // + private: + // + Double_t fBz; // value of the magnetic field + 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 + + Bool_t fDebug; // debug option + // + Double_t CalculateTOFmass(Double_t mom,Double_t length,Double_t time) const; + void MakeTracksRefFile(Int_t evFirst=0,Int_t evLast=0) const; + void MakeTracksRefFileESD() const; + Bool_t SelectInvMass(const Double_t p[6]) const; + void SelectTracks(TTree &trkTree, + TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP, + TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const; + void SelectTracksESD(AliESD &event, + TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP, + TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const; + void SelectTracksESDvtx(AliESD &event,TTree *trkTree, + 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; return; } + void SimulationInfo(TTree *treeD0in,TTree *treeD0out) const; + Bool_t SingleTrkCuts(const AliITStrackV2& trk) const; + // + ClassDef(AliD0toKpiAnalysis,1) // Reconstruction of D0 candidates class +}; + + +#endif + + + + + + + + diff --git a/ANALYSIS/AliD0toKpiPlots.C b/ANALYSIS/AliD0toKpiPlots.C new file mode 100644 index 00000000000..131ac179a5b --- /dev/null +++ b/ANALYSIS/AliD0toKpiPlots.C @@ -0,0 +1,232 @@ +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@pd.infn.it + //-------------------------------------------------------------------------- + + // 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/ANALYSIS/AliD0toKpiTest.C b/ANALYSIS/AliD0toKpiTest.C new file mode 100644 index 00000000000..7f7484f3313 --- /dev/null +++ b/ANALYSIS/AliD0toKpiTest.C @@ -0,0 +1,83 @@ +//-------------------------------------------------------------------------- +// Test macro for reconstruction and analysis of D0->Kpi +// +// Andrea Dainese, andrea.dainese@pd.infn.it +//-------------------------------------------------------------------------- + +void AliD0toKpiReco() { + + //============== R E C O N S T R U C T I O N ============================== + + // Look for field value in galice.root + Double_t field = 0.4; + if (gAlice) { + delete gAlice->GetRunLoader(); + delete gAlice; + gAlice=0; + } + if(!gSystem->AccessPathName("galice.root",kFileExists)) { + AliRunLoader *rl = AliRunLoader::Open("galice.root"); + rl->LoadgAlice(); + field=(Double_t)(gAlice->Field()->SolenoidField())/10.; + printf(" B = %3.1f T read from gAlice and set\n",field); + delete gAlice->GetRunLoader(); + delete gAlice; + gAlice=0; + } else { + printf(" File galice.root not found: default %3.1f T being used!\n",field); + } + + AliD0toKpiAnalysis *analysis = new AliD0toKpiAnalysis(); + //--- set magnetic field + analysis->SetBz(field); + // set simulation to take info on PDG codes from kine tree + //analysis->SetSimulation(); + //--- set this is you want only signal candidates in output file + //analysis->SetOnlySignal(); + //--- set this if you want to compute primary vertex D0 by D0 using + // other tracks in the event (for pp, broad interaction region); + // it is time-consuming procedure, so it can be done after a + // preselection on invariant mass + //analysis->SetVertexOnTheFly(); + //analysis->SetMassCut(0.1); // GeV + //--- set single-track preselections + analysis->SetPtCut(0.5); // GeV + analysis->Setd0Cut(50.); // micron + //--- set cuts on D0 candidates to be written to file + // (see AliD0toKpiAnalysis.h for a description and for the defaults) + analysis->SetD0Cuts(0.1,1000.,1.1,0.,0.,10000.,10000.,0.,.5); + //analysis->SetD0Cuts(); + + //--- check the settings + analysis->PrintStatus(); + + Int_t evFirst = 0; + Int_t evLast = 0; + //analysis->SetDebug(); + //analysis->FindCandidates(evFirst,evLast,"AliD0toKpi.root"); + analysis->FindCandidatesESD(evFirst,evLast,"AliD0toKpi.root"); + delete analysis; + + return; +} +//========================================================================== +void AliD0toKpiSele() { + + //======================== S E L E C T I O N ============================ + + AliD0toKpiAnalysis *analysis = new AliD0toKpiAnalysis(); + analysis->SetSimulation(); + analysis->SetOnlySignal(); + analysis->SetD0Cuts(0.1,1000.,1.1,0.,0.,10000.,10000.,0.,.5); + analysis->ApplySelection("AliD0toKpi.root","AliD0toKpi_sele.root"); + delete analysis; + + return; +} +//========================================================================== + + + + + + diff --git a/ANALYSIS/libANALYSIS.pkg b/ANALYSIS/libANALYSIS.pkg new file mode 100644 index 00000000000..7207543d3ec --- /dev/null +++ b/ANALYSIS/libANALYSIS.pkg @@ -0,0 +1,6 @@ +SRCS= AliD0toKpi.cxx AliD0toKpiAnalysis.cxx + + +HDRS= $(SRCS:.cxx=.h) + +DHDR:=ANALYSISLinkDef.h