Page MenuHomeHEPForge

PhaseSpacePlots.C
No OneTemporary

Size
3 KB
Referenced Files
None
Subscribers
None

PhaseSpacePlots.C

#include <string>
using std::string;
void PhaseSpacePlots() {
string fileName1 = "rootFiles/B0JpsiKpi.root";
string fileName2 = "rootFiles/B0KpiJpsi.root";
string decay1 = "B^{0} #rightarrow J/#psi K #pi";
string decay2 = "B^{0} #rightarrow K #pi J/#psi";
string xLabel = "m^{2}(K#pi)";
string yLabel = "m^{2}(J/#psi #pi)";
int xInt1 = 1;
int yInt1 = 2;
int xInt2 = 3;
int yInt2 = 1;
TLatex latex;
latex.SetTextSize(0.05);
latex.SetNDC(kTRUE);
gROOT->SetStyle("Plain");
gStyle->SetOptStat(0);
TCanvas* theCanvas = new TCanvas("theCanvas", "", 1500, 600);
theCanvas->UseCurrentStyle();
theCanvas->Divide(2,1);
theCanvas->cd(1);
TH2D* plot1 = getDalitzPlot(fileName1, xInt1, yInt1, xLabel, yLabel);
plot1->Draw("colz");
latex.DrawLatex(0.3, 0.95, decay1.c_str());
theCanvas->cd(2);
TH2D* plot2 = getDalitzPlot(fileName2, xInt2, yInt2, xLabel, yLabel);
plot2->Draw("colz");
latex.DrawLatex(0.3, 0.95, decay2.c_str());
theCanvas->Print("PhaseSpacePlots.png");
}
TH2D* getDalitzPlot(string fileName, int xInt = 1, int yInt = 2,
string xLabel = "m^{2}_{23}",
string yLabel = "m^{2}_{13}") {
TFile* theFile = new TFile(fileName.c_str(), "read");
if (!theFile) {return 0;}
TTree* theTree = dynamic_cast<TTree*>(theFile->Get("dalitzTree"));
if (!theTree) {return 0;}
double m12Sq(0.0), m23Sq(0.0), m13Sq(0.0);
theTree->SetBranchAddress("invMass12Sq", &m12Sq);
theTree->SetBranchAddress("invMass23Sq", &m23Sq);
theTree->SetBranchAddress("invMass13Sq", &m13Sq);
// Turn off unneeded branches
theTree->SetBranchStatus("invMass12", 0);
theTree->SetBranchStatus("invMass23", 0);
theTree->SetBranchStatus("invMass13", 0);
double m12SqMin = theTree->GetMinimum("invMass12Sq");
double m12SqMax = theTree->GetMaximum("invMass12Sq");
double m23SqMin = theTree->GetMinimum("invMass23Sq");
double m23SqMax = theTree->GetMaximum("invMass23Sq");
double m13SqMin = theTree->GetMinimum("invMass13Sq");
double m13SqMax = theTree->GetMaximum("invMass13Sq");
// Default: xAxis = m23Sq, yAxis = m13Sq, for xInt = 1, yInt = 2
double xMin(m23SqMin), xMax(m23SqMax);
double yMin(m13SqMin), yMax(m13SqMax);
if (xInt == 2) {
xMin = m13SqMin; xMax = m13SqMax;
} else if (xInt == 3) {
xMin = m12SqMin; xMax = m12SqMax;
}
if (yInt == 1) {
yMin = m23SqMin; yMax = m23SqMax;
} else if (yInt == 3) {
yMin = m12SqMin; yMax = m12SqMin;
}
// Round off limits to nearest "integer"
xMin = roundMin(xMin);
xMax = roundMax(xMax);
yMin = roundMin(yMin);
yMax = roundMax(yMax);
TH2D* DPHist = new TH2D("DPHist", "", 100, xMin, xMax, 100, yMin, yMax);
DPHist->SetDirectory(0);
DPHist->SetXTitle(xLabel.c_str());
DPHist->SetYTitle(yLabel.c_str());
int N = theTree->GetEntries();
int i(0);
for (i = 0; i < N; i++) {
theTree->GetEntry(i);
if (i%100000 == 0) {cout<<"i = "<<N-i<<endl;}
double xVal = m23Sq;
double yVal = m13Sq;
if (xInt == 2) {
xVal = m13Sq;
} else if (xInt == 3) {
xVal = m12Sq;
}
if (yInt == 1) {
yVal = m23Sq;
} else if (yInt == 3) {
yVal = m12Sq;
}
DPHist->Fill(xVal, yVal);
}
return DPHist;
}
double roundMin(double x) {
double value = floor(x)*0.98;
return value;
}
double roundMax(double x) {
double value = ceil(x)*1.02;
return value;
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 4:40 AM (13 h, 26 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6457109
Default Alt Text
PhaseSpacePlots.C (3 KB)

Event Timeline