#include <Riostream.h>  // for cout, iostream
#include <TFile.h>
#include <TStyle.h>
#include <TH2.h>
#include <TTree.h>
#include <TH1.h>
#include <TF1.h>
#include <TGaxis.h>
#include <TPad.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TLegendEntry.h>
#include<TColor.h>
#include<TMath.h>

constexpr Double_t Pi()
{
   return 3.14159265358979323846;
}

void Re_Smear_3gauss()
{


gStyle->SetOptStat(1002211);
gStyle->SetOptFit(1111);

/////////////////// Tree opening \\\\\\\\\\\\\\\

TFile *file1=TFile::Open("Analysis.root");
TTree *tree1=(TTree*)file1->Get("Analysis");

TH1D*h1 = new TH1D("h1","h1", 101.0,39.5,140.5);

//Defining a TF1Convolution
    TF1Convolution *f_conv = new TF1Convolution("[0]*exp(-((x-[1])/[2])**2) + [3]*exp(-((x-([1]-[4]))/[5])**2) + [3]*exp(-((x-([1]+[4]))/[5])**2)", "gaus", 46.,130., true);
    
    
f_conv->SetRange(30.,150.);
f_conv->SetNofPointsFFT(1000);


    TF1 *f        = new TF1("f", *f_conv, 40., 140., f_conv->GetNpar());
    TF1 *f_no_res = new TF1("f_no_res", "[0]*exp(-((x-[1])/[2])**2) + [3]*exp(-((x-([1]-[4]))/[5])**2) + [3]*exp(-((x-([1]+[4]))/[5])**2)", 46.,130.);
    
    f->SetParameters(1e3,88.,13., 1e3 ,30.,30., 1.,0.,3.5);
    
    f->SetParNames("A_{1}","#mu_{1}","#sigma_{1}","A_{2}","#mu_shift","#sigma_{2}","Amplitude Convolution","#mu_conv","#sigma_res");

    // I set the amplitude of the 2 gaussian to be positive
    f->SetParLimits(0,0.,1e6);
    f->SetParLimits(3,0.,1e6);
    
    // I set the mean  of the 2 gaussian to be between 60 and 120
    f->SetParLimits(1,60.,120.);
    f->SetParLimits(4,0.,20.);
    
    // I set the sigma  of the 2 gaussian to be between 0 and 100
    f->SetParLimits(2,1.,40.);
    f->SetParLimits(5,1.,40.);
    
// i fix the parameters of the convoluted gaussian
    f->FixParameter(6,1./(3.5*TMath::Sqrt(2*TMath::Pi())) );
    f->FixParameter(7,0.);
    f->FixParameter(8,3.5);
    
    //  f_no_res->SetParameters(f->GetParameters());

 TCanvas *c1 = new TCanvas("c1","Without constraint",800,600);
c1->Divide(2,2);
c1->cd(1);
tree1->Draw("M_Tot>>h1");
h1->Draw("E1");

c1->cd(2);
TH1D *h1_clone = (TH1D*)h1->Clone();
h1_clone->Draw("E1");
h1_clone->Fit(f,"ME","",46.,130.);
    
f_no_res->SetParameters(f->GetParameters());
    
auto h2 = new TH1D("h2","h2_Mass_resolution",101.0,39.5,140.5);
h2 = (TH1D*)f_no_res->GetHistogram();
h2->SetStats(kTRUE);
h2->SetLineColorAlpha(kBlue, 0.35);

TF1 *f12 = new TF1("f12","[0]*exp(-0.5*((x-[1])/[2])**2) + [3]*exp(-0.5*((x-([1]-[4]))/[5])**2) + [3]*exp(-0.5*((x-([1]+[4]))/[5])**2)",46.,130.);


   f12->SetFillColor(19);
   f12->SetFillStyle(0);
   f12->SetLineColor(2);
   f12->SetLineWidth(2);
  

   f12->SetParameter(0,1000);
   f12->SetParLimits(0,0,100000);

   f12->SetParameter(1, 88);
   f12->SetParLimits(1,0,150);

   f12->SetParameter(2, 12);
   f12->SetParLimits(2,0,20);

   f12->SetParameter(3, 3000);
   f12->SetParLimits(3,0,80000);

   f12->SetParameter(4, 13.8);
   f12->SetParLimits(4,0,50);

   f12->SetParameter(5, 8);
   f12->SetParLimits(5,2,40);


c1->cd(3);
h2->Fit("f12","ME", "",46.,130.);
h2->Draw("E1");   

Double_t * p1 = f12->GetParameters();
cout<<"Mode Composition"<<"\n";


   TF1 *f11 = new TF1("f11","[0]*exp(-0.5*((x-[1])/[2])**2)",46.,130.);
   f11->SetFillColor(19);
   f11->SetFillStyle(0);

   f11->SetLineColor(kBlack);
   f11->SetLineStyle(1);
   f11->SetLineWidth(2);

   f11->GetXaxis()->SetLabelFont(42);
   f11->GetXaxis()->SetTitleOffset(1);
   f11->GetXaxis()->SetTitleFont(42);
   f11->GetYaxis()->SetLabelFont(42);
   f11->GetYaxis()->SetTitleFont(42);

   f11->SetParameter(0,p1[0]);
   f11->SetParameter(1, p1[1]);
   f11->SetParameter(2,p1[2]);
   f11->Draw("same");

TF1 *f122 = new TF1("f122","[3]*exp(-0.5*((x-(88-[4]))/[5])**2)",46.,130.);
   f122->SetFillColor(19);
   f122->SetFillStyle(0);
   f122->SetLineColor(kGreen+1);
   f122->SetLineStyle(8);
   f122->SetLineWidth(5);

   f122->GetXaxis()->SetLabelFont(42);
   f122->GetXaxis()->SetTitleOffset(1);
   f122->GetXaxis()->SetTitleFont(42);
   f122->GetYaxis()->SetLabelFont(42);
   f122->GetYaxis()->SetTitleFont(42);


   f122->SetParameter(3,p1[3]);
   f122->SetParameter(4,p1[4]);
   f122->SetParameter(5,p1[5]);
   f122->Draw("same");

   TF1 *f13 = new TF1("f13","[3]*exp(-0.5*((x-(88+[4]))/[5])**2)",46.,130.);
   f13->SetFillColor(19);
   f13->SetFillStyle(0);
   f13->SetLineColor(kGreen+1);
   f13->SetLineStyle(8);
   f13->SetLineWidth(5);

   f13->GetXaxis()->SetLabelFont(42);
   f13->GetXaxis()->SetTitleOffset(1);
   f13->GetXaxis()->SetTitleFont(42);
   f13->GetYaxis()->SetLabelFont(42);
   f13->GetYaxis()->SetTitleFont(42);


   f13->SetParameter(3,p1[3]);
   f13->SetParameter(4,p1[4]);
   f13->SetParameter(5,p1[5]);
   f13->Draw("same");


Double_t Sym_S12, Asym_A12, Asym_A22, Tot_Area2;
Double_t Mode_S12, Mode_A12, Mode_A22;

// Area of modes 

 Sym_S12  = p1[0]*TMath::Sqrt(2*Pi()*pow(p1[2],2));
 Asym_A12 = 2*(p1[3]*TMath::Sqrt(2*Pi()*pow(p1[5],2)));
 Tot_Area2 = Sym_S12 + Asym_A12;

// Calculating mode composition

Mode_S12 = (Sym_S12/Tot_Area2)*100;
Mode_A12 = (Asym_A12/Tot_Area2)*100;


cout<<"Mode_S1 : "<<Mode_S12<<" %"<< "\t" << "Symmetric peak position   : "<<p1[1]<<"\n";
cout<<"Mode_A1 : "<<Mode_A12<<" %"<< "\t" << "Displacement Asym A1 mode : "<<p1[4]<<"\n";


}