I’ve managed to reproduce the results in this paper, I got their data from this website, and now the last thing I have to do is to dive the histograms to see how close my data is to theirs.
The thing is that when I divide them something strange happens:
Even though I get no warning nor error from the Canvas only the Transverse momentum histogram shows up, in the wrong pad mind you, should be in the second but appears in the first.
I don’t know what’s going on, I am this close to finishing. Here is my code, I’ve signaled with a “@@@@@@@@@@” the part where I get the bin edges and with a “#########################” the part where try to do the ratios
#include <iostream>
using namespace std;
#include "Pythia8/Pythia.h"
#include "TROOT.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TF1.h"
#include "TFile.h"
using namespace Pythia8;
int main()
{
//get pro data
TFile* proData1 = TFile::Open("/Users/Fer/Documents/Pythia/myhists/HEPData-ins1395253-v1-Table1.root");
TFile* proData2 = TFile::Open("/Users/Fer/Documents/Pythia/myhists/HEPData-ins1395253-v1-Table2.root");
//get pro histograms
TH1F* proHist_eta_inel0= (TH1F*)proData1->Get("Table 1/Hist1D_y1");
TH1F* proHist_eta_inel= (TH1F*)proData1->Get("Table 1/Hist1D_y2");
//don't confuse, variables with same name from different files
TH1F* proHist_ptz_inel0= (TH1F*)proData2->Get("Table 2/Hist1D_y1");
//Save pro histograms
TCanvas* procanvas = new TCanvas("ProCanvas", "Eta y pTz", 600, 600);
procanvas->Divide(2,2);
procanvas->cd(1);
proHist_eta_inel->Draw();
procanvas->Update();
procanvas->cd(2);
proHist_ptz_inel0->Draw();
procanvas->Update();
procanvas->cd(3);
proHist_eta_inel0->Draw();
procanvas->Update();
procanvas->SaveAs("prohistos.pdf");
procanvas->SaveAs("prohistos.root");
procanvas->Close();
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//Get the data from the pro histograms
//Get bin number
int proBinNumber_eta_inel0 = proHist_eta_inel0->GetNbinsX();
int proBinNumber_eta_inel = proHist_eta_inel->GetNbinsX();
int proBinNumber_ptz_inel0 = proHist_ptz_inel0->GetNbinsX();
//this are used later
double lastEdge = 0;
//For to get the edges for Eta at INEL<0
double proEdges_eta_inel0[proBinNumber_eta_inel0+1];
for (int bin = 1; bin < proBinNumber_eta_inel0+1; ++bin)
{
proEdges_eta_inel0[bin-1]=proHist_eta_inel0->GetXaxis()->GetBinLowEdge(bin);
}
//Get the last edge outside the for
lastEdge= proHist_eta_inel0->GetXaxis()->GetBinUpEdge(proBinNumber_eta_inel0);
proEdges_eta_inel0[proBinNumber_eta_inel0]=lastEdge;
//Eta inel edges
double proEdges_eta_inel[proBinNumber_eta_inel+1];
for (int bin = 1; bin < proBinNumber_eta_inel+1; ++bin)
{
proEdges_eta_inel[bin-1]=proHist_eta_inel->GetXaxis()->GetBinLowEdge(bin);
}
lastEdge= proHist_eta_inel->GetXaxis()->GetBinUpEdge(proBinNumber_eta_inel);
proEdges_eta_inel[proBinNumber_eta_inel]=lastEdge;
//Ptz inel0 edges
double proEdges_ptz_inel0[proBinNumber_ptz_inel0+1];
for (int bin = 1; bin < proBinNumber_ptz_inel0+1; ++bin)
{
proEdges_ptz_inel0[bin-1]=proHist_ptz_inel0->GetXaxis()->GetBinLowEdge(bin);
}
lastEdge= proHist_ptz_inel0->GetXaxis()->GetBinUpEdge(proBinNumber_ptz_inel0);
proEdges_ptz_inel0[proBinNumber_ptz_inel0]=lastEdge;
//Create my histograms with the number of bins and edges of the pro histograms
TH1F *hist_eta_inel0 = new TH1F("Pseudorapidez", "Pseudorapidez", proBinNumber_eta_inel0, proEdges_eta_inel0);
TH1F *hist_eta_inel = new TH1F("Pseudorapidez", "Pseudorapidez", proBinNumber_eta_inel, proEdges_eta_inel);
TH1F *hist_ptz_inel0 = new TH1F("Momento Transverso", "Momento Transverso", proBinNumber_ptz_inel0, proEdges_ptz_inel0);
TCanvas *canvas = new TCanvas("canvas", "Eta y pTz", 600, 600);
canvas->Divide(2,2);
Pythia pythia;
//Beam propoerties
int Nev=10000;
pythia.readString("SoftQCD:inelastic = on");
pythia.readString("Beams:eCM = 13000.");
pythia.readString("Tune:pp = 14");
pythia.init();
//events are created, selected, and the histograms filled
for (int iEvent = 0; iEvent < Nev; ++iEvent)
{
if (!pythia.next()) continue; // sirve para evtar eventos vacios
for (int i = 0; i < pythia.event.size(); ++i) //anaizamos particulas en el eventos
{
//INEL eta
if (pythia.event[i].isFinal() && pythia.event[i].isCharged() && abs(pythia.event[i].eta())<=1.8 )
{
hist_eta_inel->Fill(pythia.event[i].eta());
}
//INEL>0 eta and pTz
if (pythia.event[i].isFinal() && pythia.event[i].isCharged() && pythia.event[i].pT()>=0 && abs(pythia.event[i].eta())<=0.8 )
{
hist_eta_inel0->Fill(pythia.event[i].eta());
hist_ptz_inel0->Fill(pythia.event[i].pT());
}
}
}
//create new histograms to be filled with the normalized data
TH1F *hist_eta_inel02 = new TH1F("Pseudorapidez INEL<0", "PseudorapidezINEL<0", proBinNumber_eta_inel0, proEdges_eta_inel0);
TH1F *hist_eta_inel2 = new TH1F("Pseudorapidez INEL", "Pseudorapidez INEL", proBinNumber_eta_inel, proEdges_eta_inel);
TH1F *hist_ptz_inel02 = new TH1F("Momento Transverso INEL<0", "Momento Transverso INEL<0", proBinNumber_ptz_inel0, proEdges_ptz_inel0);
double c1=0;
double c2=0;
double c3=0;
double dN=0;
double dptz=0;
double deta=0;
//Normalize each bin by the proper constants c1 c2 and c3
//INEL>0 eta
for (int bin = 1; bin <= proBinNumber_eta_inel0; ++bin)
{
deta=hist_eta_inel0->GetBinWidth(bin);
dN= hist_eta_inel0->GetBinContent(bin);
c1=dN/(deta*Nev);
hist_eta_inel02->SetBinContent(bin,c1);
}
//INEL eta
for (int bin = 1; bin <= proBinNumber_eta_inel; ++bin)
{
deta=hist_eta_inel->GetBinWidth(bin);
dN= hist_eta_inel->GetBinContent(bin);
c2=dN/(deta*Nev);
hist_eta_inel2->SetBinContent(bin,c2);
}
//INEL>0 pTz
deta=1.6;
double ptz=0;
for (int bin = 1; bin <= proBinNumber_ptz_inel0; ++bin)
{
ptz= hist_ptz_inel0->GetBinCenter(bin);
dptz= hist_ptz_inel0->GetBinWidth(bin);
dN= hist_ptz_inel0->GetBinContent(bin);
c3=(dN/(Nev*deta*dptz))*(1/(2*3.1459*ptz));
hist_ptz_inel02->SetBinContent(bin,c3);
}
pythia.stat();
canvas->cd(1);
hist_eta_inel2->Draw();
canvas->Update();
canvas->cd(2);
hist_ptz_inel02->Draw();
canvas->Update();
canvas->cd(3);
hist_eta_inel02->Draw();
canvas->Update();
canvas->SaveAs("histos.pdf");
canvas->SaveAs("histos.root");
canvas->Close();
//###################################################################################
//Comparison
//divide each histogram by its pro counterpart
hist_eta_inel02->Divide(proHist_eta_inel0,hist_eta_inel02);
hist_eta_inel2->Divide(proHist_eta_inel,hist_eta_inel2);
hist_ptz_inel02->Divide(proHist_ptz_inel0,hist_ptz_inel02);
TCanvas* comCanvas = new TCanvas("ComCanvas", "Comparaciones", 600, 600);
comCanvas->Divide(2,2);
comCanvas->cd(1);
hist_eta_inel2->Draw();
comCanvas->Update();
canvas->cd(2);
hist_eta_inel02->Draw();
comCanvas->Update();
canvas->cd(3);
hist_ptz_inel02->Draw();
comCanvas->Update();
comCanvas->SaveAs("comHistos.pdf");
comCanvas->SaveAs("comHistos.root");
comCanvas->Close();
return 0;
}