// Run with program with
//
//		root DrawHiIntegral.cpp



#include "TFile.h"
#include "TH1F.h"
#include "TVector3.h"
#include "TTreeReader.h"
#include "TTreeReaderValue.h"

#include "TCanvas.h"
#include "TGraph.h"
#include "TAxis.h"
#include "TPad.h"
#include "TSpectrum.h"
#include "TRandom.h"

#include <TString.h>
#include "TSystem.h"

#include <fstream>
#include <iostream>
using namespace std;

double peakIntegral;
double e_peakIntegral;
// Set range of peak
	// Attention: // if you specifiy a lower edge of a bin for the variable peakRight, the full bin will be included, if you want to exclude this bin, pick a value slightly below the edge
	double peakLeft = 4.3;
	double peakRight = 4.6;
	// Set range for background fit
	double backgoundLeft = 4.1;
	double backgoundRight = 4.8;

	Double_t background(Double_t *x, Double_t *par){
	// Define gap with peak to skip
	if (x[0]>par[3] && x[0]<par[4]){
		TF1::RejectPoint();
		return 0;
	}
	// return a second order polynomial
	return par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
}

// Function to fit the background around a peak, and then calculate the peak integral without the background
// @param peakIntegral		: variable to hold the final result (passed by reference)
// @param e_peakIntegral	: variable to hold the error of the result (passed by reference)
// @param h					: histogram to fit to
// @param peakLeft			: left edge of peak
// @param peakRight			: right edge of peak
// @param backgoundLeft		: left point to start background fit
// @param backgoundRight	: right point to end background fit
void ComputePeakIntegral(double &peakIntegral, double &e_peakIntegral, TH1F* h, double peakLeft, double peakRight, double backgoundLeft, double backgoundRight){

       TCanvas* cPeak = new TCanvas("cPeak","cPeak");
       h->Draw();
       h->GetXaxis()->SetRangeUser(2*backgoundLeft-backgoundRight,2*backgoundRight-backgoundLeft);

 /* ////////////// New Part Root Forum
    TF1 *f = new TF1("f", "gausn(0) + pol1(3)"); // normalized gaussian peak + linear background
    f->SetParNames("Area", "Mean", "Sigma", "p0", "p1");
    Double_t xmin = 4.2, xmax = 4.7;
    // set "reasonable" initial values for all parameters
    f->SetParameters(h->Integral(h->FindFixBin(xmin), h->FindFixBin(xmax)), // "Area"
                 (xmax + xmin) / 2., (xmax - xmin) / 10., // "Mean", "Sigma"
                 0., 0.); // "p0", "p1"
    h->Fit(f, "", "", xmin, xmax); // fit from "xmin" to "xmax"
    //////////////////////// */

	// Define a function to fit the backgound
	TF1* f = new TF1("f",background,backgoundLeft,backgoundRight,5); // xmin, xmax, number of parameters
	f->SetLineColor(kRed);
	// Give some starting parameters for the fit
	f->SetParameters(0,h->GetBinContent(peakLeft)-peakLeft*(h->GetBinContent(peakRight)-h->GetBinContent(peakLeft))/(peakRight-peakLeft));
	f->SetParameters(1,(h->GetBinContent(peakRight)-h->GetBinContent(peakLeft))/(peakRight-peakLeft));
	f->SetParameters(2,0);
	// Fix the position of the gap
	f->FixParameter(3,peakLeft);
	f->FixParameter(4,peakRight);
	// Fit the function in the given range (option R=apply range)
	TFitResultPtr fitResult = h->Fit("f","RS");
	TMatrixDSym covMatrix = fitResult->GetCovarianceMatrix();
	fitResult->Print("V");
	f->Draw("same");
	cPeak->Update();

	// Get the full peak integral
	double fullPeak = h->Integral(h->FindBin(peakLeft),h->FindBin(peakRight));
	cout<<"fullPeak = "<<fullPeak<<endl ;
	// Calculate the backgound integral from the fit and divide by binwidth w to get number of entries in this area
	//backgroundIntegral 	= int (p0  +  p1 * x  +  p2 * x^2) dx from b1 to b2 / w
							//= (p0 * x  +  1/2 * p1 * x^2  +  1/3 * p2 * x^3) from b1 to b2 / w
							//= p0 * (b2-b1)/w  +  1/2 * p1 * (b2^2/b1^2)/w  +  1/3 * p2 * (b2^3/b1^3)/w
							//= p0 * a  +  p1 * b  +  p2 * c
	double b1 = h->GetXaxis()->GetBinLowEdge(h->FindBin(peakLeft)); // left edge of first bin to have full bins and thus same range as integral over histogram
	double b2 = h->GetXaxis()->GetBinUpEdge(h->FindBin(peakRight)); // right edge of last bin
	double a = (b2-b1)/h->GetBinWidth(1);
	double b = 0.5*(b2*b2-b1*b1)/h->GetBinWidth(1);
	double c = 1./3.*(b2*b2*b2-b1*b1*b1)/h->GetBinWidth(1);
	double backgroundIntegral = f->GetParameter(0)*a+f->GetParameter(1)*b+f->GetParameter(2)*c;
	// Calculate squared error on background integral, as the parameters are correlated, we have to consider their covariance
	double e2_backgroundIntegral = a*a*covMatrix(0,0) + b*b*covMatrix(1,1) + c*c*covMatrix(2,2) + 2.*a*b*covMatrix(0,1) + 2.*a*c*covMatrix(0,2) + 2.*b*c*covMatrix(1,2);
	cout<<"backgroundIntegral from "<<b1<<" to "<<b2<<" = "<<backgroundIntegral<<" +- "<<sqrt(e2_backgroundIntegral)<<endl;

	// Calculate the peak as difference from the two
	peakIntegral = fullPeak - backgroundIntegral;

	// Get the error through error propagation
	e_peakIntegral = sqrt(fullPeak+e2_backgroundIntegral);

	// Keep the window open so that the user can check the fit until he/she presses a key
	std::cout<<"Insert any number and press enter to continue"<<std::endl;
	int dummy;
	std::cin >> dummy;

	// Delete canvas
	delete cPeak;

}

//Making fits.

void DrawPeakIntegral() {


// PRESETTINGS

	// Set path to root file
	std::string myPath = "../Integral/";

	// Set name of root file
	std::string myFileName = "Out_run_1.root";

	// Open the file containing the tree
	TFile *myFile = TFile::Open(std::string(myPath+myFileName).c_str());

	// Create a TTreeReader for the tree by passing the TTree's name and the TFile it is in
	TTreeReader myReader("Secondaries", myFile);

	// Set variables to access the data in the tree

	//The branch "ParticleID" contains a vector<int>; access them as vID
   TTreeReaderValue<vector<int>> vPDGcode (myReader, "ParticleID"); //New ADD

   // A histogram for Particle ID ranging from 1 to 100  with 200 bins
	TH1F* hPDGcode = new TH1F("hPDGcode", ";Number;#entries", 2500, -30, 30); //New ADD


	// The branch "ParticleEnergy" contains a vector<double>; access them as vEnergies
	TTreeReaderValue<vector<double>> vEnergies(myReader, "ParticleEnergy");

	// Create histograms for the values we read

	// A histogram for the energies ranging from 0 MeV to 20 MeV with 200 bins
	TH1F* hEnergyGamma = new TH1F("hEnergyGamma", ";E in MeV;#entries", 600, 0,10);

	TH1F* hEnergyAll = new TH1F("hEnergyAll", ";E in MeV;#entries", 600, 0, 10);


// READ THE TREE

	// Loop over all entries of the TTree.

    while (myReader.Next()) {
    int energySize = vEnergies->size();

		for(int i = 0; i<energySize; ++i){

			double myEnergy= vEnergies->at(i);
			int myPDG= vPDGcode->at(i);

			hEnergyAll->Fill(myEnergy);
			if (myPDG==22){
				hEnergyGamma->Fill(myEnergy);
				}


			}
    }

  ComputePeakIntegral(peakIntegral, e_peakIntegral, hEnergyGamma, peakLeft, peakRight, backgoundLeft, backgoundRight);


  std::cout << "#####$$$$$   peakIntegral = " << peakIntegral << std::endl;

  std::cout << "#####$$$$$  e_peakIntegral = " << e_peakIntegral << std::endl;




}