Transforming bivariate Gaussians instead of generating

ROOT
10

This should do the job. You just need to take care of normalizing correctly the area / amplitude (integral of counts per beamlet).

TF2* f2 = new TF2("f2","bigaus");
Double_t meanX = -3;
Double_t sigmaX = 0.2;
Double_t meanY = 5;
Double_t sigmaY = 0.5;
Double_t correlation = 0.01;
Double_t amplitude = 1; // Change this with Area/sqrt(2pi sigmax sigmay) to normalize to Area. You can also divide by binWidthX*binWidthY if you want a density
f2->SetParameters(amplitude,meanX,sigmaX,meanY,sigmaY,correlation);// “const”, “meanx”, “sigmax”, “meany”, “sigmay”, “correlation”
Double_t maxSigmaFactor = 3; // stop after +/-3*sigma
TH2F* h2 = new TH2F("h2", "h2", 500, -10, 10, 500, -5, 8);
Int_t startX = h2->GetXaxis()->FindBin(meanX - maxSigmaFactor * sigmaX);
Int_t endX = h2->GetXaxis()->FindBin(meanX + maxSigmaFactor * sigmaX);
Int_t startY = h2->GetYaxis()->FindBin(meanY - maxSigmaFactor * sigmaY);
Int_t endY = h2->GetYaxis()->FindBin(meanY + maxSigmaFactor * sigmaY);

for(Int_t i = std::max(1,startX); i<=std::min(h2->GetNbinsX(),endX); ++i)
    for(Int_t j = std::max(1,startY); j<=std::min(h2->GetNbinsY(),endY); ++j)
        h2->SetBinContent(i,j,f2->Eval(h2->GetXaxis()->GetBinCenter(i), h2->GetYaxis()->GetBinCenter(j)));
auto c = new TCanvas();
h2->Draw("COLZ");

image

1 Like