//This can be run in ROOT by doing .x convolution.C[++]
#include <TMath.h>
#include <TStyle.h>
#include <TF1.h>
#include <TGraph.h>
#include <TLegend.h>
#include <TCanvas.h>
//Breit-Wigner Properties
Double_t bwMean = 0.0;
Double_t bwGamma = 2.5;
//Gaussian Properties
Double_t gWidth = 2.0;
//Set Precision of the Integrals (number of steps)
Int_t precision = 200;
//Create an animated gif?
Bool_t createGIF = kFALSE; // kFALSE or kTRUE
//Instantiate the functions and graphs to be used
//Their properties will be set later
TF1 *func1; //The Function to be reflected ("stationary")
TF1 *func2; //The function to be translated ("moving")
TF1 *func3; //The product of the two functions (func1 and func2)
TF1 *func4; //The overlap of the two functions (func1 and func2)
TGraph *graph4; //The Resulting Convolution Graph
//Does the product of the two functions (func1 and func2)
Double_t multipliedFuncs(Double_t *x, Double_t *par) {
#if !defined(__CINT__)
par = par; // get rid of compiler's warning: unused parameter ‘par’
#endif /* !defined(__CINT__) */
Double_t xx = x[0];
return func1->Eval(xx) * func2->Eval(xx);
}//End multipliedFuncs
//Find the smaller value of the functions (func1 and func2)
//This is used to shade the overlap region
Double_t smallerValue(Double_t *x, Double_t *par) {
#if !defined(__CINT__)
par = par; // get rid of compiler's warning: unused parameter ‘par’
#endif /* !defined(__CINT__) */
Double_t xx = x[0];
return TMath::Min(func1->Eval(xx), func2->Eval(xx));
}//End smallerValue
void convolution() {
//Define the Breit-Wigner (func1)
func1 = new TF1("BreitWigner","TMath::BreitWigner(x,[0],[1])",-10,10);
func1->SetTitle("Convolution of a Breit-Wigner and a Gaussian");
func1->SetParameters(bwMean,bwGamma);
func1->SetNpx(precision);
func1->SetLineColor(kBlue);
func1->SetLineWidth(3);
//Create the Gaussian (func2)
func2 = new TF1("Gaussian","TMath::Gaus(x,[0],[1],1)",-10,10);
func2->SetParameters(-10,gWidth);
func2->SetNpx(precision);
func2->SetLineColor(kBlack);
func2->SetLineWidth(3);
//Define the Product of the two functions (func1 and func2)
func3 = new TF1("MultipliedFuncs",multipliedFuncs,-10,10,0);
func3->SetNpx(precision);
func3->SetLineColor(kOrange);
func3->SetLineWidth(3);
func3->SetFillColor(kOrange);
func3->SetFillStyle(1001);
//Define the Overlap of the two functions (func1 and func2)
func4 = new TF1("OverlapFuncs",smallerValue,-10,10,0);
func4->SetNpx(precision);
func4->SetLineColor(kYellow);
func4->SetLineWidth(3);
func4->SetFillColor(kYellow);
func4->SetFillStyle(1001);
//Create Convolution Graph
graph4 = new TGraph();
graph4->SetName("Convolution");
graph4->SetLineColor(kRed);
graph4->SetLineWidth(3);
//Just a trick ... we will draw a line, so we need at least two points ...
// graph4->SetPoint((graph4->GetN()), -10, (0.9 * func3->Integral(-10, 10)));
graph4->SetPoint((graph4->GetN()), -10, func3->Integral(-10, 10));
//Create a Legend
TLegend *legend = new TLegend(0.55, 0.72, 0.88, 0.88, "", "brNDC");
legend->SetFillColor(gStyle->GetCanvasColor());
legend->AddEntry(func1,
TString::Format("Breit-Wigner ( #mu = %.1f, #Gamma = %.1f )", bwMean, bwGamma),
"l");
legend->AddEntry(func2,
TString::Format("Gaussian ( #sigma = %.1f )", gWidth),
"l");
legend->AddEntry(graph4, "Convolution of BW and Gaussian", "l");
legend->AddEntry(func4, "Overlap Area of BW and Gaussian", "f");
legend->AddEntry(func3, "Product of BW and Gaussian", "f");
//Create a Canvas to Draw to
TCanvas *canvas = new TCanvas("canvas","canvas");
//Draw the Functions and Graphs, Legend
func1->Draw("L");
func4->Draw("L SAME");
func3->Draw("L SAME");
func1->Draw("L SAME");
func2->Draw("L SAME");
graph4->Draw("L");
legend->Draw();
gPad->RedrawAxis(); // Redraw axes
// gPad->RedrawAxis("g"); // Redraw gridx and/or gridy
gPad->Modified(); gPad->Update();
//Create an animated gif
if (createGIF){
canvas->SaveAs("BW_Gaus_Convolution.gif");
}
//Amount by which to move the second function every step
Double_t moveOver = 0.05;
//Loop to move the second function toward the first
for (Int_t i = 0; i < 400; i++) {
//Set the Mean of the Guassian - This moves the function
Double_t xx = func2->GetParameter(0) + moveOver;
func2->SetParameter(0, xx);
graph4->SetPoint((graph4->GetN()), xx, func3->Integral(-10, 10));
//Update the picture
gPad->Modified(); gPad->Update();
//Update an animated gif
if (createGIF){
canvas->SaveAs("BW_Gaus_Convolution.gif+50");
}
}
}
BTW. In your original post, you were explicitly asking to draw the “product” of both functions as the “yellow area”, not the “overlap” (and so, the “yellow area” was correct, but it was simply the “product”, not the “overlap”).