#include "TChain.h"
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TObjArray.h"
#include "TString.h"
#include "TSystem.h"

Double_t getDecayBR(const TString &name) {

    // Most branching ratios are taken from 2012 PDG.

    // Values for Delta+ are from QED calculations:
    // Probing Hadronic Structure with the decay Delta -> N l+ l-
    // Malcolm N. Butler
    // http://arxiv.org/pdf/nucl-th/9303015.pdf

    if(name.Contains("Delta+.")) {
	if(name.Contains(".pepem"))
	    return 4.2e-5;
	if(name.Contains(".pmpmm"))
	    return 3.e-7;
    }

    if(name.Contains("eta.")) {
        if(name.Contains(".gepem"))
           return 6.9e-3;
        if(name.Contains(".gmpmm"))
            return 3.1e-4;
        if(name.Contains(".epem"))
           return 5.6e-6;
        if(name.Contains(".mpmm"))
            return 5.8e-6;
    }

    if(name.Contains("eta'.") || name.Contains("etap")) {
	if(name.Contains(".gmpmm"))
            return 1.07e-4;
    }

    if(name.Contains("omega.")) {
        if(name.Contains(".pi0epem"))
            return 7.7e-4;
        if(name.Contains(".pi0mpmm"))
            return 1.3e-4;
        if(name.Contains(".epem"))
            return 7.28e-5;
        if(name.Contains(".mpmm"))
            return 9.0e-5;
    }

    if(name.Contains("pi0.")) {
        if(name.Contains(".gepem"))
            return 1.17e-2;
    }

    if(name.Contains("phi.")) {
        if(name.Contains(".epem"))
            return 2.97e-4;
        if(name.Contains(".mpmm"))
            return 2.87e-4;
    }

    if(name.Contains("rho0.")) {
        if(name.Contains(".epem"))
            return 4.7e-5;
        if(name.Contains(".mpmm"))
            return 4.55e-5;
    }

    if(name.Contains("jpsi.")) {
        if(name.Contains(".epem"))
            return 0.0594;
        if(name.Contains(".mpmm"))
            return 0.0594;
        if(name.Contains(".gepem"))
	    return 8.8e-3;
    }

    cout<<"No branching ratio found for file "<<name.Data()<<endl;

    return 0.;
}

Double_t getDecayMultiplicity(const TString &name) {

    if(name.Contains("25gev")) {
	if(name.Contains("Delta+."))
	    return 506.; // 337 * 3/2

	if(name.Contains("eta."))
	    return 36.;

	if(name.Contains("omega."))
	    return 38.;

	if(name.Contains("pi0."))
	    return 337.;

	if(name.Contains("phi."))
	    return 1.28;

	if(name.Contains("rho0."))
	    return 23.;

	if(name.Contains("jpsi."))
	    return 1.92e-5;
    } else {
	if(name.Contains("Delta+."))
	    return 506.; //?

	if(name.Contains("eta."))
	    return 16.;

	if(name.Contains("omega."))
	    return 19.;

	if(name.Contains("pi0."))
	    return 221.;

	if(name.Contains("phi."))
	    return 0.12;

	if(name.Contains("rho0."))
	    return 9.;

	if(name.Contains("jpsi."))
	    return 1.74e-7;
    }
    cout<<"No multiplicity found for file "<<name.Data()<<endl;
    return 0.;
}

void scaleHist(TH1F* const h, TString const& fileName) {
    // h: histogram to scale (will be modified)
    // fileName: Name of Pluto QA file. Used to derive branching ratio
    // and multiplicity from.

    if(!h) {
        cout<<"No histogram pointer."<<endl;
        return;
    }

    Double_t br    = getDecayBR(fileName);
    Double_t mult  = getDecayMultiplicity(fileName);
    if(br != 0. && mult != 0.) {
	h->Sumw2();
        h->Scale(br * mult / ((Double_t)h->GetEntries() * (Double_t)h->GetBinWidth(1)));
    }
}


Int_t pluto_qa_scale_hists(TString inFilePattern="./qa/*root") {
    // Scale histograms by branching ratio, multiplicity, number of events
    // and bin width.
    // inFilePattern: name of QA file(s).

    // Get list of input files
    TString inFiles = gSystem->GetFromPipe("ls -1 " + inFilePattern);
    TObjArray *inFilesList = inFiles.Tokenize("\n");
    if(inFilesList->GetEntries() == 0) {
	cout<<"Input files: "<<inFilePattern.Data()<<endl;
        return 1;
    }

    Int_t nEntries = inFilesList->GetEntries();

    TFile *f = NULL;
    TH1F  *h = NULL;
    const Int_t nHists = 14;
    // Name of histograms to scale.
    TString arrHists[nHists] = { "h_invmass", "h_p", "h_pt", "h_mt",
    "h_p1", "h_p2", "h_px", "h_py", "h_pz", "h_px_ll", "h_py_ll", "h_pz_ll",
    "h_y", "h_opang"};

    for (Int_t iFile = 0; iFile < nEntries; iFile++) {
        // Nomalize histograms by number of entries and bin width.
        TString fileName = ((TObjString*)inFilesList->At(iFile))->GetString();
        f = new TFile(fileName.Data(), "UPDATE");

        for(Int_t iHist = 0; iHist < nHists; iHist++) {
            h = (TH1F*)gDirectory->Get(arrHists[iHist]+";1");
            scaleHist(h, fileName);
            h->Write("", TObject::kOverwrite);
        }

        // Scale m_t - pole mass histogram.
        TH1F *h_mt_m0 = (TH1F*)gDirectory->Get("h_mt_m0;1");
        TH1F *h_y = (TH1F*)gDirectory->Get("h_y;1");
        Double_t normMtM0 = getDecayBR(fileName) * getDecayMultiplicity(fileName) /
            ((Double_t)h_mt_m0->GetEntries() * (Double_t)h_mt_m0->GetBinWidth(1) *
             (Double_t)h_y->GetEntries()  * (Double_t)h_y->GetBinWidth(1));
        h_mt_m0->Scale(normMtM0);

        f->Close();
    }

    return 0;
}