#include <iostream>
#include "TCanvas.h"
#include "TH1.h"

void macro()
{
    //   TFile* f = new("file.root","RECREATE");
    TCanvas *c1 = new TCanvas("c1", "c1",16,30,700,500);

    Float_t xAxis[41] = {-5, -4, -3.5, -3, -2.5, -2, -1.8, -1.6, -1.4, -1.2, -1, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.5, 3, 3.5, 4., 5.};

    TH1F* histPtr = new TH1F("histPtr","90B_MRS #pi^{-},#pi^{+}",40,xAxis);
    histPtr->SetBinContent(5,0.03429431);
    histPtr->SetBinContent(6,0.09132375);
    histPtr->SetBinContent(7,0.1702631);
    histPtr->SetBinContent(8,0.3497895);
    histPtr->SetBinContent(9,0.7220107);
    histPtr->SetBinContent(10,1.50242);
    histPtr->SetBinContent(11,2.822891);
    histPtr->SetBinContent(12,4.279733);
    histPtr->SetBinContent(13,6.717525);
    histPtr->SetBinContent(14,10.26024);
    histPtr->SetBinContent(15,16.24549);
    histPtr->SetBinContent(16,26.02419);
    histPtr->SetBinContent(25,28.30753);
    histPtr->SetBinContent(26,17.15854);
    histPtr->SetBinContent(27,10.76037);
    histPtr->SetBinContent(28,6.739475);
    histPtr->SetBinContent(29,4.234743);
    histPtr->SetBinContent(30,2.826327);
    histPtr->SetBinContent(31,1.513767);
    histPtr->SetBinContent(32,0.6685444);
    histPtr->SetBinContent(33,0.3211848);
    histPtr->SetBinContent(34,0.1572088);
    histPtr->SetBinContent(35,0.08562638);
    histPtr->SetBinContent(36,0.02806619);
    histPtr->SetBinError(5,0.001139351);
    histPtr->SetBinError(6,0.003293322);
    histPtr->SetBinError(7,0.004831763);
    histPtr->SetBinError(8,0.007429496);
    histPtr->SetBinError(9,0.0116407);
    histPtr->SetBinError(10,0.01890519);
    histPtr->SetBinError(11,0.04140646);
    histPtr->SetBinError(12,0.0543397);
    histPtr->SetBinError(13,0.07988723);
    histPtr->SetBinError(14,0.1125551);
    histPtr->SetBinError(15,0.1766601);
    histPtr->SetBinError(16,0.3250688);
    histPtr->SetBinError(25,0.36238);
    histPtr->SetBinError(26,0.1904849);
    histPtr->SetBinError(27,0.1236194);
    histPtr->SetBinError(28,0.08403891);
    histPtr->SetBinError(29,0.0583256);
    histPtr->SetBinError(30,0.04308339);
    histPtr->SetBinError(31,0.02006726);
    histPtr->SetBinError(32,0.01166102);
    histPtr->SetBinError(33,0.007285574);
    histPtr->SetBinError(34,0.004740945);
    histPtr->SetBinError(35,0.003277002);
    histPtr->SetBinError(36,0.001057489);
    histPtr->SetEntries(92);
    histPtr->SetLineColor(kGreen);
    histPtr->SetMarkerColor(kGreen);
    histPtr->Draw("");


    TH1F *histPtr2 = new TH1F("histPtr2","90A_MRS #pi^{-},#pi^{+}",40, xAxis);
    histPtr2->SetBinContent(5,0.02775054);
    histPtr2->SetBinContent(6,0.08447972);
    histPtr2->SetBinContent(7,0.1681275);
    histPtr2->SetBinContent(8,0.3388171);
    histPtr2->SetBinContent(9,0.7020755);
    histPtr2->SetBinContent(10,1.51282);
    histPtr2->SetBinContent(11,2.799845);
    histPtr2->SetBinContent(12,4.267647);
    histPtr2->SetBinContent(13,6.618573);
    histPtr2->SetBinContent(14,10.28068);
    histPtr2->SetBinContent(15,17.05969);
    histPtr2->SetBinContent(16,27.1398);
    histPtr2->SetBinContent(26,17.3341);
    histPtr2->SetBinContent(27,10.64911);
    histPtr2->SetBinContent(28,6.701062);
    histPtr2->SetBinContent(29,4.437294);
    histPtr2->SetBinContent(30,2.912697);
    histPtr2->SetBinContent(31,1.561351);
    histPtr2->SetBinContent(32,0.715322);
    histPtr2->SetBinContent(33,0.3551669);
    histPtr2->SetBinContent(34,0.1833069);
    histPtr2->SetBinContent(35,0.09623345);
    histPtr2->SetBinContent(36,0.03405941);
    histPtr2->SetBinError(5,0.001310508);
    histPtr2->SetBinError(6,0.004185124);
    histPtr2->SetBinError(7,0.006635421);
    histPtr2->SetBinError(8,0.01065523);
    histPtr2->SetBinError(9,0.0169384);
    histPtr2->SetBinError(10,0.02752402);
    histPtr2->SetBinError(11,0.05911794);
    histPtr2->SetBinError(12,0.07964745);
    histPtr2->SetBinError(13,0.1126821);
    histPtr2->SetBinError(14,0.1711308);
    histPtr2->SetBinError(15,0.2865658);
    histPtr2->SetBinError(16,0.7144528);
    histPtr2->SetBinError(26,1.210666);
    histPtr2->SetBinError(27,0.3716564);
    histPtr2->SetBinError(28,0.2040013);
    histPtr2->SetBinError(29,0.1453484);
    histPtr2->SetBinError(30,0.09371448);
    histPtr2->SetBinError(31,0.03963412);
    histPtr2->SetBinError(32,0.02005007);
    histPtr2->SetBinError(33,0.01128522);
    histPtr2->SetBinError(34,0.006727023);
    histPtr2->SetBinError(35,0.004328359);
    histPtr2->SetBinError(36,0.001364499);
    histPtr2->SetEntries(87);
    histPtr2->SetLineColor(kBlue);
    histPtr2->SetMarkerColor(kBlue);

    histPtr2->Draw("same");



    TH1F*  hResult = new TH1F("hResult","hResult",40,xAxis);


    int simple_way = 0;
    if (simple_way == 1)
    {

// simple way    read about Sumw2() and options in TH1 class
//        hResult->Sumw2();
        hResult->SetBit(TH1::kIsAverage);
        histPtr->SetBit(TH1::kIsAverage);
        histPtr2->SetBit(TH1::kIsAverage);
        hResult->Add(histPtr,histPtr2,1,1);
    }
    else
    {
// sophisticated way, treat each bin individually


        if (histPtr->GetNbinsX() != histPtr2->GetNbinsX())
        {
            cout << " MUST HAVE same number of bins" << endl;
            return;
        }
        for (int i = 1; i < histPtr->GetNbinsX();i++)
        {
            double binContent1 =   histPtr->GetBinContent(i);
            double binContent2 =   histPtr2->GetBinContent(i);
            double newBinContent = 0;

            if (binContent1 == 0)
            {
                newBinContent = binContent2;
            }
            else
            {
                if (binContent2 == 0)
                {
                    newBinContent = binContent1;
                }
                else
                {
                    newBinContent = 0.5*(binContent1+binContent2);

                }
            }
// set content
            hResult->SetBinContent(i,newBinContent);

// set errors
            binContent1 = histPtr->GetBinError(i);
            binContent2 = histPtr2->GetBinError(i);
            newBinContent = sqrt(binContent1*binContent1 + binContent2*binContent2);
            hResult->SetBinError(i,newBinContent);
        }

        hResult->SetEntries(0.5*(histPtr->GetEntries()+histPtr2->GetEntries()));
    }

    hResult->SetLineColor(kRed);
    hResult->SetMarkerColor(kRed);

    hResult->Draw("same");

    c1->Modified();
    c1->cd();
    c1->SetSelected(c1);

//f->Close();
}