--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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 <Riostream.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TCanvas.h>
+#include <TPaveLabel.h>
+
+#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<<fWgtAD0<<" "<<fWgtAD0bar<<endl;
+ cerr<<fWgtBD0<<" "<<fWgtBD0bar<<endl;
+ cerr<<fWgtCD0<<" "<<fWgtCD0bar<<endl;
+
+ if(fWgtAD0<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
+ if(fWgtAD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
+ if(fWgtBD0<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
+ if(fWgtBD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
+ if(fWgtCD0<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
+ if(fWgtCD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
+ */
+
+ return;
+}
+//----------------------------------------------------------------------------
+void AliD0toKpi::CorrectWgt4BR(Double_t factor) {
+ // correct weights of background from charm
+
+ fWgtAD0 *= factor;
+ fWgtAD0bar *= factor;
+ fWgtBD0 *= factor;
+ fWgtBD0bar *= factor;
+ fWgtCD0 *= factor;
+ fWgtCD0bar *= factor;
+ fWgtDD0 *= factor;
+ fWgtDD0bar *= factor;
+
+ return;
+}
+//----------------------------------------------------------------------------
+Double_t AliD0toKpi::CosPointing() const {
+ // cosine of pointing angle in space
+
+ TVector3 mom(Px(),Py(),Pz());
+ TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
+
+ Double_t pta = mom.Angle(flight);
+
+ return TMath::Cos(pta);
+}
+//----------------------------------------------------------------------------
+Double_t AliD0toKpi::CosPointingXY() const {
+ // cosine of pointing angle in transverse plane
+
+ TVector3 momXY(Px(),Py(),0.);
+ TVector3 flightXY(fV2x-fV1x,fV2y-fV1y,0.);
+
+ Double_t ptaXY = momXY.Angle(flightXY);
+
+ return TMath::Cos(ptaXY);
+}
+//----------------------------------------------------------------------------
+void AliD0toKpi::CosThetaStar(Double_t &ctsD0,Double_t &ctsD0bar) const {
+ // cosine of decay angle in the D0 rest frame
+
+ Double_t pStar = TMath::Sqrt(TMath::Power(kMD0*kMD0-kMK*kMK-kMPi*kMPi,2.)-4.*kMK*kMK*kMPi*kMPi)/(2.*kMD0);
+
+ Double_t beta = P()/Energy();
+ Double_t gamma = Energy()/kMD0;
+
+ ctsD0 = (Ql(1)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
+ // if(ctsD0 > 1.) { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0<<"!\n"; }
+ // if(ctsD0 < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0<<"!\n"; }
+
+ ctsD0bar = (Ql(0)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
+ // if(ctsD0bar > 1.) { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0bar<<"!\n"; }
+ // if(ctsD0bar < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0bar<<"!\n";}
+
+ return;
+}
+//----------------------------------------------------------------------------
+Double_t AliD0toKpi::Eta() const {
+ // pseudorapidity of the D0
+
+ Double_t theta = TMath::Pi()/2.-TMath::ATan2(Pz(),Pt());
+ Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
+ return eta;
+}
+//----------------------------------------------------------------------------
+Double_t AliD0toKpi::EtaChild(Int_t child) const {
+ // pseudorapidity of the decay tracks
+
+ Double_t theta = TMath::Pi()/2.-TMath::ATan2(fPz[child],PtChild(child));
+ Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
+ return eta;
+}
+//----------------------------------------------------------------------------
+void AliD0toKpi::GetWgts(Double_t &WgtD0,Double_t &WgtD0bar,TString sample)
+ const {
+ // returns the weights for pid
+
+ const char *sampleA = strstr(sample.Data(),"A");
+ const char *sampleB = strstr(sample.Data(),"B");
+ const char *sampleC = strstr(sample.Data(),"C");
+ const char *sampleD = strstr(sample.Data(),"D");
+ const char *sampleABCD = strstr(sample.Data(),"ABCD");
+ const char *sampleABC = strstr(sample.Data(),"ABC");
+ const char *sampleBC = strstr(sample.Data(),"BC");
+
+ if(sampleA) { WgtD0 = fWgtAD0; WgtD0bar = fWgtAD0bar; }
+ if(sampleB) { WgtD0 = fWgtBD0; WgtD0bar = fWgtBD0bar; }
+ if(sampleC) { WgtD0 = fWgtCD0; WgtD0bar = fWgtCD0bar; }
+ if(sampleD) { WgtD0 = fWgtDD0; WgtD0bar = fWgtDD0bar; }
+ if(sampleABCD) {
+ WgtD0 = fWgtAD0+fWgtBD0+fWgtCD0+fWgtDD0;
+ WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar+fWgtDD0bar;
+ }
+ if(sampleABC) {
+ WgtD0 = fWgtAD0+fWgtBD0+fWgtCD0;
+ WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar;
+ }
+ if(sampleBC) {
+ WgtD0 = fWgtBD0+fWgtCD0;
+ WgtD0bar = fWgtBD0bar+fWgtCD0bar;
+ }
+
+
+ if(fSignal) {
+ if(fMum[0]==421) WgtD0bar = 0.;
+ if(fMum[0]==-421) WgtD0 = 0.;
+ }
+
+ return;
+}
+//----------------------------------------------------------------------------
+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;
+}
+//-----------------------------------------------------------------------------
+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; i<nBins; i++) {
+ if(p<(i+0.5)*Bin) {
+ slope = (values[i]-values[i-1])/Bin;
+ value = values[i-1]+slope*(p-Bin*(i-0.5));
+ break;
+ }
+ }
+ }
+
+ if(value<0.) value=0.;
+ if(value>1.) value=1.;
+
+ return value;
+}
+//----------------------------------------------------------------------------
+
+
+
+
+
+
+/*
+//____________________________________________________________________________
+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;
+}
+//____________________________________________________________________________
+*/
+
+
+
+
+
+
+
+
+
--- /dev/null
+#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 <TMath.h>
+#include <TCanvas.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TPaveLabel.h>
+#include <TVector3.h>
+
+//----------------------------------------------------------------------------
+// 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
+
+
+
+
+
+
+
+
--- /dev/null
+/**************************************************************************
+ * 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 <Riostream.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TKey.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TString.h>
+#include <TSystem.h>
+#include <TParticle.h>
+#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; i<entries; i++) {
+ // get event from tree
+ treeD0in->GetEvent(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; iTrkP<nTrksP; iTrkP++) {
+ if(iTrkP%10==0) printf(" Processing positive track number %d of %d\n",iTrkP,nTrksP);
+
+ // get track from track array
+ postrack = (AliITStrackV2*)trksP.At(iTrkP);
+ trkNum[0] = trkEntryP[iTrkP];
+
+ // loop on negative tracks
+ for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
+ // get track from tracks array
+ negtrack = (AliITStrackV2*)trksN.At(iTrkN);
+ trkNum[1] = trkEntryN[iTrkN];
+
+ AliITStrackV2 nt(*negtrack), pt(*postrack), *pnt=&nt, *ppt=&pt;
+
+ //
+ // ----------- DCA MINIMIZATION ------------------
+ //
+ // find the DCA and propagate the tracks to the DCA
+ dca = vertexer2->PropagateToDCA(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(ev<evFirst || ev>evLast) { 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; iTrkP<nTrksP; iTrkP++) {
+ if((iTrkP%10==0) || fDebug) printf(" Processing positive track number %d of %d\n",iTrkP,nTrksP);
+
+ // get track from track array
+ postrack = (AliITStrackV2*)trksP.UncheckedAt(iTrkP);
+ trkNum[0] = trkEntryP[iTrkP];
+
+ // loop on negative tracks
+ for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
+
+ // get track from tracks array
+ negtrack = (AliITStrackV2*)trksN.UncheckedAt(iTrkN);
+ trkNum[1] = trkEntryN[iTrkN];
+
+ AliITStrackV2 nt(*negtrack), pt(*postrack), *pnt=&nt, *ppt=&pt;
+
+ //
+ // ----------- DCA MINIMIZATION ------------------
+ //
+ // find the DCA and propagate the tracks to the DCA
+ dca = vertexer2->PropagateToDCA(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; i<entr; i++) {
+
+ AliITStrackV2 *track = new AliITStrackV2;
+ trkTree.SetBranchAddress("tracks",&track);
+
+ trkTree.GetEvent(i);
+
+ // single track selection
+ if(!SingleTrkCuts(*track)) { delete track; continue; }
+
+ if(track->Get1Pt()>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; i<entr; i++) {
+ //if(fDebug) printf(" SelectTracksESD: %d/%d\n",i,entr);
+
+ AliESDtrack *esdtrack = (AliESDtrack*)event.GetTrack(i);
+ UInt_t status = esdtrack->GetStatus();
+
+ 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; i<entr; i++) {
+
+ AliESDtrack *esdtrack = (AliESDtrack*)event.GetTrack(i);
+ UInt_t status = esdtrack->GetStatus();
+
+ 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; i<nentr; i++) {
+ tracktree->GetEvent(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; i<nentr; i++) {
+ AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(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; i<entries; i++) {
+ if(i%100==0) printf(" done %d candidates of %d\n",i,entries);
+
+ treeD0in->GetEvent(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;
+}
+
+
+
+
+
+
--- /dev/null
+#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 <Riostream.h>
+#include <TString.h>
+#include <TNamed.h>
+#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
+
+
+
+
+
+
+
+
--- /dev/null
+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; i<entries; i++) {
+ if(i%10000==0) printf(" candidate %d of %d\n",i,entries);
+
+ // get event from tree
+ treeD0->GetEvent(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;
+}
+
+
+
+
--- /dev/null
+//--------------------------------------------------------------------------
+// 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;
+}
+//==========================================================================
+
+
+
+
+
+
--- /dev/null
+SRCS= AliD0toKpi.cxx AliD0toKpiAnalysis.cxx
+
+
+HDRS= $(SRCS:.cxx=.h)
+
+DHDR:=ANALYSISLinkDef.h
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff