#include<iostream>
#include<cmath>
#include<fstream>
#include<string>
#include<TCanvas.h>
#include <TFile.h>
#include<TApplication.h>
#include<TGraphErrors.h>
#include<TH1D.h>
#include<TF1.h>
#include<TStyle.h>
#include<TAxis.h>
#include<TLatex.h>
#include<TMath.h>
#include<TTree.h>
#include<cstdlib>
#include<stdlib.h>
#include<TRandom.h>
#include "TH2F.h"
#include "TCanvas.h"
#include "TPad.h"
#include<cstring>
#include<sstream>  
#include<iomanip>

using namespace std;

double SIGMA(double E, int j) {

	double A[2] = {0,     0};
  	double B[2] = {0.05, 0.05};
  	double C[2] = {0,     0};
  	return A[j] + B[j]*sqrt(E) + C[j]*E;
}

void sigmatree() {

  	TRandom *r = new TRandom();
  	const int Ndet = 2;

  		std::vector<double> *e0  = new vector <double>();
  			std::vector<double> *e  = new vector <double>();
  
  	TFile *RootIn = TFile::Open("Sim.root","READ");
  	TTree *hIn = (TTree*) (RootIn->Get("Tree1"));
  	TFile *RootOut = new TFile("SimRes.root","RECREATE");
  	TTree *hOut = new TTree("Tree1", "fwhm on");
    		hIn->SetBranchAddress("Edep", &e0);
    		hOut->Branch("Edep", &e);
  	

  	long int NEvents = hIn->GetEntries();
  	double sigma;
  	cout << "Number of Events: " << NEvents << endl;
  	long int Step = (long int)(NEvents/10);
  	for(long int i=0; i<NEvents; i++) {
    		if (i%Step==0) cout << (int)(i*1e2/NEvents) << "% events processed.." << endl;
    		hIn->GetEntry(i);
    		if (e0->size()>0){
    		for (int j=0; j<e0->size(); j++) {
      			sigma = SIGMA(e0->at(j),j);
      			e->push_back(e0->at(j)+r->Gaus(0,sigma)); 
      		    }
    		}
    	
    	hOut->Fill();
  	}
  	hOut->Write();
  	RootIn->Close();
  	RootOut->Close();
}