Contour lines and colz option in TH2D


Please read tips for efficient and successful posting and posting code

ROOT Version: 6.16/00
Platform: Ubuntu 20.04
Compiler: linuxx8664gcc


Hello co-rooters!

I am trying to plot a TH2D histogram using the colz option.
I also want the contour lines for z=1, 2, 3, … and later on I’d like to add labels (i.e. z=1, z=2, etc).

I firstly draw using `colz’ and then, after I set the contours, I redraw it using 'cont3 same". Here are the differences when I only draw with ‘colz’

Screenshot from 2024-12-23 14-42-31

and with colz and then with cont3 same

Screenshot from 2024-12-23 14-43-37

You see that around the red mark, there is white gap which is unclear where it comes from.
Any idea on why the white gap is occurring and how can I overlay the colz options with only the contour lines? Thanks in advance!

My code is (it won’t run because I couldn’t upload the dat file) :

#include "TF1.h"
#include "TF2.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TPad.h"
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TGraph.h"
#include "TText.h"
#include "TGraphQQ.h"
#include "TCanvas.h"
#include "TFrame.h"
#include "TLegend.h"
#include "TLine.h"
#include "Math/PdfFuncMathCore.h"
#include <TFitResult.h>
#include "TGraphErrors.h"
#include <TMultiGraph.h>
#include "Math/Functor.h"
#include "Math/WrappedFunction.h"
#include "Math/IFunction.h"
#include "Math/Integrator.h"

#include <algorithm>


#include <vector>
#include <iostream>
#include <fstream>

using namespace std;

// Define the threshold function
double Threshold(double *x, double *par) {
    double E = x[0];
    double a = par[0];
    double b = par[1];
    return a * pow(E, b);
}

// Define the normalized chi-square PDF
double h0(double *x, double *par) {
    double Gn0     = x[0];
    double C       = par[0];
    double avg     = par[1];
    double nu      = par[2];
    double chi_pdf = C*ROOT::Math::chisquared_pdf(Gn0 / avg, nu);
    return chi_pdf;
}

void MLA() {
    // Step 1: Load data
    std::ifstream file("En_Gn0.dat");
    if (!file.is_open()) {
        std::cerr << "Error: Unable to open data file." << std::endl;
    }

    std::vector<double> energies, widths;
    double energy, width;
    while (file >> energy >> width) {
        energies.push_back(energy);
        widths.push_back(width);
    }
    file.close();

    // Step 2: Plot data with threshold    
    TGraph *g = new TGraph(energies.size(), energies.data(), widths.data());
    g->SetMarkerStyle(20);
    g->SetTitle("Energy vs Neutron Widths;Energy (eV);Gn0 (meV)");
    //g->Draw("AP");

	double a = 1.8e-2;
    double b = 0.5;

	TF1 *threshold = new TF1("threshold", Threshold, *std::min_element(energies.begin(), energies.end()), *std::max_element(energies.begin(), energies.end()), 2);
    threshold->SetParameters(a, b);
    threshold->SetLineColor(kRed);
    threshold->SetTitle("Threshold");
    //threshold->Draw("LSAME");
    
     // Step 3: Filter data above threshold
    std::vector<double> filtered_widths;
    std::vector<double> filtered_energies;
    for (size_t i = 0; i < energies.size(); ++i) {
        if (widths[i] > Threshold(&energies[i], threshold->GetParameters())) {
            filtered_widths.push_back(widths[i]);
            filtered_energies.push_back(energies[i]);
        }
    }

    // Step 4: Define grid for nu and avg
    int    nu_points = 300 , avg_points = 500;
    double nu_min    = 0.1, nu_max     = 3.1;
    double avg_min   = 0.1, avg_max    = 5.1;

    TH2D *h_likelihood      = new TH2D("h_likelihood"     , "Likelihood L;Average Gn0;Degrees of Freedom (nu)", avg_points, avg_min, avg_max, nu_points, nu_min, nu_max);
    TH2D *h_likelihood_Paul = new TH2D("h_likelihood_Paul", "Likelihood L;Average Gn0;Degrees of Freedom (nu)", avg_points, avg_min, avg_max, nu_points, nu_min, nu_max);

	TF1 *f_h0 = new TF1("f_h0", h0, *std::min_element(filtered_widths.begin(), filtered_widths.end()), 3.e3, 3);

	double bin_content = 1.;
	double binCenterX  = 1.;
	double binCenterY  = 1.;
	double C           = 1.;
	double I           = 1.;
	double Lmax        = 0.;
	
	for (int i=0; i<avg_points; ++i){
		for (int j=0; j<nu_points; ++j){
			binCenterX = h_likelihood->GetXaxis()->GetBinCenter(i+1);
			binCenterY = h_likelihood->GetYaxis()->GetBinCenter(j+1);
			f_h0->SetParameter(0, 1.);
			f_h0->SetParameter(1, binCenterX);
			f_h0->SetParameter(2, binCenterY);
			//f_h0->SetTitle(TString::Format("f_h0_%f_%f", binCenterX, binCenterY));
			
			
			C = f_h0->Integral(*std::min_element(filtered_widths.begin(), filtered_widths.end()), 3.e3 , 1.e-10);
			f_h0->SetParameter(0, 1./C);
			
			I = f_h0->Integral(*std::min_element(filtered_widths.begin(), filtered_widths.end()), 3.e3 , 1.e-10);
			
			bin_content = 1.;
			//bin_content = TMath::Log( f_h0->Eval(filtered_widths[0]) );
			
			for (int k=0; k<filtered_energies.size(); ++k){
				bin_content *= f_h0->Eval( filtered_widths[k] );
			}
			h_likelihood->SetBinContent(i+1, j+1, bin_content);
		}
	}
	
	Lmax = h_likelihood->GetMaximum();
	// Find the bin indices of the maximum bin
	int binX, binY, binZ;
	h_likelihood->GetMaximumBin(binX, binY, binZ);

	// Get the bin centers in X and Y
	double avg_Gn0 = h_likelihood->GetXaxis()->GetBinCenter(binX);
	double nu      = h_likelihood->GetYaxis()->GetBinCenter(binY);
	
	cout << "Lmax = " << Lmax  << ", Avg_Gn0 = " << avg_Gn0 << ", nu = " << nu << endl;
	
	for (int i=0; i<avg_points; ++i){
		for (int j=0; j<nu_points; ++j){
			binCenterX = h_likelihood->GetXaxis()->GetBinCenter(i+1);
			binCenterY = h_likelihood->GetYaxis()->GetBinCenter(j+1);
			
			bin_content = TMath::Sqrt(2.)*TMath::Sqrt(TMath::Log(Lmax) - TMath::Log( h_likelihood->GetBinContent(i+1, j+1) ) );
			h_likelihood_Paul->SetBinContent(i+1, j+1, bin_content);
		}
	}
	
	h_likelihood_Paul->Draw("colz");
	
	// Find the maximum bin content
    double zMax = h_likelihood_Paul->GetMaximum();

    //Define sigma levels for contour lines up to zMax
    const int nLevels = zMax;  // Number of levels
    double sigmaLevels[nLevels];
    for (int i = 0; i < nLevels; ++i) {
        sigmaLevels[i] = zMax * (i + 1.0) / nLevels;
    }

    //Set contour levels and draw the contours
    h_likelihood_Paul->SetContour(nLevels, sigmaLevels);
    h_likelihood_Paul->Draw("CONT3 same"); 
	
	TGraph *g_Lmax = new TGraph(1);
	g_Lmax->SetPoint(0, avg_Gn0, nu);
	g_Lmax->SetMarkerStyle(20);
	g_Lmax->SetMarkerSize(2.);
	g_Lmax->SetMarkerColor(kRed);
	g_Lmax->Draw("Psame");
	
	TF1 *f_nu1 = new TF1("f_nu1", "1", avg_min, avg_max);
	f_nu1->SetLineColor(kBlack);
	f_nu1->SetLineWidth(2);
	f_nu1->SetLineStyle(9);
	f_nu1->Draw("Lsame");

}

Looking at the first plot, the data inside the white area fall in the lowest colour bin, which is from 0 to 1.xx (as there are around 3.5 bins from 0 to 5).
In the second plot you can see that the lowest values (~ 0-1) have the colour white. Since you are redefining the levels of the same histogram to do the second plot (with nLevels and sigmaLevels), it looks like you are leaving out the lowest values. So, you can try adjusting the levels to ensure you include the lowest values in the data (maybe from zero?) or, alternatively, maybe use DrawCopy the first time, before modifying the levels, and then Draw for the second:

//...
	h_likelihood_Paul->DrawCopy("colz");	
	// Find the maximum bin content, etc.
//...
    h_likelihood_Paul->Draw("CONT3 same");

Thank you for your help!
I didn’t think of using DrawCopy which seems to do the job just fine!!