Generate TH2 from TF1 with a scatter around its trend

[vlisovsk]

Hello,

I am trying to generate a TH2 histogram from a given TF1 function that provides the trend y:x (TrendFunction), and another TF1 function that provides the y resolution (i.e. the scatter of points around this trend; ResolutionFunction). This plot gives an idea of what I mean:

Right now, my implementation is the following (the actual code is below):
(1) loop over every X bin;
(2) evaluate the TrendFunction at each X bin center;
(3) generate a 1D histogram from the ResolutionFunction with its mean set to the value obtained from the previous step (2);
(4) fill each X slice of my 2D histogram with the contents of the 1D histograms obtained in the step (3) by looping over its Y bins

This works ok (and produces the plot shown above), but I see problems in case the TrendFunction is rapidly rising and the binning is coarse, when the steps 2-3 of my procedure produce biased results. This is because I only rely on the bin center in my second step, while there can be huge variations of the TrendFunction value within one bin.

Therefore I am wondering if there is a better (simpler?) way to handle this generation problem?
Thanks.

My current code:

void test2Dgen(){
gStyle->SetOptStat(0);
int m_numXBins = 100;
int m_numYBins = 100;
int m_toGenerate = 1000;
TH2D* generated_histogram = new TH2D("generated_histogram", ";x;y", m_numXBins, 0, 1, m_numYBins, 0, 1);

TF1* TrendFunction = new TF1("TrendFunction", "[0] + [1] * x", 0, 1);
TrendFunction->SetParameters(0.25, 0.25);

TF1 *ResolutionFunction = new TF1("ResolutionFunction", "[0]*TMath::Gaus(x, [1], [2])", 0, 1);
ResolutionFunction->SetParameters(1, 1, 0.05);

for (int xbin = 1; xbin <= m_numXBins + 1; xbin++) {

    double mean_value = TrendFunction->Eval(generated_histogram->GetXaxis()->GetBinCenter(xbin));
    ResolutionFunction->FixParameter(1, mean_value);

    // create a projection (1D histogram) in a given x bin
    TH1D* x_slice = (TH1D*)generated_histogram->ProjectionY("slice", xbin, xbin);
    x_slice->FillRandom(ResolutionFunction->GetName(), m_toGenerate);

    // normalise each x slice to unity
    if (x_slice->Integral() > 0) {
      x_slice->Scale(1. / x_slice->Integral());
    }
    // fill back the 2D histo with the result
    for (int ybin = 0; ybin <= m_numYBins + 1; ybin++) {
      generated_histogram->SetBinContent(xbin, ybin, x_slice->GetBinContent(ybin));
      generated_histogram->SetBinError(xbin, ybin, x_slice->GetBinError(ybin));
    }
  }
generated_histogram->Draw("COLZ");
TrendFunction->Draw("SAME");
}

ROOT Version: 6.26/14
Platform: linuxx8664gcc
Compiler: g++ (GCC) 14.1.0

---

[Danilo]

Hi,

I don’t think ROOT provides anything to automate this kind of operation.

Best,
D

[Dilicus]

Hi @vlisovsk,
below Iyou found your code with my changes. Instead of generating a random gaussian for each slice, I generated uniform random number between x_min and x_max. Then I evaluated the mean value and generated a random gaussian around the mean value, and then I filled the 2D histogram. I renormalized the histogram to the numbers of x bins like you did in your code.

Attached you find a picutre of the histogram I generated with the code below.

 void test2Dgen(){
gStyle->SetOptStat(0);
TRandom3* rand=new TRandom3();
int m_numXBins = 100;
int m_numYBins = 100;
int m_toGenerate = 1000;
TH2D* generated_histogram = new TH2D("generated_histogram", ";x;y", m_numXBins, 0, 1, m_numYBins, 0, 1);

TF1* TrendFunction = new TF1("TrendFunction", "[0] + [1] * x", 0, 1);
TrendFunction->SetParameters(-1, 3);

TF1 *ResolutionFunction = new TF1("ResolutionFunction", "[0]*TMath::Gaus(x, [1], [2])", 0, 1);
ResolutionFunction->SetParameters(1, 1, 0.05);
double x,mean_value,y;
for (int i=0;i<m_toGenerate*m_numXBins;i++) {

  x = rand->Uniform(0,1);
  mean_value = TrendFunction->Eval(x);
  ResolutionFunction->FixParameter(1, mean_value);

  y = rand->Gaus(mean_value,0.05);

  generated_histogram->Fill(x,y);

  }
 
generated_histogram->Scale(100./generated_histogram->Integral());
generated_histogram->Draw("COLZ");
TrendFunction->Draw("SAME");
}