/*
*	Last update: 14.04.2009 (minor changes)
*
*	TDjResEvt class: Djangoh EVENTS file container.
*	This class is used for DEBUG PURPOSES ONLY.
*	Use TDjRes class as Djangoh radiative correction results container.
*	See webpage for details: http://cern.ch/~makarenko/dj_res/
*
*	Files:
*	1. "djangoh_r.dat": precalculared cross section ratio table.
*		Use TDjRes class to operate with it.
*		Created either by Djangoh User Procedure (or using TDjResEvt class).
*	2. "djangoh_r.evts": tables of events generated at Djangoh sampling stage.
*		Used by TDjResEvt class for DEBUG purposes only (e.g. to regenerate "djangoh_r.dat" file).
*		Created by Djangor User Procedure.
*		The events nubmer is counted in sertain bins around the set of {x,y} points.
*		Usually contains event tables for various variables sets.
*		The leptonic- and JB-variables data may be used for internal checks.
*		Refer to Djangor User Procedure for 'additional tables' details & order.
*		It also can contain Born events (sampled separatly).
*		Djangoh User Procedure automatically combines RC and Born sampling data into single file if possible.
*	3. "djangoh_r.born": precalculared born cross section table.
*		Used by TDjResEvt class for DEBUG purposes only (e.g. to regenerate "djangoh_r.dat" file).
*		Created by Djangor User Procedure after every sampling procedure.
*	See webpage for details: http://cern.ch/~makarenko/dj_res/
*
*	TDjResEvt usage:
*	1. Get delta value from events file:
*		TDjResEvt data;
*		data.Load("djangoh_r.evts","djangoh_r.born");
*		Bool_t valueOK=data.GetValue(x,y,&corr,&error);
*	The value is calculated as 'number-of-events'*'c/s per event'/{born c/s}-1.
*	The Born value is taken from "djangoh_r.born" (preferred) or from Born-events sampling data.
*
*	2. Update the results file ("djangoh_r.dat") by current events data
*	Note: this task is usually performed in DJANGOH user procedure.
*		data.LoadDjRes("djangoh_r.dat");
*		data.UpdateDjRes();
*		data.SaveDjRes("djangoh_r.dat");
*
*	3. Get RC value in other variables (e.g. from events pane #3(JB-variables))
*		valueOK=data.GetValue(x,y,&corr,&error,3);
*
*	See also DjResEvtExample_() functions at the end of this file.
*/
#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 (200 is default in Djangor user procedure)
#define NMAX_TDjResEvt		1000
// Max. number of additional data tables
#define MAXADDS_TDjResEvt	15
// Max. number of parameters in "djangoh_r.evts"-file
#define MAXDPARS_TDjResEvt	20

class TDjResEvt {
public:
	Double_t fRPTSX[NMAX_TDjResEvt],fRPTSY[NMAX_TDjResEvt]; //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
	Int_t fNEVTOT_BORN,fNEVTOT_CORR; // total number of events generated (usually LESS then sum of events in table)
	Int_t fNEVTS_BORN[NMAX_TDjResEvt][NMAX_TDjResEvt],fNEVTS_CORR[NMAX_TDjResEvt][NMAX_TDjResEvt]; // number of events tables
	// Delta & Error values (operated by LoadDjRes,UpdateDjRes,SaveDjRes functions)
	Bool_t fRESArrayFilled;
	Double_t fRESDELTA[NMAX_TDjResEvt][NMAX_TDjResEvt];
	Double_t fRESERRSTA[NMAX_TDjResEvt][NMAX_TDjResEvt];
	Double_t fRESERRSYS[NMAX_TDjResEvt][NMAX_TDjResEvt];
	// Additional data tables (correction in other variables, see dhangoh_u_r.f code for details)
	Int_t fNEVTS_AddN;
	Int_t fNEVTS_Add[NMAX_TDjResEvt][NMAX_TDjResEvt][MAXADDS_TDjResEvt];
	// DJANGOH parameters
	Int_t fNDPARS;
	Double_t fDPARS[MAXDPARS_TDjResEvt];
	Double_t fFBINWX,fFBINWY; // width of event counting region at sampling, x(1-FBINWX)<x<x(1+FBINWX)
	Int_t fCutMode;
	Double_t fXmin,fXmax,fYmin,fYmax,fQ2min,fQ2max,fW2min; // cuts
	Double_t fS;
	Double_t fCS_Evt, fCS_EvtBorn; // cross section per event
	// Precalculated Born
	Bool_t fBornOK;
	Double_t fBorn[NMAX_TDjResEvt][NMAX_TDjResEvt];
	// 'unknown values'
	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:
	TDjResEvt() { fNSTATUS=0; fNEVTS_AddN=0; fBornOK=false; fDJ_RES_UNKNOWN=1e6; fDJ_ERR_UNKNOWN=1.; fNDPARS=0; fFBINWX=0.1; fFBINWY=0.1; fRESArrayFilled=0; fCS_Evt=-1; fCS_EvtBorn=-1; }
	Bool_t Load(const char* fn, const char* fn_born=NULL);
	Bool_t IsReady(); // returns  (fNSTATUS==1)
	Bool_t GetValue(Double_t x, Double_t y, Double_t* res, Double_t* err, Int_t nAddTbl=-1);
	Bool_t GetValue(Double_t x, Double_t y, Double_t* res, Double_t* errSta, Double_t* errSys, Int_t nAddTbl);
	Bool_t GetValueI(Int_t nx, Int_t ny, Double_t* res, Double_t* errSta, Double_t* errSys, Int_t nAddTbl);
	Double_t GetBornCS(Double_t x, Double_t y);
	Bool_t CalcDelta(Int_t nx, Int_t ny, Double_t* delta, Int_t nAddTbl=-1);
	Bool_t CalcError(Int_t nx, Int_t ny, Double_t* errSta, Double_t* errSys, Double_t delta, Int_t nAddTbl=-1);
	Int_t GetEvNumBorn(Double_t x, Double_t y);
	Int_t GetEvNum(Double_t x, Double_t y, Int_t nAddTbl=-1);
	Bool_t GetLeftPoint(Double_t x, Double_t y, Int_t *nxleft, Int_t *nyleft, Bool_t bNeedRightOK);
	Bool_t LoadDjRes(const char* fn);
	Bool_t UpdateDjRes(Int_t nAddTbl=-1); //nAddTbl=-1 means correction in ZEUS variables
	Bool_t SaveDjRes(const char* fn);
};

Bool_t TDjResEvt::Load(const char* fn, const char* fn_born) {
	if ((fNSTATUS != 0) || (!fn)) return false;
	fNSTATUS = 2;
	cout << "Loading file: " << fn << "\n";
/*
C. ----- Data file format ----- ver.1
C. (int) IVER :file format version (currently 1, reserved)
C. (int) fNDPARS :size of parameters array
C. (double) fDPARS(NPARAMS) : array of parameters (cuts and other DJANGOH settings)
C. #0-12: FBINWX,FBINWY,ICUT,XMIN,XMAX,YMIN,YMAX,Q2MIN,Q2MAX,WMIN,GSP,CS_EVT,CS_EVT_BORN
C. (int) NPTSX,NPTSY,NMODEX,NMODEY
C. (double) XPTMIN,XPTMAX,YPTMIN,YPTMAX
C. (int) NEVTS_BORN,IEVTS_BORN(NPTSX*NPTSY)
C. (int) NEVTS_CORR,IEVTS_CORR(NPTSX*NPTSY)
C. (int) fNEVTS_AddN (number of additional data tables)
C. (int) fNEVTS_Add(NPTSX*NPTSY*fNEVTS_AddN) (additional data)
C. (int) NCTRLSUM: currently NPTSX*NPTSY
*/
	Int_t i1,i2,i3,iTmp1,iTmp2,iTmp3,iTmp4;
	Double_t dTmp1,dTmp2,dTmp3,dTmp4;
	fFBINWX=-1; fFBINWY=-1;
	fCutMode=-1; fXmin=-1; fXmax=-1; fYmin=-1; fYmax=-1; fQ2min=-1; fQ2max=-1; fW2min=-1;
	fS=-1.; fCS_Evt=-1.; fCS_EvtBorn=-1,; fBornOK=false;

	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<<"Error loading events file. Unknown format."; ifs.close(); return false; }
	ifs>>fNDPARS; if(!ifs.good()){ifs.close();return false;} // Load DJANGOH parameters
	if (fNDPARS<0) { cout<<"Invalid number of DJANGOH parameters: "<<fNDPARS<<endl; ifs.close(); return false;  }
	if (fNDPARS>MAXDPARS_TDjResEvt) { cout<<"Too many DJANGOH parameters. Change MAXDPARS_TDjResEvt parameter to "<<fNDPARS<<", current: "<<MAXDPARS_TDjResEvt<<endl; ifs.close(); return false;  }
	for(i1=0; i1<fNDPARS; i1++) {
		ifs>>dTmp1; if(!ifs.good()){ifs.close();return false;}
		fDPARS[i1]=dTmp1;
		if ((i1==0) && (dTmp1>0) && (dTmp1<1)) fFBINWX=dTmp1;
		if ((i1==1) && (dTmp1>0) && (dTmp1<1)) fFBINWY=dTmp1;
		if ((i1==2) && (dTmp1>0) && (dTmp1<4)) fCutMode=(dTmp1+0.01); // 0.1 gap for precision error
		if ((i1==3) && (dTmp1>=0) && (dTmp1<=1)) fXmin=dTmp1;
		if ((i1==4) && (dTmp1>=0) && (dTmp1<=1)) fXmax=dTmp1;
		if ((i1==5) && (dTmp1>=0) && (dTmp1<=1)) fYmin=dTmp1;
		if ((i1==6) && (dTmp1>=0) && (dTmp1<=1)) fYmax=dTmp1;
		if ((i1==7) && (dTmp1>=0)) fQ2min=dTmp1;
		if ((i1==8) && (dTmp1>=fQ2min)) fQ2max=dTmp1;
		if ((i1==9) && (dTmp1>=0)) fW2min=dTmp1;
		if ((i1==10) && (dTmp1>=0)) fS=dTmp1;
		if ((i1==11) && (dTmp1>=0)) fCS_Evt=dTmp1;
		if ((i1==12) && (dTmp1>=0)) fCS_EvtBorn=dTmp1;
	}
	ifs>>fNPTSX>>fNPTSY>>fNMODEX>>fNMODEY; if(!ifs.good()){ifs.close();return false;} // Load data table size
	if ((fNPTSX>NMAX_TDjResEvt)||(fNPTSX>NMAX_TDjResEvt)) { cout<<"Can't allocate events file. Change NMAX_TDjResEvt parameter to "<<fNPTSX<<endl; ifs.close(); return false;  }
	if ((fNPTSX<=0)||(fNPTSY<=0)) { cout << "Invalid events 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;
	ifs >> fNEVTOT_BORN; if(!ifs.good()){ifs.close();return false;}
	for(i1=0; i1<fNPTSX; i1++) {
		for(i2=0; i2<fNPTSY; i2++) {
			ifs>>iTmp1; if(!ifs.good()){ifs.close();return false;}
			fNEVTS_BORN[i1][i2]=iTmp1;
		}
	}
	ifs >> fNEVTOT_CORR; if(!ifs.good()){ifs.close();return false;}
	for(i1=0; i1<fNPTSX; i1++) {
		for (i2=0; i2<fNPTSY; i2++) {
			ifs>>iTmp1; if(!ifs.good()){ifs.close();return false;}
			fNEVTS_CORR[i1][i2]=iTmp1;
		}
	}
	// load additional (debug) tables if exist
	ifs >> fNEVTS_AddN; if(!ifs.good()){ifs.close();return false;}
	if(fNEVTS_AddN<0)fNEVTS_AddN=0;
	Int_t NAddInFile=fNEVTS_AddN;
	if (fNEVTS_AddN>MAXADDS_TDjResEvt) { fNEVTS_AddN=MAXADDS_TDjResEvt; cout<<"Can't allocate "<<NAddInFile<<" additional tables. The only "<<fNEVTS_AddN<<" will be loaded"<<endl; }
	for(i3=0; i3<NAddInFile; i3++) {
		for(i1=0; i1<fNPTSX; i1++) {
			for (i2=0; i2<fNPTSY; i2++) {
				ifs>>iTmp1; if(!ifs.good()){ifs.close();return false;}
				if (i3<fNEVTS_AddN) fNEVTS_Add[i1][i2][i3]=iTmp1;
			}
		}
	}
	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 events file: unexpected grid mode\n"; return false; }
	if ((fXPTMIN>fXPTMAX)||(fYPTMIN>fYPTMAX)) { cout<<"Invalid events file: unexpected grid margins"; return false; }
	if (((fXPTMIN<=0.)||(fXPTMAX<=0.))&&(fNMODEX==1)) { cout<<"Invalid events file: NEGATIVE X MARGINS + LOG-SCALE"; return false; }
	if (((fYPTMIN<=0.)||(fYPTMAX<=0.))&&(fNMODEY==1)) { cout<<"Invalid events file: NEGATIVE Y MARGINS + LOG-SCALE"; return false; }
	if ((fNPTSX<2)||(fNPTSY<2)) { cout<<"Invalid fNPTSX or fNPTSY"; return false; }
	if (fNEVTOT_CORR==0) { cout<<"Error: 0 events in table. fNEVTOT_CORR:" << fNEVTOT_CORR << endl; 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));
		}
	}

// ".born" file format
//  #(int) IVER: file id/version (reserved, currently 1)
//  #(int) S: S-Mp2-Ml2
//  #(int) NPTSX,NPTSY,NMODEX,NMODEY
//	#(double) XPTMIN,XPTMAX,YPTMIN,YPTMAX
//  #(double) DSIGMA/DXDY x(NPTSX*NPTSY) in [pb]
//  #(int) NCTRLSUM: currently NPTSX*NPTSY
	if (fn_born!=NULL) {
		ifstream ifsB; ifsB.open(fn_born); Bool_t BrdOK=true;
		if (ifsB) {
			if(BrdOK) {
				Int_t iverB=0; ifsB>>iverB; if(!ifsB.good()){ BrdOK=false; }
				if(iverB!=1) { cout<<"Error loading Born data file. Unknown format."; BrdOK=false; }
			}
			if(BrdOK) {
				ifsB>>dTmp1;
				if(!ifsB.good()){ BrdOK=false; }
				else { if ((dTmp1<=0.)||(TMath::Abs(dTmp1-fS)>dTmp1*1e-5)) { cout<<"Can't load Born data file: different energy."; BrdOK=false; } }
			}
			if(BrdOK) {
				ifsB>>iTmp1>>iTmp2>>iTmp3>>iTmp4;
				if(!ifsB.good()){ BrdOK=false; }
				else { if ((iTmp1!=fNPTSX) || (iTmp2!=fNPTSY) || (iTmp3!=fNMODEX) || (iTmp4!=fNMODEY)) { cout<<"Can't load Born data file: different grid."; BrdOK=false; } }
			}
			if(BrdOK) {
				ifsB>>dTmp1>>dTmp2>>dTmp3>>dTmp4;
				if(!ifsB.good()){ BrdOK=false; }
				else {
					if ((dTmp1<=0.) || (TMath::Abs(dTmp1-fXPTMIN)>dTmp1*1e-5) ||
						(dTmp2<=0.) || (TMath::Abs(dTmp2-fXPTMAX)>dTmp2*1e-5) ||
						(dTmp3<=0.) || (TMath::Abs(dTmp3-fYPTMIN)>dTmp3*1e-5) ||
						(dTmp4<=0.) || (TMath::Abs(dTmp4-fYPTMAX)>dTmp4*1e-5)) { cout<<"Can't load Born data file: different grid margins."; BrdOK=false; }
				}
			}
			for(i1=0; ((i1<fNPTSX) && (BrdOK)); i1++) {
				for(i2=0; ((i2<fNPTSY) && (BrdOK)); i2++) {
					ifsB>>dTmp1; if(!ifsB.good()){ BrdOK=false; break; }
					fBorn[i1][i2]=dTmp1;
				}
			}
			if(BrdOK) {
				ifsB>>iTmp1;
				if(!ifsB.good()){ BrdOK=false; }
				else if(iTmp1!=fNCTRLSUM) {cout<<"Can't load Born data file: invalid control sum."; BrdOK=false;}
			}
			ifsB.close();
			if(BrdOK) fBornOK=true;
		} else {
			cout<<"Can't open Born data file: "<<fn_born<<endl;
		}
	}
	if ((!fBornOK) && (fNEVTOT_BORN==0)) { cout<<"No Born data loaded. fNEVTOT_BORN: " <<fNEVTOT_BORN << endl; return false; }

	fNSTATUS=1;
	cout << "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 << "fNEVTOT_BORN="<<fNEVTOT_BORN<<", fNEVTOT_CORR="<<fNEVTOT_CORR<<"\n";
	cout << "fFBINWX="<<fFBINWX<<", fFBINWY="<<fFBINWY<<"\n";
	cout << "C/s per event [pb]: (RC)"<<fCS_Evt<<", (Born)"<<fCS_EvtBorn<<"\n";
	cout << "Born cross section table loaded: "<<fBornOK<<"\n";
	cout << "Additional panes loaded: "<<fNEVTS_AddN<<"\n";
	return true;
}

Bool_t TDjResEvt::IsReady() { return (fNSTATUS==1); }

Bool_t TDjResEvt::GetValue(Double_t x, Double_t y, Double_t* res, Double_t* err, Int_t nAddTbl) {
	Double_t errSta=0.; Double_t errSys=fDJ_ERR_UNKNOWN; if(err)*err=fDJ_ERR_UNKNOWN;
	Bool_t bRes=GetValue(x, y, res, &errSta, &errSys, nAddTbl);
	if((bRes)&&(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 TDjResEvt::GetValue(Double_t x, Double_t y, Double_t* res, Double_t* errSta, Double_t* errSys, Int_t nAddTbl) {
	if(res)*res=fDJ_RES_UNKNOWN; if(errSta)*errSta=0.; if(errSys)*errSys=fDJ_ERR_UNKNOWN; // default values
	if (fNSTATUS != 1) return false;
	if ((nAddTbl<-1) || (nAddTbl>=fNEVTS_AddN)) return false;
	Int_t nxleft,nyleft;
	if(!GetLeftPoint(x,y,&nxleft,&nyleft,true)) return false;

	// Get values & errors of the four nearest points
	Double_t valLL,valRL,valLR,valRR;
	if (!CalcDelta(nxleft,nyleft,&valLL,nAddTbl)) valLL=fDJ_RES_UNKNOWN;
	if (!CalcDelta(nxleft+1,nyleft,&valRL,nAddTbl)) valRL=fDJ_RES_UNKNOWN;
	if (!CalcDelta(nxleft,nyleft+1,&valLR,nAddTbl)) valLR=fDJ_RES_UNKNOWN;
	if (!CalcDelta(nxleft+1,nyleft+1,&valRR,nAddTbl)) valRR=fDJ_RES_UNKNOWN;

	Double_t staLL,staRL,staLR,staRR,sysLL,sysRL,sysLR,sysRR;
	CalcError(nxleft,nyleft,&staLL,&sysLL,valLL,nAddTbl);
	CalcError(nxleft+1,nyleft,&staRL,&sysRL,valRL,nAddTbl);
	CalcError(nxleft,nyleft+1,&staLR,&sysLR,valLR,nAddTbl);
	CalcError(nxleft+1,nyleft+1,&staRR,&sysRR,valRR,nAddTbl);

	// 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 TDjResEvt::GetValueI(Int_t nx, Int_t ny, Double_t* res, Double_t* errSta, Double_t* errSys, Int_t nAddTbl) {
	if(res)*res=fDJ_RES_UNKNOWN; if(errSta)*errSta=0.; if(errSys)*errSys=fDJ_ERR_UNKNOWN; // default values
	if (fNSTATUS != 1) return false;
	if ((nAddTbl<-1) || (nAddTbl>=fNEVTS_AddN)) return false;
	if (!CalcDelta(nx,ny,res,nAddTbl)) return false;
	if (!CalcError(nx,ny,errSta,errSys,*res,nAddTbl)) return false;
	return true;
}

Double_t TDjResEvt::GetBornCS(Double_t x, Double_t y) {
	Double_t res=-1;
	if (fNSTATUS != 1) return res;
	Int_t nxleft,nyleft;
	if(!GetLeftPoint(x,y,&nxleft,&nyleft,true)) return res;

	// Get values & errors of the four nearest points
	Double_t valLL,valRL,valLR,valRR,valAL,valAR,valAA;
	valLL=fBorn[nxleft][nyleft]; valRL=fBorn[nxleft+1][nyleft]; valLR=fBorn[nxleft][nyleft+1]; valRR=fBorn[nxleft+1][nyleft+1];

	// Check points for 'unknown'
	if ((valLL<0)&&(valLR<0)&&(valRL<0)&&(valRR<0)) return res;
	if (valLL<0) { if(valLR>=0){valLL=valLR; } else if(valRL>=0){valLL=valRL; } else if(valRR>=0){valLL=valRR; } }
	if (valLR<0) { if(valLL>=0){valLR=valLL;} else if(valRR>=0){valLR=valRR;} else if(valRL>=0){valLR=valRL;} }
	if (valRL<0) { if(valLL>=0){valRL=valLL;} else if(valRR>=0){valRL=valRR;} else if(valLR>=0){valRL=valLR;} }
	if (valRR<0) { if(valLR>=0){valRR=valLR;} else if(valRL>=0){valRR=valRL;} else if(valLL>=0){valRR=valLL;} }

	// 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;
	if(dx>1e-14){ valAL=(valLL*dxr+valRL*dxl)/dx; valAR=(valLR*dxr+valRR*dxl)/dx; }else{ valAL=(valLL+valRL)/2; valAR=(valLR+valRR)/2; }
	if(dy>1e-14){ valAA=(valAL*dyr+valAR*dyl)/dy; }else{ valAA=(valAL+valAR)/2; }

	res=valAA;
	return res;
}

Bool_t TDjResEvt::CalcDelta(Int_t nx, Int_t ny, Double_t* delta, Int_t nAddTbl) {
	if (delta)*delta=fDJ_RES_UNKNOWN;
	if ((fNSTATUS != 1) || (nAddTbl<-1) || (nAddTbl>=fNEVTS_AddN)) return false;
	Int_t nevRC=0;
	if(nAddTbl<0){ nevRC=fNEVTS_CORR[nx][ny]; }else{ nevRC=fNEVTS_Add[nx][ny][nAddTbl]; }
	if(nevRC<=0)return false;
	Double_t val=fDJ_RES_UNKNOWN;
	if (fBornOK) {
		Double_t x=fRPTSX[nx]; Double_t y=fRPTSY[ny];
		Double_t born=fBorn[nx][ny]; if(born<=0.) return false;
		if ((fFBINWX<=0.) || (fFBINWY<=0.)) return false;
		if (fCS_Evt<=0.) return false;
		Double_t x2=x*(1.+fFBINWX); Double_t x1=x*(1.-fFBINWX); Double_t y2=y*(1.+fFBINWY); Double_t y1=y*(1.-fFBINWY);
		// Check for borders vicinity
		if(x2>fXmax)x2=fXmax; if(x1<fXmin)x1=fXmin;
		if (fCutMode==3) { if(y2>fYmax)y2=fYmax; if(y1<fYmin)y1=fYmin; if((fS*x2*y2)>fQ2max)return false; }
		if((fS*x1*y1)<fQ2min)return false;
		Double_t dx=x2-x1; Double_t dy=y2-y1; if((dx<0)||(dy<0))return false;
		val=nevRC; val=val*fCS_Evt/dx/dy/born-1.;
	} else {	
		Int_t nb=fNEVTS_BORN[nx][ny]; if(nb<=0) return false;
		val=nevRC; val=val*fNEVTOT_BORN/fNEVTOT_CORR/nb-1.;
	}
	if (delta)*delta=val;
	return (val!=fDJ_RES_UNKNOWN);
}

Int_t TDjResEvt::GetEvNumBorn(Double_t x, Double_t y) {
	if (fNSTATUS!=1) return 0;
	Int_t nxleft,nyleft;
	if(!GetLeftPoint(x,y,&nxleft,&nyleft,false))return 0;
	Int_t res=fNEVTS_BORN[nxleft][nyleft];
	return res;
}

Int_t TDjResEvt::GetEvNum(Double_t x, Double_t y, Int_t nAddTbl) {
	if (fNSTATUS != 1) return 0;
	if ((nAddTbl<-1) || (nAddTbl>=fNEVTS_AddN)) return 0;
	Int_t nxleft,nyleft;
	if(!GetLeftPoint(x,y,&nxleft,&nyleft,false))return 0;
	Int_t res=0;
	if (nAddTbl<0) { res=fNEVTS_CORR[nxleft][nyleft]; }
	else { res=fNEVTS_Add[nxleft][nyleft][nAddTbl]; }
	return res;
}

Bool_t TDjResEvt::CalcError(Int_t nx, Int_t ny, Double_t* errSta, Double_t* errSys, Double_t delta, Int_t nAddTbl) {
	if(errSta)*errSta=0.; if(errSys)*errSys=fDJ_ERR_UNKNOWN; //default error values
	if ((delta<-1.) || (delta==fDJ_RES_UNKNOWN)) return false;
	if ((fNSTATUS != 1) || (nAddTbl<-1) || (nAddTbl>=fNEVTS_AddN)) return false;
	Double_t x=fRPTSX[nx]; Double_t y=fRPTSY[ny];
	Int_t Ne=0; if(nAddTbl==-1){ Ne=fNEVTS_CORR[nx][ny]; }else{ Ne=fNEVTS_Add[nx][ny][nAddTbl]; }
	Int_t Nb=0; if(!fBornOK)Nb=fNEVTS_BORN[nx][ny];
	Bool_t bRes=false;
	if ((Ne>0)&&(fNEVTOT_CORR>Ne)&&( (fBornOK) || ((Nb>0)&&(fNEVTOT_BORN>Nb)) )) {
		Double_t sigStaC2=0.; Double_t sigStaB2=0.; Double_t sigSysC2=0.; Double_t sigSysB2=0.;

		// Statistical error
		// Nt: SigmaC=Corr*Sqrt(1/Nbin-1/Nevts)~Corr*Sqrt(1/Nbin)
		// Nt: SigmaB=Born*Sqrt(1/Nbin-1/Nevts)~Born*Sqrt(1/Nbin)
		sigStaC2=sigStaC2+1./Ne-1./fNEVTOT_CORR;
		if(!fBornOK) { sigStaB2=sigStaB2+1./Nb-1./fNEVTOT_BORN; }

		//Systematical error (due to fixed bin width)
		//For correction table
		Double_t err_linx=0; Double_t err_liny=0; Double_t err_lin=0;
		Int_t neXL=GetEvNum(x*(1.-fFBINWX), y, nAddTbl); Int_t neXR=GetEvNum(x*(1.+fFBINWX), y, nAddTbl);
		Int_t neYL=GetEvNum(x, y*(1.-fFBINWY), nAddTbl); Int_t neYR=GetEvNum(x, y*(1.+fFBINWY), nAddTbl);
		Int_t nx_tmp,ny_tmp;
		if ((neXL>0) && (neXR>0)) { err_linx=((Double_t)(neXL+neXR-2*Ne))/2.; }
		else if (neXL>0) { err_linx=neXL-Ne; }
		else if (neXR>0) { err_linx=neXR-Ne; }
		else { err_linx=0; }
		if ((neYL>0) && (neYR>0)) { err_liny=((Double_t)(neYL+neYR-2*Ne))/2.; }
		else if (neYL>0) { err_liny=neYL-Ne; }
		else if (neYR>0) { err_liny=neYR-Ne; }
		else { err_liny=0; }
		//err_lin=(err_linx+err_liny)/2/Ne;
		if (err_linx<0)err_linx=-err_linx; if (err_liny<0)err_liny=-err_liny;
		if (err_linx>err_liny){ err_lin=err_linx/Ne; }else{ err_lin=err_liny/Ne; } //err_lin=max(abs(err_linx),abs(err_liny))
		if (err_lin<0)err_lin=-err_lin; if(err_lin>1)err_lin=1;
		sigSysC2=sigSysC2+err_lin*err_lin;

		if(!fBornOK) {
			// Sys. error for Born table
			err_linx=0; err_liny=0;
			Int_t nbXL=GetEvNumBorn(x*(1.-fFBINWX), y); Int_t nbXR=GetEvNumBorn(x*(1.+fFBINWX), y);
			Int_t nbYL=GetEvNumBorn(x, y*(1.-fFBINWY)); Int_t nbYR=GetEvNumBorn(x, y*(1.+fFBINWY));
			if ((nbXL>0) && (nbXR>0)) { err_linx=((Double_t)(nbXL+nbXR-2*Nb))/2; }
			else if (nbXL>0) { err_linx=nbXL-Nb; }
			else if (nbXR>0) { err_linx=nbXR-Nb; }
			else { err_linx=0; }
			if ((nbYL>0) && (nbYR>0)) { err_liny=((Double_t)(2*Nb-nbYL-nbYR))/2; }
			else if (nbYL>0) { err_liny=Nb-nbYL; }
			else if (nbYR>0) { err_liny=Nb-nbYR; }
			else { err_liny=0; }
			//err_lin=(err_linx+err_liny)/2/Nb;
			if (err_linx<0)err_linx=-err_linx; if (err_liny<0)err_liny=-err_liny;
			if (err_linx>err_liny){ err_lin=err_linx/Nb; }else{ err_lin=err_liny/Nb; } //err_lin=max(abs(err_linx),abs(err_liny))
			if(err_lin<0)err_lin=-err_lin; if(err_lin>1)err_lin=1;
			sigSysB2=sigSysB2+err_lin*err_lin;
		}

		// Total error
		// Nt: SigmaDelta=Corr/Born*Sqrt((SigmaC/Corr)^2+(SigmaB/Born)^2)
		Double_t csRatio=delta+1.; //=((Double_t)Ne)*fNEVTOT_BORN/fNEVTOT_CORR/Nb;
		Double_t resSta=TMath::Sqrt(sigStaC2+sigStaB2)*csRatio;
		Double_t resSys=TMath::Sqrt(sigSysC2+sigSysB2)*csRatio;

		if(errSta)*errSta=resSta; if(errSys)*errSys=resSys;
		bRes=true;
	}
	return bRes;
}

Bool_t TDjResEvt::GetLeftPoint(Double_t x, Double_t y, Int_t *nxleft, Int_t *nyleft, Bool_t bNeedRightOK) {
	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;
	}
	if (bNeedRightOK) { // The right point (nxleft+1) should also be valid.
		if ((nxlTmp==(fNPTSX-1))&&(nxlTmp>0))nxlTmp=nxlTmp-1;
	}
	//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 (bNeedRightOK) { // The right point (nyleft+1) should also be valid.
		if ((nylTmp==(fNPTSY-1))&&(nylTmp>0))nylTmp=nylTmp-1;
	}
	if (nxleft) *nxleft=nxlTmp;
	if (nyleft) *nyleft=nylTmp;
	return true;
}

/*
*	Loads pre-calculated radiative correction (delta) results.
*	To be loaded after event file loading only.
*	The points array should be already filled.
*	Use TDjRes class to work with results file only.
*/
Bool_t TDjResEvt::LoadDjRes(const char* fn) {
	if (fNSTATUS!=1) { cout << "WARNING: can't load RESULTS file before EVENTS file is loaded\n"; return false; }
	if (!fn) return false;
	cout << "Loading pre-calculated results file: " << fn << "\n";
	if (fRESArrayFilled) cout << "WARNING: previous DjRes data will be overwritten\n";
	fRESArrayFilled=false;
/*
	----- Results file format ----- ver.2
   #(int) IVER: file version (reserved, currently 2)                 *
   #(int) I2ERRS: 1(0) the stat. & sys. errors are separated         *
   #(int) S: S-Mp2-Ml2                                               *
   #(int) NPTSX,NPTSY,NMODEX,NMODEY: number of points and grid mode  *
   #(double) XPTMIN,XPTMAX,YPTMIN,YPTMAX: min/max point values        *
   #(int) DELTA,ERRSTA[,ERRSYS] x(NPTSX*NPTSY): results array        *
   #(int) NCTRLSUM: currently NPTSX*NPTSY                            *
*/
	Int_t i1,i2;
	Int_t iTmp1,iTmp2,iTmp3,iTmp4;
	Double_t dTmp1,dTmp2,dTmp3,dTmp4;

	Int_t iver=0;
	ifstream ifs; ifs.open(fn); if(!ifs){cout<<"Error loading results file. No file found.\n"; return false;}
	ifs>>iver; if(!ifs.good()){cout<<"Error loading results file.\n"; ifs.close();return false;}
	if (iver!=2) { cout<<"Error loading results file. Unknown format."; ifs.close(); return false; }
	// 2-errors flag
	Int_t iKeep2Errs=1;
	ifs>>iKeep2Errs; if(!ifs.good()){ cout<<"Can't load results file. iKeep2Errs parameter expected.\n"; ifs.close();return false; }
	if ((iKeep2Errs<0) || (iKeep2Errs>1)) { cout<<"Can't load results file: unexpected iKeep2Errs parameter value, "<<iKeep2Errs<<"\n"; ifs.close(); return false; }
	// energy
	ifs>>dTmp1; if(!ifs.good()){ cout<<"Can't load results file. Stopped on parameter #"<<i1<<"\n"; ifs.close();return false; }
	if ((dTmp1<=0.)||(TMath::Abs(dTmp1-fS)>dTmp1*1e-5)) { cout<<"Can't load results file: different energy\n"; ifs.close(); return false; }
	// points grid
	ifs>>iTmp1>>iTmp2>>iTmp3>>iTmp4; if(!ifs.good()){cout<<"Can't load results file. Stopped on points grid data.\n";ifs.close();return false;} // Load data table size
	if ((iTmp1!=fNPTSX) || (iTmp2!=fNPTSY) || (iTmp3!=fNMODEX) || (iTmp4!=fNMODEY)) {
		cout<<"Can't load results file: different points grid settings to current one\n"; ifs.close(); return false;
	}
	ifs>>dTmp1>>dTmp2>>dTmp3>>dTmp4; if(!ifs.good()){cout<<"Can't load results file. Stopped on points margins.\n";ifs.close();return false;}
	if ((TMath::Abs(dTmp1-fXPTMIN)>1E-6) || (TMath::Abs(dTmp2-fXPTMAX)>1E-6) ||
		(TMath::Abs(dTmp3-fYPTMIN)>1E-6) || (TMath::Abs(dTmp4-fYPTMAX)>1E-6)) {
		cout<<"Can't load results file: different points grid settings to current one\n"; ifs.close(); return false;
	}
	// data
	for(i1=0; i1<fNPTSX; i1++) {
		for(i2=0; i2<fNPTSY; i2++) {
			ifs>>dTmp1>>dTmp2>>dTmp3; if(!ifs.good()){cout<<"Can't load results file. Table loaded up to i1="<<i1<<", i2="<<i2<<"\n";ifs.close();return false;}
			fRESDELTA[i1][i2]=dTmp1;
			fRESERRSTA[i1][i2]=dTmp2;
			if (iKeep2Errs==1) { fRESERRSYS[i1][i2]=dTmp3; } else { fRESERRSYS[i1][i2]=0.; }
		}
	}
	ifs >> iTmp1;
	ifs.close();
	// Check data control sum (currently number of numbers read).
	if (iTmp1 != (fNPTSX*fNPTSY)) { cout << "Error results file loading: Control sum failed\n"; return false; }
	fRESArrayFilled=true;
	cout << "Results file loaded.\n";
	return true;
}

/*
*	Update previously loaded results or create new results array
*/
Bool_t TDjResEvt::UpdateDjRes(Int_t nAddTbl) {
	if (fNSTATUS!=1) return false;
	Int_t i1,i2;
	for(i1=0; i1<fNPTSX; i1++) {
		for(i2=0; i2<fNPTSY; i2++) {
			Double_t delta=fDJ_RES_UNKNOWN; Double_t err=fDJ_ERR_UNKNOWN; Double_t errSta=fDJ_ERR_UNKNOWN; Double_t errSys=fDJ_ERR_UNKNOWN;
			Bool_t cuts_OK=true;
			Double_t Q2=fS*fRPTSX[i1]*fRPTSY[i2];
			if (Q2*(1.0-fFBINWX)*(1.0-fFBINWY) < fQ2min) cuts_OK=false;
			if (fRPTSX[i1]*(1.0+fFBINWX)>fXmax) cuts_OK=false;
			if (fCutMode==3) { if (fRPTSY[i2]*(1.0-fFBINWY) < fYmin) cuts_OK=false; }
			if (cuts_OK) {
				if (GetValueI(i1,i2,&delta,&errSta,&errSys,nAddTbl)) {
				//if (GetValue(fRPTSX[i1],fRPTSY[i2],&delta,&errSta,&errSys,nAddTbl)) {
					if ((errSta<fDJ_ERR_UNKNOWN) && (errSys<fDJ_ERR_UNKNOWN)) { err=TMath::Sqrt(errSta*errSta+errSys*errSys); }
					else { err=fDJ_ERR_UNKNOWN; }
				} else {
					delta=fDJ_RES_UNKNOWN; err=fDJ_ERR_UNKNOWN; errSta=fDJ_ERR_UNKNOWN; errSys=fDJ_ERR_UNKNOWN;
				}
			}
			// Note: if no previous data exist, we compare error to 'old_error=1.'
			Double_t errOld=fDJ_ERR_UNKNOWN;
			if (fRESArrayFilled) {
				if((fRESERRSTA[i1][i2]<fDJ_ERR_UNKNOWN) && (fRESERRSYS[i1][i2]<fDJ_ERR_UNKNOWN)) {
					errOld=TMath::Sqrt(fRESERRSTA[i1][i2]*fRESERRSTA[i1][i2]+fRESERRSYS[i1][i2]*fRESERRSYS[i1][i2]);
				}
			}
			if (err<errOld) {
				fRESDELTA[i1][i2]=delta;
				fRESERRSTA[i1][i2]=errSta;
				fRESERRSYS[i1][i2]=errSys;
			} else if (!fRESArrayFilled) { // fill 'no data' values
				fRESDELTA[i1][i2]=fDJ_RES_UNKNOWN;
				fRESERRSTA[i1][i2]=fDJ_ERR_UNKNOWN;
				fRESERRSYS[i1][i2]=fDJ_ERR_UNKNOWN;
			}
		}
	}
	fRESArrayFilled=true;
	return true;
}

/*
*	Saves radiative correction (delta) results.
*	The same procedure is embeeded into DJANGOH user procedure (djangoh_u_r.f)
*	Use this function to counter-check of DJANGOH user procedure only.
*/
Bool_t TDjResEvt::SaveDjRes(const char* fn) {
	if ((!fRESArrayFilled) || (!fn)) return false;
	cout<<"Saving results file.\n";
/*
	----- Results file format ----- ver.2
   #(int) IVER: file version (reserved, currently 2)                 *
   #(int) I2ERRS: 1(0) the stat. & sys. errors are separated         *
   #(int) S: S-Mp2-Ml2                                               *
   #(int) NPTSX,NPTSY,NMODEX,NMODEY: number of points and grid mode  *
   #(double) XPTMIN,XPTMAX,YPTMIN,YPTMAX: min/max point values        *
   #(int) DELTA,ERRSTA[,ERRSYS] x(NPTSX*NPTSY): results array        *
   #(int) NCTRLSUM: currently NPTSX*NPTSY                            *
*/
	Int_t i1,i2;
	ofstream ofs(fn); //ios::trunc
	if(!ofs){cout<<"Error: can't create results file.\n";return false;}
	ofs<<2<<endl; //version
	ofs<<1<<endl; //keep errors separated
	ofs<<fS<<endl; //parameter
	ofs<<fNPTSX<<endl;
	ofs<<fNPTSY<<endl;
	ofs<<fNMODEX<<endl;
	ofs<<fNMODEY<<endl;
	ofs<<fXPTMIN<<endl;
	ofs<<fXPTMAX<<endl;
	ofs<<fYPTMIN<<endl;
	ofs<<fYPTMAX<<endl;
	for(i1=0; i1<fNPTSX; i1++) {
		for(i2=0; i2<fNPTSY; i2++) {
			ofs<<fRESDELTA[i1][i2]<<endl;
			ofs<<fRESERRSTA[i1][i2]<<endl;
			ofs<<fRESERRSYS[i1][i2]<<endl;
		}
	}
	ofs<<(fNPTSX*fNPTSY)<<endl; //control sum
	ofs.close();
	cout << "Results file saved.\n";
	return true;
}

/******************************************************************************
***                                                                         ***
***                               Examples                                  ***
***                                                                         ***
******************************************************************************/

/*
*	Example #1
*	Plot 1D histograms for ZEUS, JB and leptonic delta.
*	'x' is fixed to 0.01,0.03,0.1 and 0.4.
*	See parameters definition in function header.
*/
#define HIST_nNy	100
int DjResEvtExample() {
//parameters
	Double_t x[]={0.01,0.03,0.1,0.4};
	Int_t numX=4; // number of plots (items in x[.], hist[][.] and errors[Ny][][.] arrays).
	Bool_t bPlotNumEvents=true; // plot events numbers instead of delta
	Bool_t bViewBorn=true; // view Born events estimation (on number-of events plot only)
// Table numbers & titles
// Note. The (default) order of data tables:
// #-1:ZEUS[ThPipe=5],#0(1):ZEUS[ThPipe=0],#1(2):ZEUS[ThPipe=3],#2(3):ZEUS[ThPipe=10],#3(4):JB,#4(5):Hadronic,#5(6):Leptonic
	Int_t tbl[3]={5,3,-1};
	const char* tblName[3]={"Lept","JB","ZEUS"};
	Int_t clrs[3]={kGreen,kRed,kBlue};
	Int_t numT=3; // size of tbl[], tblName[], clrs[], hist[.][] and errors[Ny][.][] arrays
//-----------------
	// init
	cout<<"\DjResEvtExample() Started\n";
	TDjResEvt data;
	Bool_t bOK = data.Load("djangoh_r.evts","djangoh_r.born");
	if (!bOK) { cout << "data.Load() failed\n"; return 0; }
	if (!data.fBornOK) { cout << "WARNING: data.Load() failed to load Born cross section table\n"; }
	Int_t i,ix,it;
	char buf[50],buf2[50];
	const char* yTitle="delta"; if(bPlotNumEvents)yTitle="N events";

	// plot layout (split canvas for nvert*nhor cells)
	Int_t nhor=2; if(numX==1)nhor=1; Int_t nvert=numX/nhor; if((nvert*nhor)<numX)nvert++;
	TCanvas *cl=new TCanvas("cl","Example",10,10,nhor*400,((nvert>3)?1200:nvert*400));

	// y-margins
	Double_t yMin=data.fYPTMIN; Double_t yMax=data.fYPTMAX; Int_t Ny=HIST_nNy;
	// Shift some plots for half-of-bin for better view:
	Double_t yMinBis=yMin+(yMax-yMin)*0.5/(Ny-1); Double_t yMaxBis=yMinBis+(yMax-yMin);

	// book histos
	Double_t errors[3][4][HIST_nNy]; //size: [points][numT][numX]
	TH1F *hist[3][4]; // size: [numT][numX]
	TH1F *histB[4]; // Born line plots, size: [numX]
	for(ix=0; ix<numX; ix++) { // x
		for(it=0; it<numT; it++) { // table
			sprintf(buf,"%s %s x=%.2f",yTitle,tblName[it],x[ix]); sprintf(buf2,"%s, x=%.2f",yTitle,x[ix]);
			hist[it][ix]=new TH1F(buf,buf2,Ny,((it<(numT-1))?yMin:yMinBis),((it<(numT-1))?yMax:yMaxBis)); // shift bins for last plot
			hist[it][ix]->GetXaxis()->SetTitle("y"); hist[it][ix]->GetYaxis()->SetTitle(yTitle); hist[it][ix]->SetStats(false);
		}
		if (bViewBorn && bPlotNumEvents) {
			sprintf(buf,"%s BORN x=%.2f",yTitle,x[ix]); sprintf(buf2,"%s, x=%.2f",yTitle,x[ix]);
			histB[ix]=new TH1F(buf,buf2,Ny,yMin,yMax);
			histB[ix]->GetXaxis()->SetTitle("y"); histB[ix]->GetYaxis()->SetTitle(yTitle); histB[ix]->SetStats(false);
		}
	}

	// fill histos
	for(i=0; i<Ny; i++) {
		Double_t y=yMin+(yMax-yMin)*(0.5+i)/Ny; // the middle of y-bin
		Double_t yBis=yMinBis+(yMaxBis-yMinBis)*(0.5+i)/Ny; // the middle of y-bin
		for(ix=0; ix<numX; ix++) { // x
			for(it=0; it<numT; it++) { // table
				if(bPlotNumEvents) { // fill with number of events
					hist[it][ix]->Fill(y,data.GetEvNum(x[ix],((it<(numT-1))?y:yBis),tbl[it]));
				} else { // fill with delta and error
					Double_t res=0.; Double_t err=1.;
					if(data.GetValue(x[ix],((it<(numT-1))?y:yBis),&res,&err,tbl[it])) {
						hist[it][ix]->Fill(((it<(numT-1))?y:yBis),res); errors[it][ix][i]=err;
					}
				}
			}
			if (bViewBorn && bPlotNumEvents) {
				Double_t val=data.GetBornCS(x[ix],y);
				if ((val>=0) && (data.fCS_Evt>0)) {
					val=val/data.fCS_Evt*x[ix]*2*data.fFBINWX*y*2*data.fFBINWY; histB[ix]->Fill(y,val);
				}
			}
		}
	}

	if(!bPlotNumEvents) {
		for(ix=0; ix<numX; ix++) { // x
			for(it=0; it<numT; it++) { // table
				hist[it][ix]->SetError(errors[it][ix]);
				hist[it][ix]->SetMinimum(-0.2); hist[it][ix]->SetMaximum(0.2);
			}
			if (bViewBorn && bPlotNumEvents) { histB[ix]->SetMinimum(-0.2); histB[ix]->SetMaximum(0.2); }
		}
	}

	// plot
	cl->Divide(nhor,nvert);
	for(ix=0; ix<numX; ix++) { // x
		cl->cd(ix+1);
		for(it=0; it<numT; it++) { // table
			hist[it][ix]->SetLineColor(clrs[it]); hist[it][ix]->Draw(((it==0)?"":"same"));
		}
		if (bViewBorn && bPlotNumEvents) { histB[ix]->Draw("C same"); }
	}
	cl->cd();
	cl->SaveAs("DjResEvtExample.eps");
	cout<<"\DjResEvtExample() Finished\n";
	return 1;
}

/*
*	Example #2
*	Combine several ".evts" files into single ".dat" file.
*	Note: the same procedure is performed automatically in Djangoh user Procedure.
*	iTable==-1: -> ZEUS variables
*	iTable==3: -> JB variables
*	iTable==5: -> Leptonic variables
*/
int DjResEvtExampleConvert() {
	cout<<"DjResEvtExampleConvert() Started\n";
	Int_t iTable=-1;

// let's combine results of 8 runs
	TDjResEvt arrData[8]; Int_t n=8;
	arrData[0].Load("djangoh_r_1.evts","djangoh_r.born");
	arrData[1].Load("djangoh_r_2.evts","djangoh_r.born");
	arrData[2].Load("djangoh_r_3.evts","djangoh_r.born");
	arrData[3].Load("djangoh_r_4.evts","djangoh_r.born");
	arrData[4].Load("djangoh_r_5.evts","djangoh_r.born");
	arrData[5].Load("djangoh_r_6.evts","djangoh_r.born");
	arrData[6].Load("djangoh_r_7.evts","djangoh_r.born");
	arrData[7].Load("djangoh_r_8.evts","djangoh_r.born");

	char fn_res[50]; sprintf(fn_res,"djangoh_r_tbl%d.dat...",iTable);
	cout << "Trying to combine event files to "<<fn_res<<endl;
	for(Int_t i=0; i<n; i++) {
		if(!arrData[i].IsReady()) {
			cout << "ERROR. Events file #"<<i+1<<" is not ready. Skipped."<<endl;
			continue;
		}
		cout << "Using events file #"<<i+1<<"..."<<endl;
		if(i>0) arrData[i].LoadDjRes(fn_res);
		arrData[i].UpdateDjRes(iTable); // Update results using events data of arrData[i]
		arrData[i].SaveDjRes(fn_res);
	}
	cout<<"DjResEvtExampleConvert() Finished\n";
	return 1;
}