+++ /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
-// for pp and PbPb interactions
-// Note: the two decay tracks are labelled: 0 (positive track)
-// 1 (negative track)
-// Origin: A. Dainese andrea.dainese@lnl.infn.it
-//----------------------------------------------------------------------------
-
-#include <TH1.h>
-#include <TH2.h>
-#include <TCanvas.h>
-#include <TPaveLabel.h>
-#include <TVector3.h>
-
-#include "AliD0toKpi.h"
-
-ClassImp(AliD0toKpi)
-
-//----------------------------------------------------------------------------
-AliD0toKpi::AliD0toKpi() {
- // Default constructor
-
- fSignal = kFALSE;
-
- fEvent = 0;
-
- fTrkNum[0] = 0;
- fTrkNum[1] = 0;
-
- fV1x = 0.;
- fV1y = 0.;
- fV1z = 0.;
- fV2x = 0.;
- fV2y = 0.;
- fV2z = 0.;
- fDCA = 0.;
-
- fPx[0] = 0.;
- fPy[0] = 0.;
- fPz[0] = 0.;
- fPx[1] = 0.;
- fPy[1] = 0.;
- fPz[1] = 0.;
-
- fd0[0] = 0.;
- fd0[1] = 0.;
-
- fPdg[0] = 0;
- fPdg[1] = 0;
- fMum[0] = 0;
- fMum[1] = 0;
-
- fTagPi[0] = 0.;
- fTagPi[1] = 0.;
- fTagKa[0] = 0.;
- fTagKa[1] = 0.;
- fTagNid[0] = 0.;
- fTagNid[1] = 0.;
-
- fWgtAD0=fWgtAD0bar=fWgtBD0=fWgtBD0bar=fWgtCD0=fWgtCD0bar=fWgtDD0=fWgtDD0bar=0;
-
-}
-//----------------------------------------------------------------------------
-AliD0toKpi::AliD0toKpi(Int_t ev,Int_t trkNum[2],
- Double_t v1[3],Double_t v2[3],
- Double_t dca,
- Double_t mom[6],Double_t d0[2]) {
- // Constructor
-
- fSignal = kFALSE;
-
- fEvent = ev;
- fTrkNum[0] = trkNum[0];
- fTrkNum[1] = trkNum[1];
-
- fV1x = v1[0];
- fV1y = v1[1];
- fV1z = v1[2];
- fV2x = v2[0];
- fV2y = v2[1];
- fV2z = v2[2];
- fDCA = dca;
-
- fPx[0] = mom[0];
- fPy[0] = mom[1];
- fPz[0] = mom[2];
- fPx[1] = mom[3];
- fPy[1] = mom[4];
- fPz[1] = mom[5];
-
- fd0[0] = d0[0];
- fd0[1] = d0[1];
-
- fPdg[0] = 0;
- fPdg[1] = 0;
- fMum[0] = 0;
- fMum[1] = 0;
-
- fTagPi[0] = 0.;
- fTagPi[1] = 0.;
- fTagKa[0] = 0.;
- fTagKa[1] = 0.;
- fTagNid[0] = 0.;
- fTagNid[1] = 0.;
-
- fWgtAD0=fWgtAD0bar=fWgtBD0=fWgtBD0bar=fWgtCD0=fWgtCD0bar=fWgtDD0=fWgtDD0bar=0;
-}
-//----------------------------------------------------------------------------
-AliD0toKpi::~AliD0toKpi() {}
-//____________________________________________________________________________
-AliD0toKpi::AliD0toKpi( const AliD0toKpi& d0toKpi):TObject(d0toKpi) {
- // dummy copy constructor
-}
-//----------------------------------------------------------------------------
-void AliD0toKpi::ApplyPID(TString pidScheme) {
- // Applies particle identification
-
- if((!pidScheme.CompareTo("TOFparamPbPb") || !pidScheme.CompareTo("TOFparamPP")) && fPdg[0]==0) {
- printf("AliD0toKpi::ApplyPID :\n Warning: TOF parameterized PID can be used only for simulation!\n");
- return;
- }
-
- if(!pidScheme.CompareTo("TOFparamPbPb")) {
- // tagging of the positive track
- if(TMath::Abs(fPdg[0])==211 || TMath::Abs(fPdg[0])==13
- || TMath::Abs(fPdg[0])==11) { // pion,muon,electron
- fTagPi[0] = LinearInterpolation(PChild(0),kPiBinsPbPb,kPiBinWidthPbPb,kPiTagPiPbPb);
- fTagNid[0] = 1.-fTagPi[0];
- fTagKa[0] = 0.;
- fTagPr[0] = 0.;
- }
- if(TMath::Abs(fPdg[0])==321) { // kaon
- fTagKa[0] = LinearInterpolation(PChild(0),kKBinsPbPb,kKBinWidthPbPb,kKTagKPbPb);
- fTagNid[0] = LinearInterpolation(PChild(0),kKBinsPbPb,kKBinWidthPbPb,kKTagNidPbPb);
- if((fTagNid[0]+fTagKa[0])>1.) fTagNid[0] = 1.-fTagKa[0];
- fTagPi[0] = 1.-fTagNid[0]-fTagKa[0];
- fTagPr[0] = 0.;
- }
- if(TMath::Abs(fPdg[0])==2212) { // proton
- fTagPr[0] = LinearInterpolation(PChild(0),kPBinsPbPb,kPBinWidthPbPb,kPTagPPbPb);
- fTagNid[0] = LinearInterpolation(PChild(0),kPBinsPbPb,kPBinWidthPbPb,kPTagNidPbPb);
- if((fTagNid[0]+fTagPr[0])>1.) fTagNid[0] = 1.-fTagPr[0];
- fTagPi[0] = 1.-fTagNid[0]-fTagPr[0];
- fTagKa[0] = 0.;
- }
- // tagging of the negative track
- if(TMath::Abs(fPdg[1])==211 || TMath::Abs(fPdg[1])==13
- || TMath::Abs(fPdg[1])==11) { // pion,muon,electron
- fTagPi[1] = LinearInterpolation(PChild(1),kPiBinsPbPb,kPiBinWidthPbPb,kPiTagPiPbPb);
- fTagNid[1] = 1.-fTagPi[1];
- fTagKa[1] = 0.;
- fTagPr[1] = 0.;
- }
- if(TMath::Abs(fPdg[1])==321) { // kaon
- fTagKa[1] = LinearInterpolation(PChild(1),kKBinsPbPb,kKBinWidthPbPb,kKTagKPbPb);
- fTagNid[1] = LinearInterpolation(PChild(1),kKBinsPbPb,kKBinWidthPbPb,kKTagNidPbPb);
- if((fTagNid[1]+fTagKa[1])>1.) fTagNid[1] = 1.-fTagKa[1];
- fTagPi[1] = 1.-fTagNid[1]-fTagKa[1];
- fTagPr[1] = 0.;
- }
- if(TMath::Abs(fPdg[1])==2212) { // proton
- fTagPr[1] = LinearInterpolation(PChild(1),kPBinsPbPb,kPBinWidthPbPb,kPTagPPbPb);
- fTagNid[1] = LinearInterpolation(PChild(1),kPBinsPbPb,kPBinWidthPbPb,kPTagNidPbPb);
- if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1];
- fTagPi[1] = 1.-fTagNid[1]-fTagPr[1];
- fTagKa[1] = 0.;
- }
- }
-
-
- if(!pidScheme.CompareTo("TOFparamPP")) {
- // tagging of the positive track
- if(TMath::Abs(fPdg[0])==211 || TMath::Abs(fPdg[0])==13
- || TMath::Abs(fPdg[0])==11) { // pion,muon,electron
- fTagPi[0] = LinearInterpolation(PChild(0),kPiBinsPP,kPiBinWidthPP,kPiTagPiPP);
- fTagNid[0] = 1.-fTagPi[0];
- fTagKa[0] = 0.;
- fTagPr[0] = 0.;
- }
- if(TMath::Abs(fPdg[0])==321) { // kaon
- fTagKa[0] = LinearInterpolation(PChild(0),kKBinsPP,kKBinWidthPP,kKTagKPP);
- fTagNid[0] = LinearInterpolation(PChild(0),kKBinsPP,kKBinWidthPP,kKTagNidPP);
- if((fTagNid[0]+fTagKa[0])>1.) fTagNid[0] = 1.-fTagKa[0];
- fTagPi[0] = 1.-fTagNid[0]-fTagKa[0];
- fTagPr[0] = 0.;
- }
- if(TMath::Abs(fPdg[0])==2212) { // proton
- fTagPr[0] = LinearInterpolation(PChild(0),kPBinsPP,kPBinWidthPP,kPTagPPP);
- fTagNid[0] = LinearInterpolation(PChild(0),kPBinsPP,kPBinWidthPP,kPTagNidPP);
- if((fTagNid[0]+fTagPr[0])>1.) fTagNid[0] = 1.-fTagPr[0];
- fTagPi[0] = 1.-fTagNid[0]-fTagPr[0];
- fTagKa[0] = 0.;
- }
- // tagging of the negative track
- if(TMath::Abs(fPdg[1])==211 || TMath::Abs(fPdg[1])==13
- || TMath::Abs(fPdg[1])==11) { // pion,muon,electron
- fTagPi[1] = LinearInterpolation(PChild(1),kPiBinsPP,kPiBinWidthPP,kPiTagPiPP);
- fTagNid[1] = 1.-fTagPi[1];
- fTagKa[1] = 0.;
- fTagPr[1] = 0.;
- }
- if(TMath::Abs(fPdg[1])==321) { // kaon
- fTagKa[1] = LinearInterpolation(PChild(1),kKBinsPP,kKBinWidthPP,kKTagKPP);
- fTagNid[1] = LinearInterpolation(PChild(1),kKBinsPP,kKBinWidthPP,kKTagNidPP);
- if((fTagNid[1]+fTagKa[1])>1.) fTagNid[1] = 1.-fTagKa[1];
- fTagPi[1] = 1.-fTagNid[1]-fTagKa[1];
- fTagPr[1] = 0.;
- }
- if(TMath::Abs(fPdg[1])==2212) { // proton
- fTagPr[1] = LinearInterpolation(PChild(1),kPBinsPP,kPBinWidthPP,kPTagPPP);
- fTagNid[1] = LinearInterpolation(PChild(1),kPBinsPP,kPBinWidthPP,kPTagNidPP);
- if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1];
- fTagPi[1] = 1.-fTagNid[1]-fTagPr[1];
- fTagKa[1] = 0.;
- }
- }
-
- return;
-}
-//----------------------------------------------------------------------------
-Double_t AliD0toKpi::ChildrenRelAngle() const {
- // relative angle between K and pi
-
- TVector3 mom0(fPx[0],fPy[0],fPz[0]);
- TVector3 mom1(fPx[1],fPy[1],fPz[1]);
-
- Double_t angle = mom0.Angle(mom1);
-
- return angle;
-}
-//----------------------------------------------------------------------------
-void AliD0toKpi::ComputeWgts() {
- // calculate the weights for PID
-
-
- // assignement of the weights from PID
- fWgtAD0 = fTagKa[1]*(fTagPi[0]+fTagNid[0]);
- fWgtAD0bar = fTagKa[0]*(fTagPi[1]+fTagNid[1]);
- fWgtBD0 = fTagPi[0]*fTagNid[1];
- fWgtBD0bar = fTagPi[1]*fTagNid[0];
- fWgtCD0 = fTagNid[0]*fTagNid[1];
- fWgtCD0bar = fTagNid[0]*fTagNid[1];
- fWgtDD0 = 1.-fWgtAD0-fWgtBD0-fWgtCD0;
- fWgtDD0bar = 1.-fWgtAD0bar-fWgtBD0bar-fWgtCD0bar;
-
- /*
- cerr<<fWgtAD0<<" "<<fWgtAD0bar<<endl;
- cerr<<fWgtBD0<<" "<<fWgtBD0bar<<endl;
- cerr<<fWgtCD0<<" "<<fWgtCD0bar<<endl;
- cerr<<fWgtDD0<<" "<<fWgtDD0bar<<endl;
- */
- /*
- if(fWgtAD0<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
- if(fWgtAD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
- if(fWgtBD0<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
- if(fWgtBD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
- if(fWgtCD0<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
- if(fWgtCD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts() Negative weight!!!\n";
- */
-
- return;
-}
-//----------------------------------------------------------------------------
-void AliD0toKpi::CorrectWgt4BR(Double_t factor) {
- // correct weights of background from charm
-
- fWgtAD0 *= factor;
- fWgtAD0bar *= factor;
- fWgtBD0 *= factor;
- fWgtBD0bar *= factor;
- fWgtCD0 *= factor;
- fWgtCD0bar *= factor;
- fWgtDD0 *= factor;
- fWgtDD0bar *= factor;
-
- return;
-}
-//----------------------------------------------------------------------------
-Double_t AliD0toKpi::CosPointing() const {
- // cosine of pointing angle in space
-
- TVector3 mom(Px(),Py(),Pz());
- TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
-
- Double_t pta = mom.Angle(flight);
-
- return TMath::Cos(pta);
-}
-//----------------------------------------------------------------------------
-Double_t AliD0toKpi::CosPointingXY() const {
- // cosine of pointing angle in transverse plane
-
- TVector3 momXY(Px(),Py(),0.);
- TVector3 flightXY(fV2x-fV1x,fV2y-fV1y,0.);
-
- Double_t ptaXY = momXY.Angle(flightXY);
-
- return TMath::Cos(ptaXY);
-}
-//----------------------------------------------------------------------------
-void AliD0toKpi::CosThetaStar(Double_t &ctsD0,Double_t &ctsD0bar) const {
- // cosine of decay angle in the D0 rest frame
-
- Double_t pStar = TMath::Sqrt(TMath::Power(kMD0*kMD0-kMK*kMK-kMPi*kMPi,2.)-4.*kMK*kMK*kMPi*kMPi)/(2.*kMD0);
-
- Double_t beta = P()/Energy();
- Double_t gamma = Energy()/kMD0;
-
- ctsD0 = (Ql(1)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
- // if(ctsD0 > 1.) { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0<<"!\n"; }
- // if(ctsD0 < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0<<"!\n"; }
-
- ctsD0bar = (Ql(0)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
- // if(ctsD0bar > 1.) { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0bar<<"!\n"; }
- // if(ctsD0bar < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0bar<<"!\n";}
-
- return;
-}
-//----------------------------------------------------------------------------
-Double_t AliD0toKpi::Eta() const {
- // pseudorapidity of the D0
-
- Double_t theta = TMath::Pi()/2.-TMath::ATan2(Pz(),Pt());
- Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
- return eta;
-}
-//----------------------------------------------------------------------------
-Double_t AliD0toKpi::EtaChild(Int_t child) const {
- // pseudorapidity of the decay tracks
-
- Double_t theta = TMath::Pi()/2.-TMath::ATan2(fPz[child],PtChild(child));
- Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
- return eta;
-}
-//----------------------------------------------------------------------------
-void AliD0toKpi::GetWgts(Double_t &WgtD0,Double_t &WgtD0bar,TString sample)
- const {
- // returns the weights for pid
-
- if(!sample.CompareTo("A")) { WgtD0 = fWgtAD0; WgtD0bar = fWgtAD0bar; }
- if(!sample.CompareTo("B")) { WgtD0 = fWgtBD0; WgtD0bar = fWgtBD0bar; }
- if(!sample.CompareTo("C")) { WgtD0 = fWgtCD0; WgtD0bar = fWgtCD0bar; }
- if(!sample.CompareTo("D")) { WgtD0 = fWgtDD0; WgtD0bar = fWgtDD0bar; }
- if(!sample.CompareTo("ABCD")) {
- WgtD0 = fWgtAD0+fWgtBD0+fWgtCD0+fWgtDD0;
- WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar+fWgtDD0bar;
- }
- if(!sample.CompareTo("ABC")) {
- WgtD0 = fWgtAD0+fWgtBD0+fWgtCD0;
- WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar;
- }
- if(!sample.CompareTo("BC")) {
- WgtD0 = fWgtBD0+fWgtCD0;
- WgtD0bar = fWgtBD0bar+fWgtCD0bar;
- }
-
- return;
-}
-//----------------------------------------------------------------------------
-Double_t AliD0toKpi::ImpPar() const {
- // D0 impact parameter in the bending plane
-
- Double_t k = -(fV2x-fV1x)*Px()-(fV2y-fV1y)*Py();
- k /= Pt()*Pt();
- Double_t dx = fV2x-fV1x+k*Px();
- Double_t dy = fV2y-fV1y+k*Py();
- Double_t absDD = TMath::Sqrt(dx*dx+dy*dy);
- TVector3 mom(Px(),Py(),Pz());
- TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
- TVector3 cross = mom.Cross(flight);
- return (cross.Z()>0. ? absDD : -absDD);
-}
-//----------------------------------------------------------------------------
-void AliD0toKpi::InvMass(Double_t &mD0,Double_t &mD0bar) const {
- // invariant mass as D0 and as D0bar
-
- Double_t energy[2];
-
- // D0 -> K- Pi+
- energy[1] = TMath::Sqrt(kMK*kMK+PChild(1)*PChild(1));
- energy[0] = TMath::Sqrt(kMPi*kMPi+PChild(0)*PChild(0));
-
- mD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
-
-
- // D0bar -> K+ Pi-
- energy[0] = TMath::Sqrt(kMK*kMK+PChild(0)*PChild(0));
- energy[1] = TMath::Sqrt(kMPi*kMPi+PChild(1)*PChild(1));
-
- mD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
-
- return;
-
-}
-//----------------------------------------------------------------------------
-Double_t AliD0toKpi::Ql(Int_t child) const {
- // longitudinal momentum of decay tracks w.r.t. to D0 momentum
-
- Double_t qL;
- TVector3 mom(fPx[child],fPy[child],fPz[child]);
- TVector3 momD(Px(),Py(),Pz());
-
- qL = mom.Dot(momD)/momD.Mag();
-
- return qL ;
-}
-//----------------------------------------------------------------------------
-Double_t AliD0toKpi::Qt() const {
- // transverse momentum of decay tracks w.r.t. to D0 momentum
-
- TVector3 mom0(fPx[0],fPy[0],fPz[0]);
- TVector3 momD(Px(),Py(),Pz());
-
- return mom0.Perp(momD);
-}
-//----------------------------------------------------------------------------
-Bool_t AliD0toKpi::Select(const Double_t* cuts,Int_t& okD0,Int_t& okD0bar)
- const {
-//
-// This function compares the D0 with a set of cuts:
-//
-// cuts[0] = inv. mass half width [GeV]
-// cuts[1] = dca [micron]
-// cuts[2] = cosThetaStar
-// cuts[3] = pTK [GeV/c]
-// cuts[4] = pTPi [GeV/c]
-// cuts[5] = d0K [micron] upper limit!
-// cuts[6] = d0Pi [micron] upper limit!
-// cuts[7] = d0d0 [micron^2]
-// cuts[8] = cosThetaPoint
-//
-// If the D0/D0bar doesn't pass the cuts it sets the weights to 0
-// If neither D0 nor D0bar pass the cuts return kFALSE
-//
- Double_t mD0,mD0bar,ctsD0,ctsD0bar;
- okD0=1; okD0bar=1;
-
- if(PtChild(1) < cuts[3] || PtChild(0) < cuts[4]) okD0 = 0;
- if(PtChild(0) < cuts[3] || PtChild(1) < cuts[4]) okD0bar = 0;
- if(!okD0 && !okD0bar) return kFALSE;
-
- if(TMath::Abs(Getd0Child(1)) > cuts[5] ||
- TMath::Abs(Getd0Child(0)) > cuts[6]) okD0 = 0;
- if(TMath::Abs(Getd0Child(0)) > cuts[6] ||
- TMath::Abs(Getd0Child(1)) > cuts[5]) okD0bar = 0;
- if(!okD0 && !okD0bar) return kFALSE;
-
- if(GetDCA() > cuts[1]) { okD0 = okD0bar = 0; return kFALSE; }
-
- InvMass(mD0,mD0bar);
- if(TMath::Abs(mD0-kMD0) > cuts[0]) okD0 = 0;
- if(TMath::Abs(mD0bar-kMD0) > cuts[0]) okD0bar = 0;
- if(!okD0 && !okD0bar) return kFALSE;
-
- CosThetaStar(ctsD0,ctsD0bar);
- if(TMath::Abs(ctsD0) > cuts[2]) okD0 = 0;
- if(TMath::Abs(ctsD0bar) > cuts[2]) okD0bar = 0;
- if(!okD0 && !okD0bar) return kFALSE;
-
- if(ProdImpParams() > cuts[7]) { okD0 = okD0bar = 0; return kFALSE; }
-
- if(CosPointing() < cuts[8]) { okD0 = okD0bar = 0; return kFALSE; }
-
- return kTRUE;
-}
-//-----------------------------------------------------------------------------
-void AliD0toKpi::SetPIDresponse(Double_t resp0[5],Double_t resp1[5]) {
- // Set combined PID detector response probabilities
-
- fPIDrespEl[0] = resp0[0];
- fPIDrespEl[1] = resp1[0];
- fPIDrespMu[0] = resp0[1];
- fPIDrespMu[1] = resp1[1];
- fPIDrespPi[0] = resp0[2];
- fPIDrespPi[1] = resp1[2];
- fPIDrespKa[0] = resp0[3];
- fPIDrespKa[1] = resp1[3];
- fPIDrespPr[0] = resp0[4];
- fPIDrespPr[1] = resp1[4];
-
- return;
-}
-//----------------------------------------------------------------------------
-Double_t AliD0toKpi::LinearInterpolation(Double_t p,Int_t nBins,Double_t Bin,
- const Double_t *values) const {
- // a linear interpolation method
-
- Double_t value=0;
- Double_t slope;
-
- if(p<0.5*Bin) {
- value = values[0];
- } else if(p>=(nBins-0.5)*Bin) {
- slope = (2*values[nBins-1]-values[nBins-2]-values[nBins-3])/Bin/2;
- value = values[nBins-2]+slope*(p-Bin*(nBins-1.5));
- } else {
- for(Int_t i=0; i<nBins; i++) {
- if(p<(i+0.5)*Bin) {
- slope = (values[i]-values[i-1])/Bin;
- value = values[i-1]+slope*(p-Bin*(i-0.5));
- break;
- }
- }
- }
-
- if(value<0.) value=0.;
- if(value>1.) value=1.;
-
- return value;
-}
-//----------------------------------------------------------------------------
-
-
-
-
-
-
-
-
-
+++ /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
-//
-// Note: the two decay tracks are labelled: 0 (positive track)
-// 1 (negative track)
-//
-// Origin: A. Dainese andrea.dainese@lnl.infn.it
-//-------------------------------------------------------------------------
-
-#include <TObject.h>
-#include <TMath.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 kPiBinsPbPb = 10;
-const Double_t kPiBinWidthPbPb = 0.250;
-const Double_t kPiTagPiPbPb[kPiBinsPbPb] = {0.211421,0.652184,0.624421,0.614727,0.610777,0.628015,0.631520,0.630324,0.637551,0.575235};
-const Double_t kPiTagNidPbPb[kPiBinsPbPb] = {0.788579,0.347816,0.375579,0.385273,0.389223,0.371985,0.368480,0.369676,0.362449,0.424765};
-// KAONS
-const Int_t kKBinsPbPb = 10;
-const Double_t kKBinWidthPbPb = 0.250;
-const Double_t kKTagKPbPb[kKBinsPbPb] = {0.000000,0.101255,0.397662,0.467586,0.517008,0.555023,0.584185,0.519029,0.464117,0.247308};
-const Double_t kKTagPiPbPb[kKBinsPbPb] = {0.102049,0.289930,0.101930,0.057771,0.040286,0.028567,0.053108,0.094369,0.066302,0.247308};
-const Double_t kKTagNidPbPb[kKBinsPbPb] = {0.897951,0.608815,0.500408,0.474643,0.442705,0.416410,0.362707,0.386603,0.469580,0.505383};
-// PROTONS
-const Int_t kPBinsPbPb = 9;
-const Double_t kPBinWidthPbPb = 0.500;
-const Double_t kPTagPPbPb[kPBinsPbPb] = {0.017940,0.350681,0.535286,0.583264,0.562935,0.560524,0.545992,0.598060,0.351245};
-const Double_t kPTagPiPbPb[kPBinsPbPb] = {0.195955,0.094949,0.039962,0.026039,0.007556,0.016986,0.030333,0.000000,0.000000};
-const Double_t kPTagNidPbPb[kPBinsPbPb] = {0.786105,0.554370,0.424751,0.390697,0.429508,0.422491,0.423675,0.401940,0.648755};
-//
-// pp PYTHIA
-//
-// *** cuts for pp ***
-//
-// PIONS
-const Int_t kPiBinsPP = 10;
-const Double_t kPiBinWidthPP = 0.250;
-const Double_t kPiTagPiPP[kPiBinsPP] = {0.194528,0.447097,0.603364,0.646413,0.647125,0.669157,0.688139,0.682564,0.689910,0.665710};
-const Double_t kPiTagNidPP[kPiBinsPP] = {0.805472,0.552903,0.396636,0.353587,0.352875,0.330843,0.311861,0.317436,0.310090,0.334290};
-// KAONS
-const Int_t kKBinsPP = 10;
-const Double_t kKBinWidthPP = 0.250;
-const Double_t kKTagKPP[kKBinsPP] = {0.000000,0.173393,0.439690,0.519423,0.587025,0.605372,0.586021,0.650139,0.444444,0.299363};
-const Double_t kKTagPiPP[kKBinsPP] = {0.000000,0.001495,0.000000,-0.000000,-0.000000,0.000000,0.032258,0.060572,0.101449,0.242038};
-const Double_t kKTagNidPP[kKBinsPP] = {1.000000,0.825112,0.560310,0.480577,0.412975,0.394628,0.381720,0.289289,0.454106,0.458599};
-// PROTONS
-const Int_t kPBinsPP = 9;
-const Double_t kPBinWidthPP = 0.500;
-const Double_t kPTagPPP[kPBinsPP] = {0.029404,0.438640,0.613710,0.665152,0.634961,0.657711,0.703704,0.685714,0.235294};
-const Double_t kPTagPiPP[kPBinsPP] = {0.000000,0.000000,0.000000,-0.000000,0.000000,0.000000,-0.000000,0.014286,-0.000000};
-const Double_t kPTagNidPP[kPBinsPP] = {0.970596,0.561360,0.386290,0.334848,0.365039,0.342289,0.296296,0.300000,0.764706};
-
-
-
-
-//-----------------------------------------------------------------------------
-class AliD0toKpi : public TObject {
- public:
- //
- AliD0toKpi();
- AliD0toKpi(Int_t ev,Int_t trkNum[2],
- Double_t v1[3],Double_t v2[3],Double_t dca,
- Double_t mom[6],Double_t d0[2]);
- virtual ~AliD0toKpi();
- AliD0toKpi(const AliD0toKpi& d0toKpi);
-
- Double_t Alpha() const { return (Ql(0)-Ql(1))/(Ql(0)+Ql(1)); }
- void ApplyPID(TString pidScheme="TOFparamPbPb");
- Double_t ChildrenRelAngle() const;
- void ComputeWgts();
- void CorrectWgt4BR(Double_t factor);
- Double_t CosPointing() const;
- Double_t CosPointingXY() const;
- void CosThetaStar(Double_t &ctsD0,Double_t &ctsD0bar) const;
- Double_t Ct() const {return Length()*kMD0/P();}
- Double_t Energy() const { return TMath::Sqrt(P()*P()+kMD0*kMD0); }
- Double_t Eta() const;
- Double_t EtaChild(Int_t child) const;
- Int_t EventNo() const {return TMath::Abs(fEvent);}
- Double_t GetDCA() const { return 10000.*fDCA; }
- Int_t GetTrkNum(Int_t child) const { return fTrkNum[child]; }
- Double_t Getd0Child(Int_t child) const { return fd0[child]; }
- Int_t GetPdgChild(Int_t child) const { return fPdg[child]; }
- Int_t GetPdgMum(Int_t child) const {return fMum[child]; }
- void GetWgts(Double_t &WgtD0,Double_t &WgtD0bar,TString sample) const;
- void GetPrimaryVtx(Double_t vtx[3]) const
- { vtx[0]=fV1x; vtx[1]=fV1y; vtx[2]=fV1z; }
- void GetSecondaryVtx(Double_t vtx[3]) const
- { vtx[0]=fV2x; vtx[1]=fV2y; vtx[2]=fV2z; }
-
- Double_t ImpPar() const;
- void InvMass(Double_t &mD0,Double_t &mD0bar) const;
- Bool_t IsSignal() const { if(fSignal) return kTRUE; return kFALSE; }
- Double_t Length() const
- { return TMath::Sqrt((fV1x-fV2x)*(fV1x-fV2x)
- +(fV1y-fV2y)*(fV1y-fV2y)+(fV1z-fV2z)*(fV1z-fV2z)); }
- Double_t P() const { return TMath::Sqrt(Pt()*Pt()+Pz()*Pz()); }
- Double_t PChild(Int_t child) const { return TMath::Sqrt(fPx[child]*fPx[child]+fPy[child]*fPy[child]+fPz[child]*fPz[child]); }
- Double_t ProdImpParams() const { return fd0[0]*fd0[1]; }
- Double_t Pt() const { return TMath::Sqrt(Px()*Px()+Py()*Py()); }
- Double_t PtChild(Int_t child) const { return TMath::Sqrt(fPx[child]*fPx[child]+fPy[child]*fPy[child]); }
- Double_t Px() const { return (fPx[0]+fPx[1]); }
- Double_t Py() const { return (fPy[0]+fPy[1]); }
- Double_t Pz() const { return (fPz[0]+fPz[1]); }
- Double_t Ql(Int_t child) const;
- Double_t Qt() const;
- Double_t Rapidity() const { return 0.5*TMath::Log((Energy()+Pz())/(Energy()-Pz()+1.e-13)); }
- Bool_t Select(const Double_t* cuts,Int_t& okD0,Int_t& okD0bar) const;
- void SetPrimaryVtx(Double_t vtx[3])
- { fV1x=vtx[0]; fV1y=vtx[1]; fV1z=vtx[2]; }
- void SetSignal() { fSignal = kTRUE; }
- void SetTOFmasses(Double_t mass[2])
- { fTOFmass[0]=mass[0]; fTOFmass[1]=mass[1]; }
- void SetPIDresponse(Double_t resp0[5],Double_t resp1[5]);
- void SetPdgCodes(Int_t pdg[2]) {fPdg[0]=pdg[0];fPdg[1]=pdg[1]; }
- void SetMumPdgCodes(Int_t mum[2]) {fMum[0]=mum[0];fMum[1]=mum[1]; }
- Double_t LinearInterpolation(Double_t p,Int_t nBins,Double_t Bin,
- const Double_t *values) const;
- //
- private:
- //
- Bool_t fSignal; // TRUE if signal, FALSE if background (for simulation)
- Int_t fEvent; // number of the event this D0 comes from
- // -1 if the D0 comes from ev. mixing
-
- Int_t fTrkNum[2]; // numbers of the two decay tracks
-
- Double_t fV1x; // X-position of the primary vertex of the event
- Double_t fV1y; // Y-position of the primary vertex of the event
- Double_t fV1z; // Z-position of the primary vertex of the event
- Double_t fV2x; // X-position of the reconstructed secondary vertex
- Double_t fV2y; // Y-position of the reconstructed secondary vertex
- Double_t fV2z; // Z-position of the reconstructed secondary vertex
- Double_t fDCA; // DCA of the two tracks
-
- Double_t fPx[2]; // X,Y,Z
- Double_t fPy[2]; // momenta of the two tracks
- Double_t fPz[2]; // at the reconstructed vertex
-
- Double_t fd0[2]; // impact parameters in the bending plane
-
- Int_t fPdg[2]; // PDG codes of the two tracks (for sim.)
- Int_t fMum[2]; // PDG codes of the mothers (for sim.)
-
- Double_t fTagPi[2]; // probability to be tagged as pion
- Double_t fTagKa[2]; // probability to be tagged as kaon
- Double_t fTagPr[2]; // probability to be tagged as proton
- Double_t fTagNid[2]; // probability to be tagged as "non-identified"
-
- Double_t fPIDrespEl[2]; // det. response to be electron
- Double_t fPIDrespMu[2]; // det. response to be muon
- Double_t fPIDrespPi[2]; // det. response to be pion
- Double_t fPIDrespKa[2]; // det. response to be kaon
- Double_t fPIDrespPr[2]; // det. response to be proton
- Double_t fTOFmass[2]; // mass estimated by the TOF (-1000. if track not reached TOF)
-
- Double_t fWgtAD0,fWgtAD0bar; // weights for the 3 samples
- Double_t fWgtBD0,fWgtBD0bar; // weights for the 3 samples
- Double_t fWgtCD0,fWgtCD0bar; // A: (K,Pi)+(K,?) B: (?,Pi) C: (?,?)
- Double_t fWgtDD0,fWgtDD0bar; // D: all other pairs
-
- ClassDef(AliD0toKpi,1) // Reconstructed D0 candidate class
-};
-
-#endif
-
-
-
-
-
-
-
-
+++ /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)
-// An example of usage can be found in the macro AliD0toKpiTest.C
-// Origin: A. Dainese andrea.dainese@lnl.infn.it
-//----------------------------------------------------------------------------
-#include <TKey.h>
-#include <TFile.h>
-#include <TTree.h>
-#include <TString.h>
-#include <TSystem.h>
-#include <TParticle.h>
-#include "AliESD.h"
-#include "AliMC.h"
-#include "AliRun.h"
-#include "AliRunLoader.h"
-#include "AliVertexerTracks.h"
-#include "AliESDVertex.h"
-#include "AliESDv0.h"
-#include "AliD0toKpi.h"
-#include "AliD0toKpiAnalysis.h"
-#include "AliLog.h"
-
-typedef struct {
- Int_t lab;
- Int_t pdg;
- Int_t mumlab;
- Int_t mumpdg;
- Int_t mumprongs;
- Float_t Vx,Vy,Vz;
- Float_t Px,Py,Pz;
-} REFTRACK;
-
-ClassImp(AliD0toKpiAnalysis)
-
-//----------------------------------------------------------------------------
-AliD0toKpiAnalysis::AliD0toKpiAnalysis() {
- // Default constructor
-
- SetPtCut();
- Setd0Cut();
- SetMassCut();
- SetD0Cuts();
- SetVertex1();
- SetPID();
- fVertexOnTheFly = kFALSE;
- fSim = kFALSE;
- fOnlySignal = kFALSE;
-}
-//----------------------------------------------------------------------------
-AliD0toKpiAnalysis::~AliD0toKpiAnalysis() {}
-//----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::ApplySelection(const Char_t *inName,const Char_t *outName) const {
- // select candidates that pass fD0Cuts and write them to a new file
-
- TFile *inFile = TFile::Open(inName);
-
- TTree *treeD0in=(TTree*)inFile->Get("TreeD0");
- AliD0toKpiAnalysis *inAnalysis = (AliD0toKpiAnalysis*)inFile->Get("D0toKpiAnalysis");
- printf("+++\n+++ I N P U T S T A T U S:\n+++\n");
- inAnalysis->PrintStatus();
-
-
- AliD0toKpi *d = 0;
- treeD0in->SetBranchAddress("D0toKpi",&d);
- Int_t entries = (Int_t)treeD0in->GetEntries();
-
- printf("+++\n+++ Number of D0 in input tree: %d\n+++\n",entries);
-
- TTree *treeD0out = new TTree("TreeD0","Tree with selected D0 candidates");
- treeD0out->Branch("D0toKpi","AliD0toKpi",&d,200000,0);
-
-
- Int_t okD0=0,okD0bar=0;
- Int_t nSel = 0;
-
- for(Int_t i=0; 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
-
-
- TString esdName="AliESDs.root";
- if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
- printf("AliD0toKpiAnalysis::FindCandidatesESD(): No ESDs file found!\n");
- return;
- }
-
- TString outName1=outName;
- TString outName2="nTotEvents.dat";
-
- Int_t nTotEv=0,nD0rec=0,nD0rec1ev=0;
- Double_t dca;
- Double_t v2[3],mom[6],d0[2];
- Int_t iTrkP,iTrkN,trkEntries;
- Int_t nTrksP=0,nTrksN=0;
- Int_t trkNum[2];
- Double_t tofmass[2];
- Int_t okD0=0,okD0bar=0;
- AliESDtrack *postrack = 0;
- AliESDtrack *negtrack = 0;
-
- // create the AliVertexerTracks object
- // (it will be used only if fVertexOnTheFly=kTrue)
- AliVertexerTracks *vertexer1 = new AliVertexerTracks;
- if(fVertexOnTheFly) {
- // open the mean vertex
- TFile *invtx = new TFile("AliESDVertexMean.root");
- AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean");
- invtx->Close();
- vertexer1->SetVtxStart(initVertex);
- delete invtx;
- }
- Int_t skipped[2];
- Bool_t goodVtx1;
-
- // create tree for reconstructed D0s
- AliD0toKpi *ioD0toKpi=0;
- TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates");
- treeD0->Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
-
- // open file with tracks
- TFile *esdFile = TFile::Open(esdName.Data());
- AliESD* event = new AliESD;
- TTree* tree = (TTree*) esdFile->Get("esdTree");
- if(!tree) {
- Error("FindCandidatesESD", "no ESD tree found");
- return;
- }
- tree->SetBranchAddress("ESD",&event);
-
- // loop on events in file
- for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
- if(iEvent > evLast) break;
- tree->GetEvent(iEvent);
- Int_t ev = (Int_t)event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
- printf("--- Finding D0 -> Kpi in event %d\n",ev);
- // count the total number of events
- nTotEv++;
-
- trkEntries = (Int_t)event->GetNumberOfTracks();
- printf(" Number of tracks: %d\n",trkEntries);
- if(trkEntries<2) continue;
-
- // retrieve primary vertex from file
- if(!event->GetPrimaryVertex()) {
- printf(" No vertex\n");
- continue;
- }
- event->GetPrimaryVertex()->GetXYZ(fV1);
-
- // call function which applies sigle-track selection and
- // separetes positives and negatives
- TObjArray trksP(trkEntries/2);
- Int_t *trkEntryP = new Int_t[trkEntries];
- TObjArray trksN(trkEntries/2);
- Int_t *trkEntryN = new Int_t[trkEntries];
- TTree *trkTree = new TTree();
- SelectTracks(event,trksP,trkEntryP,nTrksP,
- trksN,trkEntryN,nTrksN);
-
- printf(" pos. tracks: %d neg .tracks: %d\n",nTrksP,nTrksN);
-
-
- nD0rec1ev = 0;
-
- // LOOP ON POSITIVE TRACKS
- for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
- if(iTrkP%1==0) printf(" Processing positive track number %d of %d\n",iTrkP,nTrksP);
-
- // get track from track array
- postrack = (AliESDtrack*)trksP.UncheckedAt(iTrkP);
- trkNum[0] = trkEntryP[iTrkP];
-
- // LOOP ON NEGATIVE TRACKS
- for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
-
- // get track from tracks array
- negtrack = (AliESDtrack*)trksN.UncheckedAt(iTrkN);
- trkNum[1] = trkEntryN[iTrkN];
-
- //
- // ----------- DCA MINIMIZATION ------------------
- //
- // find the DCA and propagate the tracks to the DCA
- Double_t b=event->GetMagneticField();
- AliESDtrack nt(*negtrack), pt(*postrack);
- dca = nt.PropagateToDCA(&pt,b);
-
- // define the AliESDv0 object
- AliESDv0 vertex2(nt,trkNum[0],pt,trkNum[1]);
-
- // get position of the secondary vertex
- vertex2.GetXYZ(v2[0],v2[1],v2[2]);
- vertex2.GetPPxPyPz(mom[0],mom[1],mom[2]);
- vertex2.GetNPxPyPz(mom[3],mom[4],mom[5]);
- // impact parameters of the tracks w.r.t. the primary vertex
- d0[0] = 10000.*pt.GetD(fV1[0],fV1[1],b);
- d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b);
- goodVtx1 = kTRUE;
-
- // no vertexing if DeltaMass > fMassCut
- if(fVertexOnTheFly) {
- goodVtx1 = kFALSE;
- if(SelectInvMass(mom)) {
- // primary vertex from *other* tracks in the event
- skipped[0] = trkEntryP[iTrkP];
- skipped[1] = trkEntryN[iTrkN];
- vertexer1->SetSkipTracks(2,skipped);
- AliESDVertex *vertex1onfly =
- (AliESDVertex*)vertexer1->FindPrimaryVertex(event);
- if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
- vertex1onfly->GetXYZ(fV1);
- //vertex1onfly->PrintStatus();
- delete vertex1onfly;
- }
- }
-
-
- // create the object AliD0toKpi
- AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
- // select D0s
- if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
- // get PID info from ESD
- AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
- Double_t esdpid0[5];
- t0->GetESDpid(esdpid0);
- if(t0->GetStatus()&AliESDtrack::kTOFpid) {
- tofmass[0] = CalculateTOFmass(t0->GetP(),
- t0->GetIntegratedLength(),
- t0->GetTOFsignal());
- } else {
- tofmass[0] = -1000.;
- }
- AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]);
- Double_t esdpid1[5];
- t1->GetESDpid(esdpid1);
- if(t1->GetStatus()&AliESDtrack::kTOFpid) {
- tofmass[1] = CalculateTOFmass(t1->GetP(),
- t1->GetIntegratedLength(),
- t1->GetTOFsignal());
- } else {
- tofmass[1] = -1000.;
- }
-
- theD0.SetPIDresponse(esdpid0,esdpid1);
- theD0.SetTOFmasses(tofmass);
-
- // fill the tree
- ioD0toKpi=&theD0;
- treeD0->Fill();
-
- nD0rec++; nD0rec1ev++;
- ioD0toKpi=0;
- }
-
- negtrack = 0;
- } // loop on negative tracks
- postrack = 0;
- } // loop on positive tracks
-
- delete [] trkEntryP;
- delete [] trkEntryN;
- delete trkTree;
-
- printf(" Number of D0 candidates: %d\n",nD0rec1ev);
- } // loop on events in file
-
-
- printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
- printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
-
- delete vertexer1;
-
- esdFile->Close();
-
- // create a copy of this class to be written to output file
- AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis");
-
- // add PDG codes to decay tracks in found candidates (in simulation mode)
- // and store tree in the output file
- if(!fSim) {
- TFile *outroot = new TFile(outName1.Data(),"recreate");
- treeD0->Write();
- copy->Write();
- outroot->Close();
- delete outroot;
- } else {
- printf(" Now adding information from simulation (PDG codes) ...\n");
- TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates");
- SimulationInfo(treeD0,treeD0sim);
- delete treeD0;
- TFile *outroot = new TFile(outName1.Data(),"recreate");
- treeD0sim->Write();
- copy->Write();
- outroot->Close();
- delete outroot;
- }
-
- // write to a file the total number of events
- FILE *outdat = fopen(outName2.Data(),"w");
- fprintf(outdat,"%d\n",nTotEv);
- fclose(outdat);
-
- return;
-}
-//-----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::PrintStatus() const {
- // Print parameters being used
-
- printf("Preselections:\n");
- printf(" fPtCut = %f GeV\n",fPtCut);
- printf(" fd0Cut = %f micron\n",fd0Cut);
- printf(" fMassCut = %f GeV\n",fMassCut);
- if(fVertexOnTheFly) printf("Primary vertex on the fly\n");
- if(fSim) {
- printf("Simulation mode\n");
- if(fOnlySignal) printf(" Only signal goes to file\n");
- }
- printf("Cuts on candidates:\n");
- printf(" |M-MD0| [GeV] < %f\n",fD0Cuts[0]);
- printf(" dca [micron] < %f\n",fD0Cuts[1]);
- printf(" cosThetaStar < %f\n",fD0Cuts[2]);
- printf(" pTK [GeV] > %f\n",fD0Cuts[3]);
- printf(" pTpi [GeV] > %f\n",fD0Cuts[4]);
- printf(" |d0K| [micron] < %f\n",fD0Cuts[5]);
- printf(" |d0pi| [micron] < %f\n",fD0Cuts[6]);
- printf(" d0d0 [micron^2] < %f\n",fD0Cuts[7]);
- printf(" cosThetaPoint > %f\n",fD0Cuts[8]);
-
- return;
-}
-//-----------------------------------------------------------------------------
-Bool_t AliD0toKpiAnalysis::SelectInvMass(const Double_t p[6]) const {
- // Apply preselection in the invariant mass of the pair
-
- Double_t mD0 = 1.8645;
- Double_t mPi = 0.13957;
- Double_t mKa = 0.49368;
-
- Double_t energy[2];
- Double_t mom2[2],momTot2;
-
- mom2[0] = p[0]*p[0] + p[1]*p[1] + p[2]*p[2];
- mom2[1] = p[3]*p[3] + p[4]*p[4] + p[5]*p[5];
-
- momTot2 = (p[0]+p[3])*(p[0]+p[3])+
- (p[1]+p[4])*(p[1]+p[4])+
- (p[2]+p[5])*(p[2]+p[5]);
-
- // D0 -> K- pi+
- energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]);
- energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]);
-
- Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
-
- // D0bar -> K+ pi-
- energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]);
- energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]);
-
- Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
-
- if(TMath::Abs(minvD0-mD0) < fMassCut) return kTRUE;
- if(TMath::Abs(minvD0bar-mD0) < fMassCut) return kTRUE;
- return kFALSE;
-}
-//-----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::SelectTracks(AliESD *event,
- TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
- TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
- // Create two TObjArrays with positive and negative tracks and
- // apply single-track preselection
-
- nTrksP=0,nTrksN=0;
-
- Int_t entr = event->GetNumberOfTracks();
-
- // transfer ITS tracks from ESD to arrays and to a tree
- for(Int_t i=0; i<entr; i++) {
-
- AliESDtrack *esdtrack = event->GetTrack(i);
- UInt_t status = esdtrack->GetStatus();
-
- if(!(status&AliESDtrack::kITSin)) continue;
-
- // single track selection
- if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue;
-
- if(esdtrack->GetSign()<0) { // negative track
- trksN.AddLast(esdtrack);
- trkEntryN[nTrksN] = i;
- nTrksN++;
- } else { // positive track
- trksP.AddLast(esdtrack);
- trkEntryP[nTrksP] = i;
- nTrksP++;
- }
-
- } // loop on ESD tracks
-
- return;
-}
-//-----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::SetD0Cuts(Double_t cut0,Double_t cut1,
- Double_t cut2,Double_t cut3,Double_t cut4,
- Double_t cut5,Double_t cut6,
- Double_t cut7,Double_t cut8) {
- // Set the cuts for D0 selection
- fD0Cuts[0] = cut0;
- fD0Cuts[1] = cut1;
- fD0Cuts[2] = cut2;
- fD0Cuts[3] = cut3;
- fD0Cuts[4] = cut4;
- fD0Cuts[5] = cut5;
- fD0Cuts[6] = cut6;
- fD0Cuts[7] = cut7;
- fD0Cuts[8] = cut8;
-
- return;
-}
-//-----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::SetD0Cuts(const Double_t cuts[9]) {
- // Set the cuts for D0 selection
-
- for(Int_t i=0; i<9; i++) fD0Cuts[i] = cuts[i];
-
- return;
-}
-//-----------------------------------------------------------------------------
-Bool_t
-AliD0toKpiAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const {
- // Check if track passes some kinematical cuts
- // Magnetic field "b" (kG)
-
- if(TMath::Abs(1./trk.GetParameter()[4]) < fPtCut)
- return kFALSE;
- if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1],b)) < fd0Cut)
- return kFALSE;
-
- return kTRUE;
-}
-//----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::MakeTracksRefFile(AliRun *gAlice,
- Int_t evFirst,Int_t evLast) const {
- // Create a file with simulation info for the reconstructed tracks
-
- TFile *outFile = TFile::Open("D0TracksRefFile.root","recreate");
- TFile *esdFile = TFile::Open("AliESDs.root");
-
- AliMC *mc = gAlice->GetMCApp();
-
- Int_t label;
- TParticle *part;
- TParticle *mumpart;
- REFTRACK reftrk;
-
- AliESD* event = new AliESD;
- TTree* tree = (TTree*) esdFile->Get("esdTree");
- tree->SetBranchAddress("ESD",&event);
- // loop on events in file
- for(Int_t iEvent=evFirst; iEvent<tree->GetEntries(); iEvent++) {
- if(iEvent>evLast) break;
- tree->GetEvent(iEvent);
- Int_t ev = (Int_t)event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
-
- gAlice->GetEvent(ev);
-
- Int_t nentr=(Int_t)event->GetNumberOfTracks();
-
- // Tree for true track parameters
- char ttname[100];
- sprintf(ttname,"Tree_Ref_%d",ev);
- TTree *reftree = new TTree(ttname,"Tree with true track params");
- reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
-
- for(Int_t i=0; i<nentr; i++) {
- AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
- label = TMath::Abs(esdtrack->GetLabel());
-
- part = (TParticle*)mc->Particle(label);
- reftrk.lab = label;
- reftrk.pdg = part->GetPdgCode();
- reftrk.mumlab = part->GetFirstMother();
- if(part->GetFirstMother()>=0) {
- mumpart = (TParticle*)gAlice->GetMCApp()->Particle(part->GetFirstMother());
- reftrk.mumpdg = mumpart->GetPdgCode();
- reftrk.mumprongs = mumpart->GetNDaughters();
- } else {
- reftrk.mumpdg=-1;
- reftrk.mumprongs=-1;
- }
- reftrk.Vx = part->Vx();
- reftrk.Vy = part->Vy();
- reftrk.Vz = part->Vz();
- reftrk.Px = part->Px();
- reftrk.Py = part->Py();
- reftrk.Pz = part->Pz();
-
- reftree->Fill();
-
- } // loop on tracks
-
- outFile->cd();
- reftree->Write();
-
- delete reftree;
- } // loop on events
-
- esdFile->Close();
- outFile->Close();
-
- return;
-}
-//-----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::SimulationInfo(TTree *treeD0in,TTree *treeD0out) const {
- // add pdg codes to candidate decay tracks (for sim)
-
- TString refFileName("D0TracksRefFile.root");
- if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) {
- printf("AliD0toKpiAnalysis::SimulationInfo: no reference file found!\n");
- return;
- }
- TFile *refFile = TFile::Open(refFileName.Data());
-
- Char_t refTreeName[100];
- Int_t event;
- Int_t pdg[2],mumpdg[2],mumlab[2];
- REFTRACK reftrk;
-
- // read-in reference tree for event 0 (the only event for Pb-Pb)
- sprintf(refTreeName,"Tree_Ref_%d",0);
- TTree *refTree0 = (TTree*)refFile->Get(refTreeName);
- refTree0->SetBranchAddress("rectracks",&reftrk);
-
- AliD0toKpi *theD0 = 0;
- treeD0in->SetBranchAddress("D0toKpi",&theD0);
- treeD0out->Branch("D0toKpi","AliD0toKpi",&theD0,200000,0);
-
- Int_t entries = (Int_t)treeD0in->GetEntries();
-
- for(Int_t i=0; 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] &&
- reftrk.mumprongs==2 &&
- ((TMath::Abs(pdg[0])==211 && TMath::Abs(pdg[1])==321) ||
- (TMath::Abs(pdg[0])==321 && TMath::Abs(pdg[1])==211))
- ) theD0->SetSignal();
-
- if(!fOnlySignal || theD0->IsSignal()) treeD0out->Fill();
-
- }
-
- delete refTree0;
-
- refFile->Close();
-
- return;
-}
-
-
-
-
-
-
+++ /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@lnl.infn.it
-//-------------------------------------------------------------------------
-
-#include <TString.h>
-#include <TNamed.h>
-#include "AliESD.h"
-#include "AliRun.h"
-
-//-----------------------------------------------------------------------------
-class AliD0toKpiAnalysis : public TNamed {
- public:
- //
- AliD0toKpiAnalysis();
- virtual ~AliD0toKpiAnalysis();
-
- void ApplySelection(const Char_t *inName="AliD0toKpi.root",
- const Char_t *outName="AliD0toKpi_sele.root") const;
- void FindCandidates(Int_t evFirst=0,Int_t evLast=0,
- const Char_t *outName="AliD0toKpi.root");
- void MakeTracksRefFile(AliRun *gAlice,Int_t evFirst=0,Int_t evLast=0) const;
- void PrintStatus() const;
- void SetVertexOnTheFly() { fVertexOnTheFly=kTRUE; }
- void SetSimulation() { fSim=kTRUE; }
- void SetOnlySignal() { fOnlySignal=kTRUE; }
- void SetPtCut(Double_t pt=0.) { fPtCut=pt; }
- void Setd0Cut(Double_t d0=0.) { fd0Cut=d0; }
- void SetMassCut(Double_t deltaM=1000.) { fMassCut=deltaM; }
- void SetD0Cuts(Double_t cut0=1000.,Double_t cut1=100000.,
- Double_t cut2=1.1,Double_t cut3=0.,Double_t cut4=0.,
- Double_t cut5=100000.,Double_t cut6=100000.,
- Double_t cut7=100000000.,Double_t cut8=-1.1);
- void SetD0Cuts(const Double_t cuts[9]);
- void SetPID(const Char_t * pid="TOFparam_PbPb") { fPID=pid; }
- //
- private:
- //
- Bool_t fVertexOnTheFly; // flag for primary vertex reco on the fly
- Bool_t fSim; // flag for the analysis of simulated events
- Bool_t fOnlySignal; // write to file only signal candidates (for sim)
- TString fPID; // PID scheme
-
- Double_t fV1[3]; // primary vertex position (in cm)
- Double_t fPtCut; // minimum track pt (in GeV/c)
- Double_t fd0Cut; // minimum track |rphi impact parameter| (in micron)
- Double_t fMassCut; // maximum of |InvMass-MD0| (in GeV)
- Double_t fD0Cuts[9]; // cuts on D0 candidates (see SetD0Cuts())
- // (to be passed to function AliD0toKpi::Select())
- // 0 = inv. mass half width [GeV]
- // 1 = dca [micron]
- // 2 = cosThetaStar
- // 3 = pTK [GeV/c]
- // 4 = pTPi [GeV/c]
- // 5 = d0K [micron] upper limit!
- // 6 = d0Pi [micron] upper limit!
- // 7 = d0d0 [micron^2]
- // 8 = cosThetaPoint
-
- //
- Double_t CalculateTOFmass(Double_t mom,Double_t length,Double_t time) const;
- Bool_t SelectInvMass(const Double_t p[6]) const;
- void SelectTracks(AliESD *event,
- TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
- TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const;
- void SetVertex1(Double_t x=0.,Double_t y=0.,Double_t z=0.)
- { fV1[0]=x;fV1[1]=y;fV1[2]=z; }
- void SimulationInfo(TTree *treeD0in,TTree *treeD0out) const;
- Bool_t SingleTrkCuts(const AliESDtrack& trk, Double_t b) const;
- //
- ClassDef(AliD0toKpiAnalysis,3) // Reconstruction of D0 candidates class
-};
-
-
-#endif
-
-
-
-
-
-
-
-
+++ /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@lnl.infn.it
- //--------------------------------------------------------------------------
-
- gSystem->Load("libANALYSIS.so");
-
- // set of cuts
- Double_t D0Cuts[9] = {0.1, // mass [GeV]
- 1000000., // dca [micron]
- 1.1, // cosThetaStar
- 0., // pT K [GeV/c]
- 0., // pT Pi [GeV/c]
- 100000., // d0K upper [micron]
- 100000., // d0Pi upper [micron]
- 10000000000., // d0d0 [micron^2]
- -1.1}; // cosThetaPointing
-
- // number of events (for normalization)
- Bool_t normalize = kFALSE;
- Double_t events = 1.;
-
-
-
- // define histograms
- TH1F *hptK = new TH1F("hptK","\"K\" p_{t} distribution",50,0,10);
- hptK->SetXTitle("p_{t} [GeV]");
-
- TH1F *hptPi = new TH1F("hptPi","\"#pi\" p_{t} distribution",50,0,10);
- hptPi->SetXTitle("p_{t} [GeV]");
-
- TH1F *hDCA = new TH1F("hDCA","DCA",50,0,1000);
- hDCA->SetXTitle("dca [#mu m]");
-
- TH1F *hptD0 = new TH1F("hptD0","D^{0} p_{t} distribution",40,0,40);
- hptD0->SetXTitle("p_{t} [GeV]");
-
- TH1F *hyD0 = new TH1F("hyD0","D^{0} rapidity distribution",50,-2,2);
- hyD0->SetXTitle("y");
-
- TH1F *hCPtaD0 = new TH1F("hCPtaD0","cosine of pointing angle distribution",100,-1,1);
- hCPtaD0->SetXTitle("cos #theta_{point}");
-
- TH1F *hCPtaXY = new TH1F("hCPtaXY","cosine of pointing angle in (x,y) plane",100,-1,1);
- hCPtaXY->SetXTitle("cos #theta_{point}");
-
- TH1F *hCts = new TH1F("hCts","cosine of decay angle",50,-1.2,1.2);
- hCts->SetXTitle("cos #theta^{*}");
-
- TH2F *hCtsVsPtK = new TH2F("hCtsVsPtK","cosine of decay angle VS \"K\" p_{t}",50,0,5,50,-1,1);
- hCtsVsPtK->SetYTitle("cos #theta^{*}");
- hCtsVsPtK->SetXTitle("p_{t} [GeV]");
-
- TH1F *hd0d0 = new TH1F("hd0d0","Product of the impact parameters",100,-100000,100000);
- hd0d0->SetXTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]");
-
- TH1F *hd0K = new TH1F("hd0K","Impact parameter of \"K\"",100,-5000,5000);
- hd0K->SetXTitle("d_{0}^{K} [#mu m]");
-
- TH1F *hd0Pi = new TH1F("hd0Pi","Impact parameter of \"#pi\"",100,-5000,5000);
- hd0Pi->SetXTitle("d_{0}^{#pi} [#mu m]");
-
- TH2F *hCPtaVsd0d0 = new TH2F("hCPtaVsd0d0","cos #theta_{point} vs d_{0}^{K} #times d_{0}^{#pi}",100,-100000,100000,100,-1,1);
- hCPtaVsd0d0->SetXTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]");
- hCPtaVsd0d0->SetYTitle("cos #theta_{point}");
-
- TH2F *hCPtaVsd0d0zoom = new TH2F("hCPtaVsd0d0zoom","cos #theta_{point} vs d_{0}^{K} #times d_{0}^{#pi}",100,-100000,0,100,.9,1);
- hCPtaVsd0d0zoom->SetXTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]");
- hCPtaVsd0d0zoom->SetYTitle("cos #theta_{point}");
-
- TH2F *hd0d0VSptD0 = new TH2F("hd0d0VSptD0","d_{0}^{K} #times d_{0}^{#pi} VS D^{0} p_{t}",50,0,25,100,-120000,120000);
- hd0d0VSptD0->SetYTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]");
- hd0d0VSptD0->SetXTitle("D^{0} p_{t} [GeV]");
-
- TH1F *hMass = new TH1F("hMass","Invariant mass distribution",50,1.765,1.965);
- hMass->SetXTitle("M[K,#pi] [GeV]");
-
- TH2F *hArm = new TH2F("hArm","Armenteros plot",50,-2,2,50,0,1);
- hArm->SetXTitle("#alpha");
- hArm->SetYTitle("q_{t}");
-
- // open input file and get tree
- TFile *inFile = TFile::Open(inName);
-
- TTree *treeD0 = (TTree*)inFile->Get("TreeD0");
- AliD0toKpi *D = 0;
- treeD0->SetBranchAddress("D0toKpi",&D);
- Int_t entries = (Int_t)treeD0->GetEntries();
-
- printf("+++\n+++ Number of D0 in tree: %d\n+++\n",entries);
-
- Double_t MD0,MD0bar,ctsD0,ctsD0bar,ctsPiD0,ctsPiD0bar;
- Double_t WgtD0,WgtD0bar;
- Double_t sampleABC=0.;
- Int_t okD0=0,okD0bar=0;
- Int_t nSel = 0;
- Int_t ptbin;
-
- // loop on D0
- for(Int_t i=0; 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@lnl.infn.it
-//--------------------------------------------------------------------------
-
-void AliD0toKpiReco() {
-
- gSystem->Load("libANALYSIS.so");
-
- //============== R E C O N S T R U C T I O N ==============================
-
- Int_t evFirst = 0;
- Int_t evLast = 1000000;
-
- // Get field from galice.root
- if (gAlice) {
- delete gAlice->GetRunLoader();
- delete gAlice;
- gAlice=0;
- }
- AliRunLoader *rl = AliRunLoader::Open("galice.root");
- if (rl == 0x0) {
- cerr<<"Can not open session"<<endl;
- return;
- }
- Int_t retval = rl->LoadgAlice();
- if (retval) {
- cerr<<"LoadgAlice returned error"<<endl;
- delete rl;
- return;
- }
- gAlice=rl->GetAliRun();
- AliMagF *fiel = (AliMagF*)gAlice->Field();
- // Set the conversion constant between curvature and Pt
- AliTracker::SetFieldMap(fiel,kTRUE);
-
- AliD0toKpiAnalysis *analysis = new AliD0toKpiAnalysis();
- // set simulation to take info on PDG codes from kine tree
- analysis->SetSimulation();
- rl->LoadKinematics();
- analysis->MakeTracksRefFile(gAlice,evFirst,evLast);
- //--- 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.); // GeV
- analysis->Setd0Cut(0.); // 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();
-
- analysis->FindCandidates(evFirst,evLast,"AliD0toKpi.root");
- delete analysis;
-
- return;
-}
-//==========================================================================
-void AliD0toKpiSele() {
-
- gSystem->Load("libANALYSIS.so");
-
- //======================== 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;
-}
-//==========================================================================
-
-
-
-
-
-
#pragma link C++ class TGliteXmlEventlist+;
-#pragma link C++ class AliD0toKpi+;
-#pragma link C++ class AliD0toKpiAnalysis+;
#pragma link C++ class AliRunAnalysis+;
#pragma link C++ class AliAnalysis+;
AliReader.cxx AliReaderAOD.cxx AliReaderKineTree.cxx \
AliReaderESD.cxx AliReaderESDTree.cxx \
AliTrackPoints.cxx AliClusterMap.cxx \
- AliD0toKpi.cxx AliD0toKpiAnalysis.cxx AliFlowAnalysis.cxx \
+ AliFlowAnalysis.cxx \
AliMuonAnalysis.cxx AliAODv0.cxx AliAODxi.cxx AliAODevent.cxx \
AliAnalysisEventCuts.cxx AliAnalysisTrackCuts.cxx
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff