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