Divided histograms don't show in Canvas

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;
}

The last two canvas->cd(...); should be comCanvas->cd(...);.

Thanks, and honestly, I’m kind of embarrased, this was too simple. You are a good person Wile_E_Coyote

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.