/*
*	Last update: 14.04.2009 (minor changes)
*
*	TDjRes class: Djangoh radiative correction results container.
*
*	Usage:
*		TDjRes data;
*		data.Load("djangoh_r.dat");
*		Bool_t valueOK=data.GetValue(x,y,&corr,&error);
*	See also DjResExample_() functions at the end of this file.
*
*	See webpage for details: http://cern.ch/~makarenko/dj_res/
*/
#include <iostream.h>
#include <Rtypes.h>
#include <TTimer.h>
#include <TCanvas.h>
#include <TProfile.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TMath.h>
#include <TROOT.h> 
#include <TGaxis.h>
#include <TStyle.h>

// Max. number of points (1000 is default in Djangor user procedure)
#define NMAX_TDjRes		1000

class TDjRes {
public:
	Double_t fRPTSX[NMAX_TDjRes],fRPTSY[NMAX_TDjRes]; //data points arrays.
	Int_t fNPTSX,fNPTSY; //data table size
	Double_t fXPTMIN,fXPTMAX,fYPTMIN,fYPTMAX; //data table margins
	Int_t fNMODEX,fNMODEY; // data grid modes (0: linear, 1: logarythmic)
	Int_t fNSTATUS; // initialization flag (0:no data read, 1:ready, 2:error)
	Int_t fNCTRLSUM; // table control sum
	Double_t fRESDELTA[NMAX_TDjRes][NMAX_TDjRes]; // Precalculated delta values
	Double_t fRESERRSTA[NMAX_TDjRes][NMAX_TDjRes]; // Precalculated error
	Double_t fRESERRSYS[NMAX_TDjRes][NMAX_TDjRes]; // Precalculated error
	Double_t fS; // s-me^2-mp^2
	Double_t fDJ_RES_UNKNOWN; // to be returned as unknown result (if GetValue()==false).
	Double_t fDJ_ERR_UNKNOWN; // to be returned as unknown error (if GetValue()==false).
public:
	TDjRes() { fDJ_RES_UNKNOWN=-1.0e10; fDJ_ERR_UNKNOWN=1e+6; fNSTATUS=0; fS=0; fNCTRLSUM=0; }
	Bool_t Load(const char* fn);
	Bool_t GetValue(Double_t x, Double_t y, Double_t* res, Double_t* err); //err=Sqrt(errSta^2+errSys^2)
	Bool_t GetValue(Double_t x, Double_t y, Double_t* res, Double_t* errSta, Double_t* errSys);
	Bool_t IsReady();
	Bool_t GetLeftPoint(Double_t x, Double_t y, Int_t *nxleft, Int_t *nyleft);
};

Bool_t TDjRes::Load(const char* fn) {
	if ((fNSTATUS != 0) || (!fn)) return false;
	fNSTATUS = 2;
	cout << "Loading file: " << fn << "\n";
/*
c. ----- Results file format -----
c.   #(int) IVER: file version (reserved, currently 2)                 *
c.   #(int) I2ERRS: 1(0) the stat. & sys. errors are separated         *
c.   #(int) S: S-Mp2-Ml2                                               *
c.   #(int) NPTSX,NPTSY,NMODEX,NMODEY: number of points and grid mode  *
c.   #(double) XPTMIN,XPTMAX,YPTMIN,YPTMAX: min/max point values        *
c.   #(int) DELTA,ERRSTA[,ERRSYS] x(NPTSX*NPTSY): results array        *
c.   #(int) NCTRLSUM: currently NPTSX*NPTSY                            *
*/
	Int_t i1,i2,i3,iTmp1;
	Double_t dTmp1,dTmp2,dTmp3,dTmp4;
	Double_t valUnknwn=1e5; // 'unknown' value in dat-file presumed to be >=valUnknwn

	fS=-1;
	Int_t iver=0;
	ifstream ifs; ifs.open(fn); if(!ifs)return false;
	ifs>>iver; if(!ifs.good()){ifs.close();return false;}
	//if (iver==1) { cout<<"Can't load results file. Unexpected events file is found."; ifs.close(); return false; }
	if (iver!=2) { cout<<"Error loading results file. Unknown format."; ifs.close(); return false; }
	Bool_t bLoadErrSys=false;
	ifs>>iTmp1; if(!ifs.good()){cout<<"Error loading results file. Error flag expected.";ifs.close();return false;}
	if ((iTmp1!=0)&&(iTmp1!=1)){cout<<"Error loading results file. Unexpected error flag."<<iTmp1<<endl;ifs.close();return false;}
	if (iTmp1>0) bLoadErrSys=true;
	ifs>>fS; if(!ifs.good()){cout<<"Can't load S-parameter.\n";ifs.close();return false;}
	if (fS<=0) { cout<<"Can't load results file. S-parameter is not initialized."; ifs.close(); return false; }
	ifs>>fNPTSX>>fNPTSY>>fNMODEX>>fNMODEY; if(!ifs.good()){cout<<"Can't load results file\n";ifs.close();return false;} // Load data table size
	if ((fNPTSX>NMAX_TDjRes)||(fNPTSX>NMAX_TDjRes)) { cout<<"Can't allocate results file. Change NMAX_TDjRes parameter to "<<fNPTSX<<endl; ifs.close(); return false;  }
	if ((fNPTSX<=0)||(fNPTSY<=0)) { cout << "Invalid results file: unexpected negative points number"<<endl; ifs.close(); return false; }
	ifs>>dTmp1>>dTmp2>>dTmp3>>dTmp4; if(!ifs.good()){ifs.close();return false;}
	fXPTMIN=dTmp1; fXPTMAX=dTmp2; fYPTMIN=dTmp3; fYPTMAX=dTmp4;
	for(i1=0; i1<fNPTSX; i1++) {
		for(i2=0; i2<fNPTSY; i2++) {
			dTmp3=0.;
			ifs>>dTmp1>>dTmp2; if(bLoadErrSys)ifs>>dTmp3; if(!ifs.good()){cout<<"Can't load results file\n";ifs.close();return false;}
			if (dTmp1>=valUnknwn) { dTmp1=fDJ_RES_UNKNOWN; dTmp2=fDJ_ERR_UNKNOWN; dTmp3=fDJ_ERR_UNKNOWN; }
			fRESDELTA[i1][i2]=dTmp1;
			fRESERRSTA[i1][i2]=dTmp2;
			fRESERRSYS[i1][i2]=dTmp3;
		}
	}
	ifs >> fNCTRLSUM;
	ifs.close();
	// Check data control sum (currently number of numbers read).
	if (fNCTRLSUM != (fNPTSX*fNPTSY)) { cout << "Warning: Control sum failed\n"; return false; }
	if ((fNMODEX<0)||(fNMODEX>1)||(fNMODEY<0)||(fNMODEY>1)) { cout << "Invalid results file: unexpected grid mode\n"; return false; }
	if ((fXPTMIN>fXPTMAX)||(fYPTMIN>fYPTMAX)) { cout<<"Invalid results file: unexpected grid margins"; return false; }
	if (((fXPTMIN<=0.)||(fXPTMAX<=0.))&&(fNMODEX==1)) { cout<<"Invalid results file: NEGATIVE X MARGINS + LOG-SCALE"; return false; }
	if (((fYPTMIN<=0.)||(fYPTMAX<=0.))&&(fNMODEY==1)) { cout<<"Invalid results file: NEGATIVE Y MARGINS + LOG-SCALE"; return false; }
	if ((fNPTSX<2)||(fNPTSY<2)) { cout<<"Error loading results file. Invalid fNPTSX or fNPTSY"; return false; }

	// Initialize points
	if (fNMODEX==1) { dTmp1=TMath::Log((Double_t)fXPTMIN); dTmp2=TMath::Log((Double_t)fXPTMAX); }
	for (i1=0;i1<fNPTSX;i1++) {
		if (fNMODEX==0) { // linear
			fRPTSX[i1]=fXPTMIN+(fXPTMAX-fXPTMIN)*i1/(fNPTSX-1);
		} else { // log
			fRPTSX[i1]=TMath::Exp(dTmp1+(dTmp2-dTmp1)*i1/(fNPTSX-1));
		}
	}
	if (fNMODEY==1) { dTmp1=TMath::Log((Double_t)fYPTMIN); dTmp2=TMath::Log((Double_t)fYPTMAX); }
	for (i1=0;i1<fNPTSY;i1++) {
		if (fNMODEY==0) { // linear
			fRPTSY[i1]=fYPTMIN+(fYPTMAX-fYPTMIN)*i1/(fNPTSY-1);
		} else { // log
			fRPTSY[i1]=TMath::Exp(dTmp1+(dTmp2-dTmp1)*i1/(fNPTSY-1));
		}
	}
	fNSTATUS=1;
	cout << "Results file " << fn << " loaded succesfully\n" << "fNPTSX="<<fNPTSX<<", fNPTSY="<<fNPTSY<<"\n";
	cout << "fXPTMIN="<<fXPTMIN<<", fXPTMAX="<<fXPTMAX<<"\n" << "fYPTMIN="<<fYPTMIN<<", fYPTMAX="<<fYPTMAX<<"\n";
	cout << "fNMODEX="<<fNMODEX<<", fNMODEY="<<fNMODEY<<"\n" << "fNCTRLSUM="<<fNCTRLSUM<<"\n";
	cout << "S="<<fS<<", Sys. error loaded: "<<bLoadErrSys<<"\n";
	return true;
}

Bool_t TDjRes::GetValue(Double_t x, Double_t y, Double_t* res, Double_t* err) {
	Double_t errSta=0.; Double_t errSys=fDJ_ERR_UNKNOWN;
	Bool_t bRes=GetValue(x, y, res, &errSta, &errSys);
	if(err) {
		if ((errSta<fDJ_ERR_UNKNOWN) && (errSys<fDJ_ERR_UNKNOWN)) {
			*err=TMath::Sqrt(errSta*errSta+errSys*errSys);
		} else {
			*err=fDJ_ERR_UNKNOWN;
		}
	}
	return bRes;
}

Bool_t TDjRes::GetValue(Double_t x, Double_t y, Double_t* res, Double_t* errSta, Double_t* errSys) {
	if(res)*res=fDJ_RES_UNKNOWN; if(errSta)*errSta=0.; if(errSys)*errSys=fDJ_ERR_UNKNOWN;
	if (fNSTATUS != 1) return false;
	Int_t nxleft,nyleft;
	
	if(!GetLeftPoint(x,y,&nxleft,&nyleft)) return false;
	if ((nxleft==(fNPTSX-1))&&(nxleft>0))nxleft=nxleft-1; // The right point (nxleft+1) must be valid.
	if ((nyleft==(fNPTSY-1))&&(nyleft>0))nyleft=nyleft-1; // The right point (nyleft+1) must be valid.
	if ((nxleft<0)||(nxleft>fNPTSX)||(nyleft<0)||(nyleft>fNPTSY)) return false;

	// Get values of the four nearest points
	Double_t valLL,valRL,valLR,valRR,staLL,staRL,staLR,staRR,sysLL,sysRL,sysLR,sysRR;
	valLL=fRESDELTA[nxleft][nyleft]; valRL=fRESDELTA[nxleft+1][nyleft];
	valLR=fRESDELTA[nxleft][nyleft+1]; valRR=fRESDELTA[nxleft+1][nyleft+1];
	staLL=fRESERRSTA[nxleft][nyleft]; staRL=fRESERRSTA[nxleft+1][nyleft];
	staLR=fRESERRSTA[nxleft][nyleft+1]; staRR=fRESERRSTA[nxleft+1][nyleft+1];
	sysLL=fRESERRSYS[nxleft][nyleft]; sysRL=fRESERRSYS[nxleft+1][nyleft];
	sysLR=fRESERRSYS[nxleft][nyleft+1]; sysRR=fRESERRSYS[nxleft+1][nyleft+1];

	// Check points for 'unknown'
	if ((valLL==fDJ_RES_UNKNOWN)&&(valLR==fDJ_RES_UNKNOWN)&&(valRL==fDJ_RES_UNKNOWN)&&(valRR==fDJ_RES_UNKNOWN)) return false;
	if (valLL==fDJ_RES_UNKNOWN) {
		if(valLR!=fDJ_RES_UNKNOWN){valLL=valLR; } else if(valRL!=fDJ_RES_UNKNOWN){valLL=valRL; } else if(valRR!=fDJ_RES_UNKNOWN){valLL=valRR; }
		staLL=0.; sysLL=1.;
	}
	if (valLR==fDJ_RES_UNKNOWN) {
		if(valLL!=fDJ_RES_UNKNOWN){valLR=valLL;} else if(valRR!=fDJ_RES_UNKNOWN){valLR=valRR;} else if(valRL!=fDJ_RES_UNKNOWN){valLR=valRL;}
		staLR=0.; sysLR=1.;
	}
	if (valRL==fDJ_RES_UNKNOWN) {
		if(valLL!=fDJ_RES_UNKNOWN){valRL=valLL;} else if(valRR!=fDJ_RES_UNKNOWN){valRL=valRR;} else if(valLR!=fDJ_RES_UNKNOWN){valRL=valLR;}
		staRL=0.; sysRL=1.;
	}
	if (valRR==fDJ_RES_UNKNOWN) {
		if(valLR!=fDJ_RES_UNKNOWN){valRR=valLR;} else if(valRL!=fDJ_RES_UNKNOWN){valRR=valRL;} else if(valLL!=fDJ_RES_UNKNOWN){valRR=valLL;}
		staRR=0.; sysRR=1.;
	}

	// Average over 4 points
	Double_t dx=fRPTSX[nxleft+1]-fRPTSX[nxleft]; Double_t dxl=x-fRPTSX[nxleft]; Double_t dxr=fRPTSX[nxleft+1]-x;
	Double_t dy=fRPTSY[nyleft+1]-fRPTSY[nyleft]; Double_t dyl=y-fRPTSY[nyleft]; Double_t dyr=fRPTSY[nyleft+1]-y;
	Double_t valAL,valAR,valAA,staAL,staAR,staAA,sysAL,sysAR,sysAA;
	if (dx > 1e-14) {
		valAL=(valLL*dxr+valRL*dxl)/dx; staAL=(staLL*dxr+staRL*dxl)/dx; sysAL=(sysLL*dxr+sysRL*dxl)/dx;
		valAR=(valLR*dxr+valRR*dxl)/dx; staAR=(staLR*dxr+staRR*dxl)/dx; sysAR=(sysLR*dxr+sysRR*dxl)/dx;
	} else {
		valAL=(valLL+valRL)/2; staAL=(staLL+staRL)/2; sysAL=(sysLL+sysRL)/2;
		valAR=(valLR+valRR)/2; staAR=(staLR+staRR)/2; sysAR=(sysLR+sysRR)/2;
	}
	if (dy > 1e-14) {
		valAA=(valAL*dyr+valAR*dyl)/dy; staAA=(staAL*dyr+staAR*dyl)/dy; sysAA=(sysAL*dyr+sysAR*dyl)/dy;
	} else {
		valAA=(valAL+valAR)/2; staAA=(staAL+staAR)/2; sysAA=(sysAL+sysAR)/2;
	}

	if (res) *res=valAA;
	if (errSta) *errSta=staAA;
	if (errSys) *errSys=sysAA;
	return true;
}

Bool_t TDjRes::IsReady() { return (fNSTATUS==1); }
Bool_t TDjRes::GetLeftPoint(Double_t x, Double_t y, Int_t *nxleft, Int_t *nyleft) {
	if ((x<fXPTMIN)||(x>fXPTMAX)||(y<fYPTMIN)||(y>fYPTMAX)) return false;
	// Find left X-point (nxleft)
	Int_t i1=0; Int_t i2=fNPTSX-1;
	Int_t nxlTmp, nylTmp;
	for (;;) {
		nxlTmp=i1+(i2-i1)/2;
		if(fRPTSX[nxlTmp]<x) {
			i1=nxlTmp; if(i2>(i1+1))continue;
			if(fRPTSX[i2]==x)nxlTmp=i2; // check for x==fXPTMAX case only
		}
		else if (fRPTSX[nxlTmp]>x) { i2=nxlTmp; if(i2>(i1+1))continue; nxlTmp=i1; }
		break;
	}
	//Find left Y-point (nyleft)
	i1=0; i2=fNPTSY-1;
	for (;;) {
		nylTmp=i1+(i2-i1)/2;
		if (fRPTSY[nylTmp]<y) {
			i1=nylTmp; if(i2>(i1+1))continue;
			if(fRPTSY[i2]==y)nylTmp=i2; // check for y==fYPTMAX case only
		}
		else if (fRPTSY[nylTmp]>y) { i2=nylTmp; if(i2>(i1+1))continue; nylTmp=i1; }
		break;
	}
	if (nxleft) *nxleft=nxlTmp;
	if (nyleft) *nyleft=nylTmp;
	return true;
}

/******************************************************************************
***                                                                         ***
***                               Examples                                  ***
***                                                                         ***
******************************************************************************/
#define HIST_nNx		100
#define HIST_nNy		100
#define HIST_nNq2		100

/***********************************************************
*	Example #1
*	Plot and compare 1D histograms for different ".dat"-files.
***********************************************************/
TH1F* DjResExNewYPlot(const char* id, TDjRes* data, Double_t x, Bool_t bLogY, Double_t* aErrs, Double_t DeltaMax, Bool_t bShift);
Int_t DjResExampleYPlots() {
	cout<<"DjResExampleYPlots()"<<endl; 
	TDjRes data1; if(!data1.Load("djangoh_r_1.dat")) { cout<<"Warning: Can't open data file #1"<<endl; }
	TDjRes data2; if(!data2.Load("djangoh_r_2.dat")) { cout<<"Warning: Can't open data file #2"<<endl; }
	TDjRes data3; if(!data3.Load("djangoh_r_3.dat")) { cout<<"Warning: Can't open data file #3"<<endl; }

	// shift one of the plots for better combining
	Bool_t shift1=false; Bool_t shift2=true; Bool_t shift3=false;

	Double_t DeltaMax=0.2; // cut of huge values

	Int_t i;
	TCanvas *cl=new TCanvas("cl","Example",10,10,1000,1400);
	TH1F *hist[3][4];
	Double_t x[]={0.01,0.03,0.1,0.4};

	Double_t err11[HIST_nNy],err12[HIST_nNy],err13[HIST_nNy],err14[HIST_nNy]; Double_t *err1[]={err11,err12,err13,err14};
	if (data1.IsReady()) {
		// Note: DjResExNewYPlot() function is implemented below
		for(i=0;i<4;i++){ hist[0][i]=DjResExNewYPlot("1",&data1,x[i],false,err1[i],DeltaMax,shift1); hist[0][i]->SetLineColor(kRed); }
	}
	Double_t err21[HIST_nNy],err22[HIST_nNy],err23[HIST_nNy],err24[HIST_nNy]; Double_t *err2[]={err21,err22,err23,err24};
	if (data2.IsReady()) {
		for(i=0;i<4;i++){ hist[1][i]=DjResExNewYPlot("2",&data2,x[i],false,err2[0],DeltaMax,shift2);hist[1][i]->SetLineColor(kGreen); }
	}
	Double_t err31[HIST_nNy],err32[HIST_nNy],err33[HIST_nNy],err34[HIST_nNy]; Double_t *err3[]={err31,err32,err33,err34};
	if (data3.IsReady()) {
		for(i=0;i<4;i++){ hist[2][i]=DjResExNewYPlot("3",&data3,x[i],false,err3[i],DeltaMax,shift3);hist[2][i]->SetLineColor(kBlue); }
	}

	cl->Divide(2,2);
	for(i=0;i<4;i++){
		cl->cd(i+1); Bool_t bFirst=true;
		if (data1.IsReady()) { hist[0][i]->Draw((bFirst?"":"same")); bFirst=false; }
		if (data2.IsReady()) { hist[1][i]->Draw((bFirst?"":"same")); bFirst=false; }
		if (data3.IsReady()) { hist[2][i]->Draw((bFirst?"":"same")); bFirst=false; }
	}
	cl->SaveAs("DjResExampleYPlots.eps");
	cout<<"DjResExampleYPlots() finished."<<endl;
	return 1;
}

/*
*	DjResExNewYPlot: create 1D histogram at fixed x.
*	Used by Example #1 (see above).
*	'bShift' flag: shift y-margins by 0.5 of the bin width
*	(used for better view with 'unshifted' plot on single canvas).
*/
TH1F* DjResExNewYPlot(const char* id, TDjRes* data, Double_t x, Bool_t bLogY, Double_t* aErrs, Double_t DeltaMax, Bool_t bShift) {
	if (data==NULL) { cout<<"Error: data==NULL"; return NULL; }
	if (!data->IsReady()) return NULL;

	Int_t i,j;
	Double_t y,lnMin,lnMax;
	char buf[100],buf2[100];

	Int_t Ny=HIST_nNy;
	Double_t yMin=data->fYPTMIN; Double_t yMax=data->fYPTMAX;

	// shift {yMin,yMax} for half of the bin
	if (bShift) {
		if(bLogY){
			lnMin=TMath::Log((Double_t)yMin); lnMax=TMath::Log((Double_t)yMax);
			Double_t tmp1=TMath::Exp(lnMin+(lnMax-lnMin)*0.5/Ny);
			Double_t tmp2=TMath::Exp(lnMin+(lnMax-lnMin)*(0.5+Ny)/Ny);
			yMin=tmp1; yMax=tmp2;
		} else {
			Double_t tmp1=yMin+(yMax-yMin)*0.5/(Ny-1);
			Double_t tmp2=yMin+(yMax-yMin)*(-0.5+Ny)/(Ny-1);
			yMin=tmp1; yMax=tmp2;
		}
	}

	// log-scale init
	Double_t binslogY[HIST_nNy+1];
	if(bLogY) { lnMin=TMath::Log((Double_t)yMin); lnMax=TMath::Log((Double_t)yMax); for(i=0;i<=Ny;i++){binslogY[i]=TMath::Exp(lnMin+(lnMax-lnMin)*i/Ny);} }

	TH1F* hist=NULL;
	sprintf(buf,"plot %s, x=%.3f",id,x); sprintf(buf2,"delta, x=%.2f",x);
	if(bLogY){ hist=new TH1F(buf,buf2,Ny,binslogY); } else { hist=new TH1F(buf,buf2,Ny,yMin,yMax); }
	hist->GetXaxis()->SetTitle("y"); hist->GetYaxis()->SetTitle("delta"); hist->SetStats(false);

	for(i=0;i<Ny;i++) {
		// get middle of y-bin
		if(bLogY){ y=(binslogY[i]+binslogY[i+1])/2; } else { y=yMin+(yMax-yMin)*(0.5+i)/Ny; }
		Double_t res=0; Double_t err=1.; Double_t errSta=0.; Double_t errSys=1.;
		Bool_t bOk=data->GetValue(x,y,&res,&errSta,&errSys);
		err=TMath::Sqrt(errSta*errSta+errSys*errSys); //err=errSta; //err=errSys;
		if (bOk) { aErrs[i]=err; hist->Fill(y,res); }
		//if((err>parEMax)||(!bOk))err=parEMax; /*if(bOk)*/aErrs[i]=err;
	}
	hist->SetError(aErrs);

	if (DeltaMax>0) { hist->SetMinimum(-DeltaMax); hist->SetMaximum(DeltaMax); }
	return hist;
}

/***********************************************************
*	Example #2
*	Plot 2D log-axis histograms {x,Q2} for delta and errors.
***********************************************************/
Int_t DjResExample2D() {
	TDjRes data; if(!data.Load("djangoh_r.dat"))return 0;

	TCanvas *cl=new TCanvas("cl","Example",10,10,1000,700);

	Int_t i,j,k;
	Double_t x,y,Q2;
	Int_t Nx=HIST_nNx; Int_t Ny=HIST_nNy; Int_t Nq=HIST_nNq2;

	// limits
	Double_t xMin=data.fXPTMIN; Double_t xMax=data.fXPTMAX;
	Double_t yMin=data.fYPTMIN; Double_t yMax=data.fYPTMAX;
	Double_t Q2min=100.0; Double_t Q2max=data.fS;

	// log axis points (binslogY is actually unused here)
	Double_t binslogX[HIST_nNx+1]; Double_t binslogY[HIST_nNy+1]; Double_t binslogQ2[HIST_nNq2+1];
	Double_t lnMin=TMath::Log((Double_t)xMin); Double_t lnMax=TMath::Log((Double_t)xMax); for(i=0;i<=Nx;i++){binslogX[i]=TMath::Exp(lnMin+(lnMax-lnMin)*i/Nx);}
	lnMin=TMath::Log((Double_t)yMin); lnMax=TMath::Log((Double_t)yMax); for(i=0;i<=Ny;i++){binslogY[i]=TMath::Exp(lnMin+(lnMax-lnMin)*i/Ny);}
	lnMin=TMath::Log((Double_t)Q2min); lnMax=TMath::Log((Double_t)Q2max); for(i=0;i<=Nq;i++){binslogQ2[i]=TMath::Exp(lnMin+(lnMax-lnMin)*i/Nq);}

	// book histos
	TH2F* h2_DEL_xQ2=new TH2F("delta","delta",Nx,binslogX,Nq,binslogQ2); h2_DEL_xQ2->GetXaxis()->SetTitle("x"); h2_DEL_xQ2->GetYaxis()->SetTitle("Q2"); h2_DEL_xQ2->SetStats(false);
	TH2F* h2_ERR_xQ2=new TH2F("error","error",Nx,binslogX,Nq,binslogQ2); h2_ERR_xQ2->GetXaxis()->SetTitle("x"); h2_ERR_xQ2->GetYaxis()->SetTitle("Q2"); h2_ERR_xQ2->SetStats(false);
	TH2F* h2_SYS_xQ2=new TH2F("err. sys.","err. sys.",Nx,binslogX,Nq,binslogQ2); h2_SYS_xQ2->GetXaxis()->SetTitle("x"); h2_SYS_xQ2->GetYaxis()->SetTitle("Q2"); h2_SYS_xQ2->SetStats(false);
	TH2F* h2_STA_xQ2=new TH2F("err. stat. ","err. stat",Nx,binslogX,Nq,binslogQ2); h2_STA_xQ2->GetXaxis()->SetTitle("x"); h2_STA_xQ2->GetYaxis()->SetTitle("Q2"); h2_STA_xQ2->SetStats(false);

	// fill {x,Q2} points
	for(i=0;i<Nx;i++) {
		x=(binslogX[i]+binslogX[i+1])/2; // get middle of the bin
		for(j=0;j<Nq;j++) {
			Q2=(binslogQ2[j]+binslogQ2[j+1])/2; y=Q2/data.fS/x;
			if(y>1)continue;
			Double_t res=0; Double_t errSta=0; Double_t errSys=0;
			if (data.GetValue(x,y,&res,&errSta,&errSys)) {
				if (res!=data.fDJ_RES_UNKNOWN) {
					h2_DEL_xQ2->Fill(x,Q2,res);
					h2_STA_xQ2->Fill(x,Q2,errSta);
					h2_SYS_xQ2->Fill(x,Q2,errSys);
					h2_ERR_xQ2->Fill(x,Q2,TMath::Sqrt(errSta*errSta+errSys*errSys));
				}
			}
		}
	}
	// set hard limits to cut unphysical results on edges
	h2_DEL_xQ2->SetMinimum(-0.25);h2_DEL_xQ2->SetMaximum(0.25);
	h2_STA_xQ2->SetMinimum(0.);h2_STA_xQ2->SetMaximum(0.2); // truncate error of more than 10% (to total c/s)
	h2_SYS_xQ2->SetMinimum(0.);h2_SYS_xQ2->SetMaximum(0.2); // truncate error of more than 10% (to total c/s)
	h2_ERR_xQ2->SetMinimum(0.);h2_ERR_xQ2->SetMaximum(0.2); // truncate error of more than 10% (to total c/s)
	
	// set palette (grayscale)
	Int_t ncolors=18; Int_t clr[18]; Bool_t bClrsOK=true;
	for(i=0; i<ncolors; i++) {
		Float_t r=1.*(ncolors-i)/(ncolors-1); Float_t g=r; Float_t b=r;
		TColor* color=gROOT->GetColor(210+i); if (color!=NULL) { color->SetRGB(r,g,b); clr[i]=210+i; } else { bClrsOK=false; break; }
	}
	if(bClrsOK)gStyle->SetPalette(ncolors,clr);

	cl->Divide(2,2);
	cl->cd(1); gPad->SetLogx(); gPad->SetLogy(); h2_DEL_xQ2->Draw("colz");
	cl->cd(2); gPad->SetLogx(); gPad->SetLogy(); h2_ERR_xQ2->Draw("colz");
	cl->cd(3); gPad->SetLogx(); gPad->SetLogy(); h2_STA_xQ2->Draw("colz");
	cl->cd(4); gPad->SetLogx(); gPad->SetLogy(); h2_SYS_xQ2->Draw("colz");
	cl->cd();
	cl->SaveAs("DjResExample2D.eps");
	return true;
}