View of /groups/HadCalib/JetLevelCalib/ReadHistogramms.C

#include "ReadHistogramms.h"
#include "TChain.h"
#include "TString.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TH2.h"
#include <iostream>
#include "TROOT.h"
#include "TApplication.h"
#include <vector>
#include <math.h>
#include <fstream>
#include <string>
#include <sstream>
#include <iomanip>
#include "TStyle.h"

// constructor

using namespace std;

ReadHistogramms::ReadHistogramms(TString prefix, TString file){

	m_prefix = prefix;
	std::cout<<"\n\n------> Root file prefix : "<< m_prefix <<std::endl;
	
	// open the histo file	
	m_file = TFile::Open( Form("%s.root",file.Data()) ,"READ");
		std::cout<<"\n\n------> Root file prefix : "<< ( Form("%s",prefix.Data()) ) <<std::endl;

	// retrieve info in the ntuple	
	TChain *t = new TChain( Form("%s",file.Data()) );
	t->Add( Form("%s.root",file.Data()) );

	//	float * EtaBins;
	//float * EnergyBins;

	int nEtaBins=0;
	int nEnergyBins=0;
	
	//	std::vector<double>  m_EtaBins;
	//	std::vector<double>  m_EnergyBins;

	t->SetBranchAddress("NumberEtaBins", &nEtaBins );
	t->SetBranchAddress("NumberEnergyBins", &nEnergyBins );
	
	t->GetEntry(0);
	const Int_t constnEtaBins=nEtaBins;
	const Int_t constnEnergyBins=nEnergyBins;

	float EtaBins[constnEtaBins];
	float EnergyBins[constnEnergyBins];

	t->SetBranchAddress("EtaBins", EtaBins );
	t->SetBranchAddress("EnergyBins", EnergyBins );

	t->GetEntry(0);

	m_nEtaBins=nEtaBins;
	m_nEnergyBins=nEnergyBins;

	for (int iEta=0; iEta<m_nEtaBins; ++iEta){
	  	 std::cout << EtaBins[iEta] << std::endl;
		 m_EtaBins.push_back( EtaBins[iEta]) ;		
		}	
	 for (int iEne=0; iEne<m_nEnergyBins; ++iEne){
	   std::cout << EnergyBins[iEne] << std::endl;
	   m_EnergyBins.push_back( EnergyBins[iEne]) ;	
	 }

	delete t;

}

float ReadHistogramms::getEtaBinUp(int bin){
    return m_EtaBins.at(bin+1);
}
float ReadHistogramms::getEtaBinDown(int bin){
    return m_EtaBins.at(bin);
}
float ReadHistogramms::getEnergyBinUp(int bin){
    return m_EnergyBins.at(bin+1);
}
float ReadHistogramms::getEnergyBinDown(int bin){
    return m_EnergyBins.at(bin);
}
//destructor
ReadHistogramms::~ReadHistogramms(){
        std::cout<<"Destruction..."<<std::endl;
        m_file->Close();
}



TGraphErrors* ReadHistogramms::plotVsEta(std::string type, int bin, bool fit) {

	TGraphErrors *g = new TGraphErrors();
	Double_t yAxis=0, dyAxis=0;
	Double_t xAxis=0, dxAxis=0;
	TString axisTitle;
	Double_t nEntries = 0;			
	int status = 0;
	double Mean[2]={0}, Reso[2]={0};
	std::string variable;
	
	int iEne =0;
	// select eta bin
	iEne=bin;
	m_title = Form("%1.2f < E_T < %1.2f", m_EnergyBins[iEne], m_EnergyBins[iEne+1] );
	std::cout<<"\tLooking at bin "<< bin <<" : "<< m_title.Data() <<std::endl;
	
	// case 1 : linearity or resolution
	if ( "linearity"==type || "resolution"==type ) variable = "RatioEnergy_eta";
	// case 2 : spatial resolution in eta 
	else if ( "spatial resolution eta" == type ) variable = "SpatialResolution_eta";	
	// case 3 : spatial resolution in phi 
	else if ( "spatial resolution phi" == type ) variable = "SpatialResolution_phi";	
	// case 4 : spatial resolution in delta R 
	else if ( "spatial resolution deltaR" == type ) variable = "SpatialResolution_deltar";	
	else if ( "DeltaEnergy" == type ) variable = "DeltaEnergy_eta";		
	else return g;
	
	std::cout<<"\n\tPlot "<< type <<" vs eta : "<< m_title.Data() <<std::endl;
	
	int iPoint=0;
	for (int iEta=0; iEta<m_nEtaBins-1; ++iEta){
		xAxis=0; 
		dxAxis=0;
		yAxis=0; 
		dyAxis=0;
							
		TString variableName = Form("%s_%s_eta%i_e%i", variable.c_str(), m_prefix.c_str(), iEta, iEne );
		TH1F *histo = (TH1F*)m_file->Get( variableName.Data() );

		nEntries = histo->GetEntries();			

		if ( nEntries>0 ) {
			
			xAxis  = 0.5*(m_EtaBins[iEta+1]+m_EtaBins[iEta]);
			dxAxis = 0.5*(m_EtaBins[iEta+1]-m_EtaBins[iEta]);
								
			if ( "resolution" == type ){
				status = fitGauss( histo, Mean, Reso, fit);
				yAxis = Reso[0];
				dyAxis = Reso[1];
				axisTitle=";|#eta_{truth}|; #sigma_{E}/E_{jet} (%)";			
			} else {
			status = fitGauss( histo, Mean, Reso, fit );
				yAxis = Mean[0];
				dyAxis = Mean[1];		
				if ("linearity" == type ) axisTitle= ";|#eta_{truth}|; E_{jet}/E_{truth}";			
				else if ("spatial resolution eta" == type ) axisTitle= ";#eta_{truth}; #eta_{rec}-#eta_{truth}";			
				else if ("spatial resolution phi" == type ) axisTitle= ";#eta_{truth}; #phi_{rec}-#phi_{truth}";
				else if ("spatial resolution phi" == type ) axisTitle= ";#eta_{truth}; #phi_{rec}-#phi_{truth}";		
				else if ( "DeltaEnergy" == type ) axisTitle= ";#eta_{truth}; E_{truth}-E_{reco} [MeV]";		
			}
			g->SetPoint(iPoint, xAxis, yAxis );
			g->SetPointError(iPoint, dxAxis, dyAxis );
			++iPoint;
		}
		TString line = Form("%s : eta [%1.2f; %1.2f]\tn=%3.0f\tx=%1.1f +/- %1.1f\ty=%1.3f +/- %1.3f\t%i", 
			variableName.Data(), m_EtaBins[iEta], m_EtaBins[iEta+1], nEntries, xAxis, dxAxis, yAxis, dyAxis, status );
		std::cout<< line.Data() <<std::endl;
	}
	g->SetTitle(m_title + axisTitle );
	return g;
}

TH1F* ReadHistogramms::plotDeltaR(std::string type, int etabin,int ebin, bool fit) {

  //TGraphErrors *g = new TGraphErrors();
	Double_t yAxis=0, dyAxis=0;
	Double_t xAxis=0, dxAxis=0;
	TString axisTitle;
	Double_t nEntries = 0;			
	int status = 0;
	std::string variable;
	
	m_title = Form("%1.2f < E_T < %1.2f", m_EnergyBins[ebin], m_EnergyBins[ebin+1] );
	std::cout<<"\tLooking at bin "<< ebin <<" : "<< m_title.Data() <<std::endl;
	std::string m_title2 = Form(" %1.2f < #eta < %1.2f", m_EtaBins[etabin], m_EtaBins[etabin+1] );
	std::cout<<"\tLooking at bin "<< etabin <<" : "<< m_title2 <<std::endl;
	m_title = "title";

	TString variableName = Form("SpatialResolution_deltaR_%s_eta%i_e%i", m_prefix.c_str(), etabin, ebin );
	
	TH1F *histo = (TH1F*)m_file->Get( variableName.Data() );

	axisTitle=";#Delta R;Probability";		

	histo->SetTitle(m_title + axisTitle );
	return histo;
}

TH2F* ReadHistogramms::plotCluster3D(TString type) {

    bool debug =false; // if we set debug to true, we will have more output

    const Int_t nEtaBins = m_nEtaBins;  // m_nEtaBins is a data member of the class ReadHistogramms. It is read out from the tree in the function ReadHistogramms::ReadHistogramms(TString prefix, TString file).
// NOTE: if you are just using a plot macro/normal programm, you can copy&paste the code from ReadHistogramm in order to read out the m_nEtaBins...

    const Int_t nEnergyBins = m_nEnergyBins;// same as above

// in the following the vectors m_EtaBins (again, datamembers of ReadHistogramms and read out in ReadHistogramms::ReadHistogramms are copied into float arays, because root is stupid and we need float arrays in order to be able to create TH2F histogramm axis. (in a pure 1-D plotting, you'd not need that).
        float xbins[nEtaBins];
        float ybins[nEnergyBins];
      double xbins2[nEtaBins];
        double ybins2[nEnergyBins];

        for (int i=0; i<nEtaBins;i++) {
	  xbins2[i] = m_EtaBins[i];
	  xbins[i] = m_EtaBins[i];
	  if(debug) cout << "xbins : " << xbins[i] << endl;
        }
        for (int j=0; j<nEnergyBins;j++){
	  ybins2[j] = m_EnergyBins[j];
	  ybins[j] = m_EnergyBins[j];
	  if(debug) cout  << "ybins : " << ybins[j] << endl;
        }


// so here we are creating the 2D histogramm, we want to fill with quantities from the 1D histogramms later
	TH2F *histo2D= new TH2F(type.Data(),type.Data(),nEtaBins-1 ,xbins ,nEnergyBins-1 ,ybins); // here use eta_reco and energy_reco;

	cout << "before loop" << endl;

// now we are looping over all 1D cluster/Jet moment histogramms (we now their number because there is always 1 histogramm for each bin in eta and energy
	for (int etabin=0;etabin<nEtaBins-1;etabin++) {
	for (int ebin=0;ebin<nEnergyBins-1;ebin++) {
	    cout << "ebin: " << ebin <<  " etabin: " << etabin <<endl;

	    TString variableName = Form("%s%s_eta%i_e%i", type.Data(),m_prefix.c_str(), etabin, ebin );// here we are creatng the name of the histogramm we want to read out, the form is TYPE of the histogramm (i.e. ClusterMoment_isol or something else..., depending on the names of the histogramms), Prefix is the name of the Jet Container, i.e. for example Cone7CaloJets. etabin and ebin are e.g. 0 and 0 for the first bin in eta/energy... i.e. the very first histogramm - HERE on would like to use only some selected 9 bins in eta/energy

	cout << variableName<<endl;
	TH1I *histo = (TH1I*)m_file->Get( variableName.Data() ); // here we read out the actual histgramm from the file

	if ( histo->GetEntries()>0) histo2D->SetBinContent(etabin,ebin,histo->GetMean()); // here we use the mean from the histogramm and set the respective bin in the 2D histogramm.
	}
	}
//	axisTitle=";#Delta R;Probability";		

	//histo->SetTitle(m_title + axisTitle );
	//return histo2D;
	
	gStyle->SetPalette(1);
        gStyle->SetPadRightMargin(0.19);
	gPad->SetLogy();
	histo2D->Draw("LEGO2Z");
	histo2D->SetTitle(Form("%s; #eta_{truth jet} ;E_{truth jet} [GeV];mean of cluster distribution" ,type.Data()));
	histo2D->GetXaxis()->SetTitleOffset(2.0);
	histo2D->GetYaxis()->SetTitleOffset(2.0);
	return histo2D;
}

TGraphErrors* ReadHistogramms::plotVsEnergy(std::string type, int bin, bool fit){
	TGraphErrors *g = new TGraphErrors();
	Double_t xAxis=0, dxAxis=0;
	Double_t yAxis=0, dyAxis=0;
	TString axisTitle;
	
	Double_t nEntries = 0;			
	int status =0;
	int iEta=0;
	Double_t Mean[2] = {0}, Reso[2] = {0};

	std::cout<<"\n\t"<< type.c_str() <<" vs E : "<<std::endl;
	
	// select eta bin
	iEta=bin;
	m_title = Form("%1.2f < #eta < %1.2f", m_EtaBins[iEta], m_EtaBins[iEta+1] );
	std::cout<<"\tLooking at bin "<< bin <<" : "<< m_title.Data() <<std::endl;
	
	// loop over all energy bins
	int nPoint = 0;
	for (int iEne=0; iEne<m_nEnergyBins-1; ++iEne){
		xAxis=0; 
		dxAxis=0;
		yAxis=0; 
		dyAxis=0;
							
		TString variableName = Form("RatioEnergy_eta_%s_eta%i_e%i", m_prefix.c_str(), iEta, iEne );
		std::cout << variableName << std::endl;
		TH1F *histo = (TH1F*)m_file->Get( variableName.Data() );

		nEntries = histo->GetEntries();			
		
		if ( nEntries>10 ) {
			status = fitGauss( histo, Mean, Reso, fit);
			xAxis = 0.5*(m_EnergyBins[iEne+1]+m_EnergyBins[iEne]);
			dxAxis = 0.5*(m_EnergyBins[iEne+1]-m_EnergyBins[iEne]);
								
			if ( "linearity" == type ){
				yAxis = Mean[0];
				dyAxis = Mean[1];
				axisTitle=";E_{truth} (GeV); E_{jet}/E_{truth}";
			} else if ( "resolution" == type  ){
			    yAxis = (Reso[0])/(Mean[0]);
				dyAxis = Reso[1];		
				axisTitle=";E_{truth} (GeV); #sigma_{E}/E_{jet} (%)";			
			} 
			else if ( "DeltaEnergy" == type ) {
				yAxis = Mean[0];
				dyAxis = Mean[1];
				axisTitle=";E_{truth} (GeV); E_{truth}-E_{reco} (MeV)";
			}
			else {
				std::cout<<"\t" <<type.c_str() <<" is not a valid option."<<std::endl;
				return g;
			}		
			g->SetPoint(nPoint, xAxis, yAxis );
			g->SetPointError(nPoint, dxAxis, dyAxis/100. );
			++nPoint;
		}
		TString line = Form("%s : E [%4.0f; %4.0f]\tn=%3.0f\tx=%1.1f +/- %1.1f\ty=%1.3f +/- %1.3f\t%i", 
			variableName.Data(), m_EnergyBins[iEne], m_EnergyBins[iEne+1], nEntries, xAxis, dxAxis, yAxis, dyAxis, status );
		std::cout<< line.Data() <<std::endl;
	}
	g->SetTitle(m_title + axisTitle );

 if ( "resolution" == type ) {

     TF1 *fitResolution = new TF1("func", "sqrt(([0]*[0])/x + [1]*[1]+ ([2]*[2])/(x*x)  )", 100, 3000);//+ ([2]*[2])/(x*x) 
                        fitResolution->SetLineColor(1);
                        if (1==bin) fitResolution->SetLineColor(1);
			fitResolution->SetLineWidth(2);
			fitResolution->SetLineStyle(2);
			fitResolution->SetParameter(0, 95.5);
			fitResolution->SetParameter(1, 4.9);
			fitResolution->SetParameter(2, 2.);
                        g->Fit(fitResolution, "q+");
			TString title = Form("a=%2.2f; b=%2.2f; c=%2.2f", fitResolution->GetParameter(0), fitResolution->GetParameter(2), fitResolution->GetParameter(1));
					std::cout << title << std::endl;
			title = Form("Sigma(a)=%2.2f; Sigma(b)=%2.2f", fitResolution->GetParError(0), fitResolution->GetParError(1));
			std::cout << title << std::endl;
			
			std::cout << "**************************" << std::endl;
			std::cout << "& " <<  Form("%2.1f \\pm %2.1f",  fitResolution->GetParameter(0), fitResolution->GetParError(0))  ;
std::cout << "& " <<  Form("%2.1f \\pm %2.1f",  fitResolution->GetParameter(1), fitResolution->GetParError(1))  ;
	std::cout << "& " <<  Form("%2.1f \\pm %2.1f",  fitResolution->GetParameter(2), fitResolution->GetParError(2))  ;		
	std::cout << "& " <<  Form("%2.1f/%d",  fitResolution->GetChisquare(), fitResolution->GetNDF())  ;
			std::cout << " \\hline " << std::endl;
			std::cout << "**************************" << std::endl;

}

	return g;
}	


TGraphErrors* ReadHistogramms::plotVsEnergy(std::string type, int bin){
	TGraphErrors *g = new TGraphErrors();
	Double_t xAxis=0, dxAxis=0;
	Double_t yAxis=0, dyAxis=0;
	TString axisTitle;
	
	Double_t nEntries = 0;			
	//Double_t mean = 0, dmean = 0;
	//Double_t rms = 0, drms = 0;
	int status =0;
	int iEta=0;
	Double_t Mean[2] = {0}, Reso[2] = {0};

	std::cout<<"\n\t"<< type.c_str() <<" vs E : "<<std::endl;
	if ( m_nEtaBins !=3 ) {
		std::cout<<"\tImpossible to plot linearity vs E for this collection of histograms !"<<std::endl;
		return g;
	}
	
	// select eta bin
	if ( bin==0 ) iEta=0;
	else if (bin==1) iEta=2;
	else return g;	
	m_title = Form("%1.2f < #eta < %1.2f", m_EtaBins[iEta], m_EtaBins[iEta+1] );
	std::cout<<"\tLooking at bin "<< bin <<" : "<< m_title.Data() <<std::endl;
	
	// loop over all energy bins
	int nPoint = 0;
	for (int iEne=0; iEne<m_nEnergyBins; ++iEne){
		xAxis=0; 
		dxAxis=0;
		yAxis=0; 
		dyAxis=0;
							
		TString variableName = Form("RatioEnergy_eta_%s_eta%i_e%i", m_prefix.c_str(), iEta, iEne );
		TH1F *histo = (TH1F*)m_file->Get( variableName.Data() );

		nEntries = histo->GetEntries();			
		
		if ( nEntries>30 ) {
			status = fitGauss( histo, Mean, Reso);
			xAxis = 0.5*(m_EnergyBins[iEne+1]+m_EnergyBins[iEne]);
			dxAxis = 0.5*(m_EnergyBins[iEne+1]-m_EnergyBins[iEne]);
								
			if ( "linearity" == type ){
				yAxis = Mean[0];
				dyAxis = Mean[1];
				axisTitle=";E_{truth} (GeV); E_{jet}/E_{truth}";
			} else if ( "resolution" == type ){
				yAxis = Reso[0]/(100*Mean[0]);
				dyAxis = Reso[1];		
				axisTitle=";E_{truth} [GeV]; #sigma_{E}/E_{jet} (%)";			
			} else {
				std::cout<<"\t" <<type.c_str() <<" is not a valid option."<<std::endl;
				return g;
			}		
			g->SetPoint(nPoint, xAxis, yAxis );
			g->SetPointError(nPoint, dxAxis, dyAxis );
			++nPoint;
		}
		TString line = Form("%s : E [%4.0f; %4.0f]\tn=%3.0f\tx=%1.1f +/- %1.1f\ty=%1.3f +/- %1.3f\t%i", 
			variableName.Data(), m_EnergyBins[iEne], m_EnergyBins[iEne+1], nEntries, xAxis, dxAxis, yAxis, dyAxis, status );
		std::cout<< line.Data() <<std::endl;
	}
	
	g->SetTitle( axisTitle );

 if ( "resolution" == type ) {

TF1 *fitResolution = new TF1("func", "sqrt( [0]*[0]/x + [1]*[1]/(x*x) + [2]*[2] )", 20, 10000);
                        fitResolution->SetLineColor(2);
                        if (1==bin) fitResolution->SetLineColor(9);
                        fitResolution->SetLineWidth(2);
                        fitResolution->SetLineStyle(2);
                        //fitResolution->SetParLimits(0, 0, 200);
                        fitResolution->SetParameter(0, 60);
			fitResolution->SetParameter(1, 600);
                        fitResolution->SetParameter(2, 3);
                        g->Fit(fitResolution, "q");
                        std::cout << "**************************" << std::endl;
			m_title = Form("a= %2.2f +/- %2.2f;  b= %2.2f +/- %2.2f;   c=%2.2f +/- %2.2f, Chi2/NDF%2.2f / %d",fitResolution->GetParameter(0),fitResolution->GetParError(0), fitResolution->GetParameter(1),fitResolution->GetParError(1), fitResolution->GetParameter(2),fitResolution->GetParError(2),fitResolution->GetChisquare(),fitResolution->GetNDF());
                        std::cout << m_title << std::endl;
			std::cout << fitResolution->GetNDF()<<  "**************************" << std::endl;
}

	
	if ( "resolution" == type ) g->GetYaxis()->SetRangeUser(2, 22);
	if ( "linearity" == type ) g->GetYaxis()->SetRangeUser(0.5, 1.3);

	return g;
}	

TGraphErrors* ReadHistogramms::plotVsEta(std::string type, int bin){
	TGraphErrors *g = new TGraphErrors();
	Double_t yAxis=0, dyAxis=0;
	Double_t xAxis=0, dxAxis=0;
	TString axisTitle;
	Double_t nEntries = 0;			
	int status = 0;
	double Mean[2]={0}, Reso[2]={0};
	std::string variable;
	
	int iEne =0;
	// select eta bin
	iEne=bin;

	m_title = Form("%3.0f < E_{T} < %3.0f", m_EnergyBins[iEne], m_EnergyBins[iEne+1] );
	std::cout<<"\tLooking at bin "<< bin <<" : "<< m_title.Data() <<std::endl;
	
	// case 1 : linearity or resolution
	if ( "linearity"==type || "resolution"==type ) variable = "RatioEnergy_eta";
	// case 2 : spatial resolution in eta 
	else if ( "spatial resolution eta" == type ) variable = "SpatialResolution_eta";	
	// case 3 : spatial resolution in phi 
	else if ( "spatial resolution phi" == type ) variable = "SpatialResolution_phi";
	else if ( "DeltaEnergy" == type ) variable = "DeltaEnergy_eta";			
	else return g;
	
	std::cout<<"\n\tPlot "<< type <<" vs eta : "<< m_title.Data() <<std::endl;
	
	int iPoint=0;
	for (int iEta=0; iEta<m_nEtaBins; ++iEta){
		xAxis=0; 
		dxAxis=0;
		yAxis=0; 
		dyAxis=0;
							
		TString variableName = Form("%s_%s_eta%i_e%i", variable.c_str(), m_prefix.c_str(), iEta, iEne );
		TH1F *histo = (TH1F*)m_file->Get( variableName.Data() );

		nEntries = histo->GetEntries();			

		if ( nEntries>30 ) {
			
			xAxis  = 0.5*(m_EtaBins[iEta+1]+m_EtaBins[iEta]);
			dxAxis = 0.5*(m_EtaBins[iEta+1]-m_EtaBins[iEta]);
								
			if ( "resolution" == type ){
				status = fitGauss( histo, Mean, Reso);
				yAxis = Reso[0];
				dyAxis = Reso[1];
				axisTitle=";|#eta_{truth}|; #sigma_{E}/E_{jet} (%)";							
			} else {
				status = fitGauss( histo, Mean, Reso, -1, 1 );
				yAxis = Mean[0];
				dyAxis = Mean[1];		
				if ("linearity" == type ) axisTitle= ";|#eta_{truth}|; E_{jet}/E_{truth}";			
				else if ("spatial resolution eta" == type ) axisTitle= ";#eta_{truth}; #eta_{rec}-#eta_{truth}";			
				else if ("spatial resolution phi" == type ) axisTitle= ";#eta_{truth}; #phi_{rec}-#phi_{truth}";			
			}
			g->SetPoint(iPoint, xAxis, yAxis );
			g->SetPointError(iPoint, dxAxis, dyAxis );
			++iPoint;
		}
		TString line = Form("%s : eta [%1.2f; %1.2f]\tn=%3.0f\tx=%1.1f +/- %1.1f\ty=%1.3f +/- %1.3f\t%i", 
			variableName.Data(), m_EtaBins[iEta], m_EtaBins[iEta+1], nEntries, xAxis, dxAxis, yAxis, dyAxis, status );
		std::cout<< line.Data() <<std::endl;
	}
	g->SetTitle( axisTitle );

	if ( "resolution" == type ) g->GetYaxis()->SetRangeUser(2, 22);
	else if ( "linearity" == type ) g->GetYaxis()->SetRangeUser(0.5, 1.3);
	else g->GetYaxis()->SetRangeUser(-0.005, 0.005);	
	return g;
}	

TGraphAsymmErrors* ReadHistogramms::plotEffVsEnergy(std::string type, int bin){
  // note: this gives back efficiency agains energy for 1 bin in eta

	TGraphAsymmErrors *g = new TGraphAsymmErrors();	
	TString axisTitle;

	/// ATTENTION: Programme apparently works with zeroth bins....
	std::cout<<"\n\t"<< type.c_str() <<" vs E : "<<std::endl;
	if ( m_nEtaBins <1 ) {
		std::cout<<"\tImpossible to plot linearity vs E for this collection of histograms !"<<std::endl;
		return g;
	}
  
	m_title = Form("%1.2f < #eta < %1.2f", m_EtaBins[bin], m_EtaBins[bin+1] );

	std::cout<<"\tLooking at bin "<< bin <<" : "<< m_title.Data() <<std::endl;
	TH2 *matched;
	TH2 *total;
	
	if ( "efficiency" == type ){		
		TString variableName = Form("%sNumberMatchedTrueJets", m_prefix.c_str() );
		matched = (TH2I*)m_file->Get( variableName.Data() );
		variableName = Form("%sNumberTrueJets", m_prefix.c_str() );
		total = (TH2I*)m_file->Get( variableName.Data() );
		axisTitle="; Energy_{truth} (GeV);Efficiency";
  } else if ( "purity" == type ){		
		TString variableName = Form("%sNumberMatchedRecoJets", m_prefix.c_str() );
		matched = (TH2*)m_file->Get( variableName.Data() );
		variableName = Form("%sNumberRecoJets", m_prefix.c_str() );
		total = (TH2*)m_file->Get( variableName.Data() );
		axisTitle="; Energy_{jet} (GeV);Purity";
	} else {
		std::cout<<"\t" <<type.c_str() <<" is not a valid option."<<std::endl;
		return g;
	}
  
  TH1D *Nmatched = matched->ProjectionY("Nmatched", bin+1, bin+1);
  TH1D *Ntotal = total->ProjectionY("Ntotal", bin+1, bin+1);		
  g->BayesDivide(Nmatched, Ntotal);
  g->SetTitle(m_title + axisTitle);

	return g;
}	

TGraphAsymmErrors* ReadHistogramms::plotEffVsEta(std::string type, int bin){
	TGraphAsymmErrors *g = new TGraphAsymmErrors();

	TString axisTitle;
	int iEne=bin;

	/// ATTENTION: Programme apparently works with zeroth bins....
	std::cout<<"\n\t"<< type.c_str() <<" vs E : "<<std::endl;
	if ( m_nEnergyBins <1 ) {
		std::cout<<"\tImpossible to plot efficiency vs E for this collection of histograms !"<<std::endl;
		return g;
	}
  
	m_title = Form("%3.0f < E_{T} < %3.0f", m_EnergyBins[iEne], m_EnergyBins[iEne+1] );
	std::cout<<"\tLooking at bin "<< bin <<" : "<< m_title.Data() <<std::endl;
	TH2 *matched;
	TH2 *total;
	
	if ( "efficiency" == type ){		
		TString variableName = Form("%sNumberMatchedTrueJets", m_prefix.c_str() );
		matched = (TH2I*)m_file->Get( variableName.Data() );
		variableName = Form("%sNumberTrueJets", m_prefix.c_str() );
		total = (TH2I*)m_file->Get( variableName.Data() );
		axisTitle="; |#eta_{truth}|;Efficiency";
	} else if ( "purity" == type ){		
		TString variableName = Form("%sNumberMatchedRecoJets", m_prefix.c_str() );
		matched = (TH2*)m_file->Get( variableName.Data() );
		variableName = Form("%sNumberRecoJets", m_prefix.c_str() );
		total = (TH2*)m_file->Get( variableName.Data() );
		axisTitle="; |#eta_{jet}|;Purity";
	} else {
		std::cout<<"\t" <<type.c_str() <<" is not a valid option."<<std::endl;
		return g;
	}
  
	TH1D *Nmatched = matched->ProjectionX("Nmatched", bin+1, bin+1);
	TH1D *Ntotal = total->ProjectionX("Ntotal", bin+1, bin+1);
		
	g->BayesDivide(Nmatched,Ntotal);
	g->SetTitle(m_title + axisTitle);
	
	TCanvas *c = new TCanvas();
	c->Divide(2,2);
	c->cd(1); total->Draw();
	c->cd(2); matched->Draw();
	c->cd(3); Nmatched->Draw();
	c->cd(4); Ntotal->Draw();
	c->Print("tmp.root");
	
 
	return g;
}	

// private function : compute the mean & resolution using a gaussian fit between +/- 2 sigma
int ReadHistogramms::fitGauss( TH1F* hist, double *mean, double *reso, double xMin, double xMax){	
	//return 0 if no entries are found in the histogram
	
	if (hist->GetEntries()==0){
		mean[0] = 0;
		mean[1] = 0;
		reso[0] = 0;            
		reso[1] = 0;
		return -1;                              	
	}	
	TCanvas *c1 = new TCanvas();
	TF1 *g = new TF1("G", "gaus", -2, 2);            
	//hist->Fit(g, "qL", "", xMin, xMax);      
	hist->Fit(g, "q");      
	double  mean0 = g->GetParameter(1);
	double sigma0 = g->GetParameter(2);
   
	g->SetParLimits(1, mean0-3*sigma0, mean0+3*sigma0 );   
	//int status = hist->Fit(g, "qL", "", mean0-2*sigma0, mean0+2*sigma0);      
	int status = hist->Fit(g, "q", "", mean0-2*sigma0, mean0+2*sigma0);      
	mean[0] = g->GetParameter(1);
	mean[1] = g->GetParError(1);	
	double sigma = g->GetParameter(2);            
	double dsigma = g->GetParError(2);                              
	
	reso[0] = 100*sigma  ;// /mean[0];
	reso[1] = reso[0]*sqrt( pow(mean[1]/mean[0], 2) + pow(dsigma/sigma, 2) );
	
	delete g;
	delete c1;     
	return status;            
}

double ReadHistogramms::GetFluctTerm(std::string type, int bin, bool fit,double *sampterm){
  
	TGraphErrors *g1 = new TGraphErrors();
	double a=0;
	Double_t xAxis=0, dxAxis=0;
	Double_t yAxis=0, dyAxis=0;
	TString axisTitle;
	
	Double_t nEntries = 0;			
	//Double_t mean = 0, dmean = 0;
	//Double_t rms = 0, drms = 0;
	int status =0;
	int iEta=0;
	Double_t Mean[2] = {0}, Reso[2] = {0};

	std::cout<<"\n\t"<< type.c_str() <<" vs E : "<<std::endl;
	//	if ( m_nEtaBins !=3 ) {
	//std::cout<<"\tImpossible to plot linearity vs E for this collection of histograms !"<<std::endl;
	//return g;
	//}
	
	// select eta bin
	iEta=bin;
	//	else if (bin==) iEta=2;
	//	else return g;	
	m_title = Form("%1.2f < #eta < %1.2f", m_EtaBins[iEta], m_EtaBins[iEta+1] );
	std::cout<<"\tLooking at bin "<< bin <<" : "<< m_title.Data() <<std::endl;
	
	// loop over all energy bins
	int nPoint = 0;
	for (int iEne=0; iEne<m_nEnergyBins; ++iEne){
		xAxis=0; 
		dxAxis=0;
		yAxis=0; 
		dyAxis=0;
							
		TString variableName = Form("RatioEnergy_eta_%s_eta%i_e%i", m_prefix.c_str(), iEta, iEne );
		TH1F *histo = (TH1F*)m_file->Get( variableName.Data() );

		nEntries = histo->GetEntries();			
		//mean = histo->GetMean();
		//rms = histo->GetRMS();
		//dmean = 0;
		//drms = 0;
		
		if ( nEntries>100 ) {
			status = fitGauss( histo, Mean, Reso, fit);
		   //dmean = rms/sqrt(nEntries);
		   //drms = rms/2*sqrt(nEntries);
			xAxis = 0.5*(m_EnergyBins[iEne+1]+m_EnergyBins[iEne]);
			dxAxis = 0.5*(m_EnergyBins[iEne+1]-m_EnergyBins[iEne]);
			//	xAxis = 0.5*(1/TMath::Sqrt(m_EnergyBins[iEne+1])+1/TMath::Sqrt(m_EnergyBins[iEne]));
			//dxAxis = 0.5*(1/TMath::Sqrt(m_EnergyBins[iEne])-1/TMath::Sqrt(m_EnergyBins[iEne+1]));
			
			yAxis = Reso[0]/(100.*xAxis);
			dyAxis = Reso[1];		
			axisTitle=";1/#sqrt{E_{truth}} [GeV^{-1/2}]; Resolution [%]";			
			
			g1->SetPoint(nPoint, xAxis, yAxis );
			g1->SetPointError(nPoint, dxAxis, dyAxis );
			++nPoint;
		}
		TString line = Form("%s : E [%4.0f; %4.0f]\tn=%3.0f\tx=%1.1f +/- %1.1f\ty=%1.3f +/- %1.3f\t%i", 
			variableName.Data(), m_EnergyBins[iEne], m_EnergyBins[iEne+1], nEntries, xAxis, dxAxis, yAxis, dyAxis, status );
		std::cout<< line.Data() <<std::endl;
	}
	g1->SetTitle(m_title + axisTitle );

 

TF1 *fitResolution = new TF1("func", "sqrt( [0]*[0]*x + [1]*[1] + [2]*[2]*x*x )", 20, 3000);
                        fitResolution->SetLineColor(2);
                        if (1==bin) fitResolution->SetLineColor(9);
                        fitResolution->SetLineWidth(2);
                        fitResolution->SetLineStyle(2);
			//fitResolution->SetParLimits(2, 1, 10);
			fitResolution->SetParameter(0, 0.50);
			fitResolution->SetParameter(1, 3);
			fitResolution->SetParameter(2, 5.);
                        g1->Fit(fitResolution, "q");
                        TString title = Form("a=%2.2f; b=%2.2f; c=%2.2f", fitResolution->GetParameter(0), fitResolution->GetParameter(1), fitResolution->GetParameter(2));
std::cout << "************ Fluctermfit**********" << std::endl;
			std::cout << title << std::endl;
			title = Form("Sigma(a)=%2.2f; Sigma(b)=%2.2f", fitResolution->GetParError(0), fitResolution->GetParError(1));
			std::cout << title << std::endl;

			std::cout << fitResolution->GetChisquare() << "/" << fitResolution->GetNDF() << std::endl;
			a=fitResolution->GetParameter(0);
			sampterm[0]=TMath::Abs(fitResolution->GetParameter(0)*100);
			sampterm[1]=TMath::Abs(fitResolution->GetParError(0)*100);
			sampterm[2]=TMath::Abs(fitResolution->GetParameter(1)*100);
			sampterm[3]=TMath::Abs(fitResolution->GetParError(1)*100);

	return a;
}
	

TGraphErrors* ReadHistogramms::plotFlutTermVsEta(std::string type, bool fit){
	TGraphErrors *g = new TGraphErrors();
	std::vector<double> resolution;

	Double_t yAxis=0, dyAxis=0;
	Double_t xAxis=0, dxAxis=0;
	TString axisTitle;
	Double_t nEntries = 0;			
	int status = 0;
	double Mean[2]={0}, Reso[2]={0};
	std::string variable;
	double Sampling[4] = {0};

	//	if ( m_nEnergyBins !=3 ) {
	//std::cout<<"\tImpossible to plot linearity vs eta for this collection of histograms !"<<std::endl;
	//return g;
	//}
	int iEne =0;
	// select eta bin
	//iEne=bin;
	//	else if (bin==1) iEne=2;
	//	else return g;	
	m_title = Form("Resolution fit" );

	//std::cout<<"\tLooking at bin "<< bin <<" : "<< m_title.Data() <<std::endl;
	// case 1 : linearity or resolution

	variable = "RatioEnergy_eta";
	
	std::cout<<"\n\tPlot "<< type <<" vs eta : "<< m_title.Data() <<std::endl;
	std::cout<< "****************** etabins " << m_nEtaBins<<std::endl;
	int iPoint=0;
	for (int iiEta=0; iiEta<m_nEtaBins; ++iiEta) {
	  std::cout<< "****************** eta " << iiEta <<std::endl;

	  double samplingterm=this->GetFluctTerm(type, iiEta,fit,Sampling);

	  std::cout << "sampling term "  << Sampling[0] <<std::endl;
		xAxis=0; 
		dxAxis=0;
		yAxis=0; 
		dyAxis=0;
							
		TString variableName = Form("%s_%s_eta%i_e%i", variable.c_str(), m_prefix.c_str(), iiEta, iEne );
		TH1F *histo = (TH1F*)m_file->Get( variableName.Data() );
		std::cout<< variableName.Data()  <<std::endl;
		nEntries = 3;			

		if ( nEntries>0 ) {
			
			xAxis  = 0.5*(m_EtaBins[iiEta+1]+m_EtaBins[iiEta]);
			dxAxis = 0.5*(m_EtaBins[iiEta+1]-m_EtaBins[iiEta]);
								
			if ( "sampling" == type ){
				status = fitGauss( histo, Mean, Reso, fit);
				yAxis = Sampling[0];
				dyAxis = Sampling[1];
				axisTitle=";#eta_{truth}; sampling term a [%]";			
			} 
			else if ( "noise" == type ){
				status = fitGauss( histo, Mean, Reso, fit);
				yAxis = TMath::Abs(Sampling[2]);
				dyAxis = TMath::Abs(Sampling[3]);
				if (dyAxis>yAxis) {yAxis=dyAxis;dyAxis=TMath::Abs(Sampling[2]);}
				axisTitle=";#eta_{truth}; constant term b [%]";			
			} 
			g->SetPoint(iPoint, xAxis, yAxis );
			g->SetPointError(iPoint, dxAxis, dyAxis );
			++iPoint;
		}
		TString line = Form("%s : eta [%1.2f; %1.2f]\tn=%3.0f\tx=%1.1f +/- %1.1f\ty=%1.3f +/- %1.3f\t%i", 
				    variableName.Data(), m_EtaBins[iiEta], m_EtaBins[iiEta+1], nEntries, xAxis, dxAxis, yAxis, dyAxis, status );
		std::cout<< line.Data() <<std::endl;
	}
	m_title = Form("Resolution fit" );
	g->SetTitle(m_title + axisTitle );
	return g;
}