#include "TCanvas.h"
#include "TF1.h"
#include "TFile.h"
#include "TH1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TLatex.h"
#include "TLine.h"
#include "TObjArray.h"
#include "TObjString.h"
#include "TPad.h"
#include "TROOT.h"
#include "TString.h"
#include "TTree.h"
#include "TArrow.h"
#include "TLegend.h"
#include "TMath.h"
#include "TGraphErrors.h"
#include "TGraphAsymmErrors.h"
#include "TRandom3.h"
#include "TPaveStats.h"
#include "TStyle.h"

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooPlot.h"

#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooArgusBG.h"
#include "RooNovosibirsk.h"
#include "RooPolynomial.h"

#include "RooProdPdf.h"
#include "RooAddPdf.h"


#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <stdio.h>

using namespace RooFit ;

int cnt=0, rooEvS=0, rooEvB=0, rooNds=0;

// ------------------------
// CONSTANTS
// ------------------------

TString savepath   = "fig/";
TString savepathqa = savepath+"scans/"; 
TString savepref   = "";

// initial event numbers e+ e-
double S0e   = 98000.;          // signal events e+e- mode
double B0hde = 9580637576.;     // hadronic (DPM) bkg events e+e- mode
double B0nre = 100000.;         // non-res (J/psi 2pi) bkg events e+e- mode

// initial event numbers mu+ mu-
double S0m   = 100000.;         // signal events mu+mu- mode
double B0hdm = 8871910789.;     // hadronic (DPM) bkg events mu+mu- mode
double B0nrm = 99000.;          // non-res (J/psi 2pi) bkg events mu+mu- mode

// Branching fractions 
double BRe   = 0.05971;         // BR(J/psi -> e+ e-)
double BRm   = 0.0596;          // BR(J/psi -> mu+ mu-)
double BRX   = 0.05;            // BR(X -> J/psi pi+ pi-) assumption

// cross section values
double SigXs   = 100;            // [nb] maximum nominal signal cross section at peak
double BkgXs   = 46e6;          // [nb] inelastic
double BkgXst  = 59e6;          // [nb] total
double BkgXsNR = 1.2;           // [nb] J/psi pi+ pi- NR

// Integrated luminosities [1/(nb·d)] at Ecm = 3.872 GeV, high R;
// https://panda.gsi.de/publication/in-ide-2015-002
double IntLHL    = 13683;  
double IntLHR    = 1368; 
double IntLHESR  = 972; 
double IntLHESRr = 1170;

double dEcmHL = 0.1678;     // [MeV] E_cm resolution in HL mode @ 3.872 GeV (dp/p = 1.0e-4)
double dEcmHR = 0.0839;     // [MeV] E_cm resolution in HR mode @ 3.872 GeV (dp/p = 0.5e-4) <==== !!

double dEBeamAbs = 0.167;   // 1E-4 energy uncertainty for first scan point
double dEBeamRel = 0.00167; // 1E-6 energy uncertainty for going from point to point and the reproduction of the interval start point 1E-6
double currEshift;

double freqFixWin = false;

// ----------------------------------------------------------
// msum : ntp2->SetAlias("msum","(f4cxd0m+f4cxd1m)");
// ----------------------------------------------------------

// efficiencies J/psi -> e+ e- (MJ2e_S98k_NR100k_HD9580M_ntp2.root)
// cut and count : cut = "abs(f4cxd0m-3.093)<0.05 && msum>3.77  &&xd0d0pide>0.95&&xd0d1pide>0.95&&f4cxm>3.867&&f4cxm<3.874&&xd0oang<2.1&&xpcm<0.4"
// ML J/psi      : cut = "abs(f4cxd0m-3.093)<0.20 && msum>3.77  &&xd0d0pide>0.95&&xd0d1pide>0.95&&f4cxm>3.867&&f4cxm<3.874&&xd0oang<2.1&&xpcm<0.4"
// ML m_sum/2D   : cut = "abs(f4cxd0m-3.093)<0.05 && msum>3.60  &&xd0d0pide>0.95&&xd0d1pide>0.95&&f4cxm>3.867&&f4cxm<3.874&&xd0oang<2.1&&xpcm<0.4"

//                     C&C            ML(J/psi)       ML(sum)     ML(2d)
//double effSige[4] = { 10967/S0e,     11915/S0e,     14038/S0e,   14038/S0e   };
//double effBhde[4] = {  0.25/B0hde,     1.0/B0hde,     1.0/B0hde,   1.0/B0hde };
//double effBnre[4] = {  2523/B0nre,    2784/B0nre,    8311/B0nre,  8311/B0nre };

// cut and count : cut = "abs(f4cxd0m-3.093)<0.025 && msum>3.77  &&xd0d0pide>0.95&&xd0d1pide>0.95&&f4cxm>3.867&&f4cxm<3.874&&xd0oang<2.1&&xpcm<0.4"
// ML J/psi      : cut = "abs(f4cxd0m-3.093)<0.20 && msum>3.77  &&xd0d0pide>0.95&&xd0d1pide>0.95&&f4cxm>3.867&&f4cxm<3.874&&xd0oang<2.1&&xpcm<0.4"
// ML m_sum/2D   : cut = "abs(f4cxd0m-3.093)<0.025 && msum>3.60  &&xd0d0pide>0.95&&xd0d1pide>0.95&&f4cxm>3.867&&f4cxm<3.874&&xd0oang<2.1&&xpcm<0.4"
double effSige[4] = {  9947/S0e,     11914/S0e,     12723/S0e,   12723/S0e   };
double effBhde[4] = {  0.125/B0hde,    1.0/B0hde,     0.5/B0hde,   0.5/B0hde };
double effBnre[4] = {  2287/B0nre,    2784/B0nre,    7524/B0nre,  7524/B0nre };

// efficiencies J/psi -> mu+ mu- (MJ2mu_S100k_NR99k_HD8870M_ntp2.root)
// cut and count : cut = "abs(f4cxd0m-3.097)<0.04 && msum>3.78  &&xd0d0pidmu>0.99&&xd0d1pidmu>0.99&&xd0oang<1.4&&essph<0.11"
// ML J/psi      : cut = "abs(f4cxd0m-3.097)<0.16 && msum>3.78  &&xd0d0pidmu>0.99&&xd0d1pidmu>0.99&&xd0oang<1.4&&essph<0.11"
// ML m_sum/2D   : cut = "abs(f4cxd0m-3.097)<0.04 && msum>3.60  &&xd0d0pidmu>0.99&&xd0d1pidmu>0.99&&xd0oang<1.4&&essph<0.11"

////                     C&C            ML(J/psi)       ML(sum)       ML(2d)
//double effSigm[4] = { 13981/S0m,     15171/S0m,     18287/S0m,   18287/S0m   };
//double effBhdm[4] = {   1.0/B0hdm,     4.0/B0hdm,     6.0/B0hdm,   6.0/B0hdm };
//double effBnrm[4] = {  2694/B0nrm,    2955/B0nrm,    9943/B0nrm,  9943/B0nrm };

// cut and count : cut = "abs(f4cxd0m-3.097)<0.02 && msum>3.78  &&xd0d0pidmu>0.99&&xd0d1pidmu>0.99&&xd0oang<1.4&&essph<0.11"
// ML J/psi      : cut = "abs(f4cxd0m-3.097)<0.16 && msum>3.78  &&xd0d0pidmu>0.99&&xd0d1pidmu>0.99&&xd0oang<1.4&&essph<0.11"
// ML m_sum/2D   : cut = "abs(f4cxd0m-3.097)<0.02 && msum>3.60  &&xd0d0pidmu>0.99&&xd0d1pidmu>0.99&&xd0oang<1.4&&essph<0.11"
double effSigm[4] = { 12864/S0m,     15171/S0m,     16778/S0m,   16778/S0m   };
double effBhdm[4] = {   0.5/B0hdm,     4.0/B0hdm,    3.75/B0hdm,  3.75/B0hdm }; //last 2 scaled from 4 times broader window
double effBnrm[4] = {  2476/B0nrm,    2955/B0nrm,    9068/B0nrm,  9068/B0nrm };


// will be set according to mode
double effS, effBhd, effBnr;


// ------------------------
// END CONSTANTS
// ------------------------

// ------------------------
// ROOFIT STUFF 
// ------------------------

	int AnaMode = 0;  // 0 = Cut and Count, 1 = fit in m(ll), 2 = fit in m(ll)+m(pipi), 3 = fit in combined 2d mode
 
// ------------------------
// ROOFIT STUFF m(pipi)+m(ll) ML
// ------------------------

	RooRealVar rx("rx","x", 3.6,        3.9);
	RooRealVar ry("ry","m(l^{+}l^{-}) [GeV/c^{2}]", 3.097-0.02, 3.097+0.02);
	
	//Signal D-Gauss : 3.09774 0.00727852 0.0338042 3.32509 
	//Signal Novo    : 3.854 0.0175444 0.840566 
	
	// *******
	// signal
	// *******
	
	// double gaus par
	RooRealVar dgmus("dgmus","", 3.09773 );
	RooRealVar dgs1s("dgs1s","", 7.27532e-03);
	RooRealVar dgs2s("dgs2s","", 4.55356e-02);
	RooRealVar dgfrs("dgfrs","", 3.57695/(1+3.57695));
   //2  #mu          3.09773e+00   8.96295e-05   7.38556e-07  -6.28637e-01
   //3  #sigma_{1}   7.27532e-03   9.63943e-05   2.01527e-07  -1.12163e+00
   //4  #sigma_{2}   4.55356e-02   1.28700e-03   2.55581e-06  -6.29619e-02
   //5  R            3.57695e+00   1.30498e-01   5.05287e-05   7.88318e-04
	//RooRealVar dgmus("dgmus","", 3.09774 );
	//RooRealVar dgs1s("dgs1s","", 0.00727852 );
	//RooRealVar dgs2s("dgs2s","", 0.0338042);
	//RooRealVar dgfrs("dgfrs","", 3.32509/(1+3.32509));
	
	// novosibirsk par
	RooRealVar nvmus("nvmus",  "", 3.854);
	RooRealVar nvsigs("nvsigs","", 0.0175444);
	RooRealVar nvtaus("nvtaus","", 0.840566);

	// double gaus pdf
	RooGaussian g1s("g1s","",ry,dgmus,dgs1s);  
	RooGaussian g2s("g2s","",ry,dgmus,dgs2s);  
	RooAddPdf   dgs("dgs","",RooArgList(g1s,g2s),dgfrs);
	
	// novosibirsk pdf
	RooNovosibirsk nvs("nvs","",rx,nvmus,nvsigs,nvtaus);
	
	// combined
	RooProdPdf pdfs("pdfs","",dgs,nvs);
	
	// *******
	// bg nr
	// *******
	//Bkg NR D-Gauss : 3.09773 0.0065123 0.0297225 3.1027 
	//Bkg NR Argus   : 3.87 0.503262 -4.91777 
	
	// double gaus par
	RooRealVar dgmub1("dgmub1","", 3.09773);
	RooRealVar dgs1b1("dgs1b1","", 0.0065123);
	RooRealVar dgs2b1("dgs2b1","", 0.0297225);
	RooRealVar dgfrb1("dgfrb1","", 3.1027/(1+3.1027));
	
	// argus par
	RooRealVar arm0b1("arm0b1","", 3.87);
	RooRealVar arcb1("arcb1","",   0.503262);
	RooRealVar arpb1("arpb1","",   -4.91777);

	// double gaus pdf
	RooGaussian g1b1("g1b1","",ry,dgmub1,dgs1b1);  
	RooGaussian g2b1("g2b1","",ry,dgmub1,dgs2b1);  
	RooAddPdf   dgb1("dgb1","",RooArgList(g1b1,g2b1),dgfrb1);
	
	// argus pdf
	RooArgusBG argb1("argb1","",rx,arm0b1,arpb1,arcb1);
	
	//combined d-gaus + argus
	RooProdPdf pdfb1("pdfb1","",dgb1,argb1);

	// *******
	// bg hd
	// *******
	//Bkg HD Pol2    : -956.559 611.339 -96.9788 
	//Bkg HD Argus   : 3.87 0.265876 4.67281 
	
	// pol 2 par
	RooRealVar a0b2("a0b2","", -172.4);
	RooRealVar a1b2("a1b2","", 111.3);
	RooRealVar a2b2("a2b2","", -16.3);
   //1  a_{0}       -1.72365e+02   5.03044e+00   3.60058e-04  -2.96012e-10
   //2  a_{1}        1.11308e+02   1.91732e+00   1.16421e-04  -2.67747e-10
   //3  a_{2}       -1.63007e+01   5.25054e-01   3.75758e-05   9.45478e-10
	//RooRealVar a0b2("a0b2","", -956.559);
	//RooRealVar a1b2("a1b2","", 611.339);
	//RooRealVar a2b2("a2b2","", -96.9788);
	
	// argus par
	RooRealVar arm0b2("arm0b2","", 3.87);
	RooRealVar arcb2("arcb2","",   0.265876);
	RooRealVar arpb2("arpb2","",   4.67281);
	
	// pol 2 pdf
	RooPolynomial p2b2("p2b2","",ry, RooArgList(a0b2, a1b2, a2b2));
	
	// argus pdf
	RooArgusBG argb2("argb2","",rx,arm0b2,arpb2,arcb2);

	//combined pol2 + argus
	RooProdPdf pdfb2("pdfb2","",p2b2,argb2);


	// *******
	// create the background pdf
	// *******
	//relb=0.31731 relS=0.0235072
	//double rel  = XSb1*effb1/(XSb1*effb1+XSb2*effb2);
	RooRealVar relbg("relbg2d","",  (BkgXsNR*(effBnre[2]*BRe + effBnrm[2]*BRm))/((effBhde[2]+effBhdm[2])*BkgXs+BkgXsNR*(effBnre[2]*BRe + effBnrm[2]*BRm)));
	
	RooAddPdf  pdfbg1dS("pdfbgS","",RooArgList(argb1,argb2),relbg);
	RooAddPdf  pdfbg2d("pdfbg","",RooArgList(pdfb1,pdfb2),relbg);
	
	// *******
	// create total model
	// *******
	
	RooRealVar sigfrac("sigfrac","",0.1,-0.5,1);
	
	RooAddPdf  model2d ("model2d",  "2D model in m(ll) and m(ll)+m(pipi)", RooArgSet(pdfs, pdfbg2d),  sigfrac);  // pdf from full 2d models
	RooAddPdf  model1dS("model1dS", "1D model in m(ll)+m(pipi)"          , RooArgSet(nvs,  pdfbg1dS), sigfrac);  // pdf in m(ll)+m(pipi) from novosib. (S) and sum of 2 argus (B)
	RooAddPdf  model1dJ("model1dJ", "1D model in m(ll)"                  , RooArgSet(dgs,  p2b2),     sigfrac);  // pdf in m(ll) from d-gaus (S) and poly2 (B); NR bkg generated with signal pdf
	

// ------------------------
// END ROOFIT STUFF 
// ------------------------


typedef std::vector<TString> StrVec;


// -------------------------------------------
//  LINE SHAPES
// -------------------------------------------

// Masses [MeV] PDG2014

const double mp    = 938.272046;
const double mpi   = 139.57018;
const double mpi0  = 134.9766;
const double mrho  = 775.26;
const double momg  = 782.65;

const double mDp   = 1869.61;
const double mD0   = 1864.84;
const double mDps  = 2010.26;
const double mD0s  = 2006.96;

const double mJpsi = 3096.916;

// -------------------------------------------
// Widths [MeV] PDG2014

const double Grho = 149.1;
const double Gomg = 8.49;

// -------------------------------------------

const double delta = mDp + mDps - (mD0 + mD0s);
const double mu1 = mD0*mD0s/(mD0+mD0s);
const double mu2 = mDp*mDps/(mDp+mDps);
const double twopi = 2.*3.1415926;

double fr   = 0.00047;//0.007;
double fo   = 0.00271;//0.036;
double g12  = 0.137; 
double Gam0 = 1.0;
double B = 1;


// -------------------------------------------

// lookup table for the energy dependent Gammas
typedef std::map<int,double> LUT;

double LUTprec = 10000;
LUT LUT_J2pi, LUT_J3pi;
bool useLUT=false;

unsigned long lookups = 0;
unsigned long calcus  = 0;


TH2F *hbias=0;

// --------------------------------------------------------------------
//      END VARIABLE DEFINITIONS
// --------------------------------------------------------------------

void setPandaStyle(TString opt="")
{
	//TStyle *pandaStyle=new TStyle("PANDA","PANDA approved plots style");

	// use plain black on white colors
	gStyle->SetFrameBorderMode(0);
	gStyle->SetCanvasBorderMode(0);
	gStyle->SetPadBorderMode(0);
	gStyle->SetPadColor(0);
	gStyle->SetCanvasColor(0);
	gStyle->SetFrameFillColor(0);
	gStyle->SetFrameFillStyle(0);
	//gStyle->SetFillColor(0);
	
	// set the paper & margin sizes
	//double rightmargin = 0.05;
	//if (opt=="colz") rightmargin=0.12;

	gStyle->SetPaperSize(20,26);
	gStyle->SetPadTopMargin(0.05);
	gStyle->SetPadRightMargin(0.05);
	gStyle->SetPadBottomMargin(0.14);
	gStyle->SetPadLeftMargin(0.13);

	gStyle->SetStatColor(0);
	gStyle->SetOptStat(1110) ; 
	gStyle->SetStatBorderSize(1);
	gStyle->SetStatFont(42);
	gStyle->SetStatFontSize(0.04);
	gStyle->SetStatX(0.95);
	gStyle->SetStatY(0.95);
	gStyle->SetStatW(0.16);
	gStyle->SetStatH(0.18);


	// use sans serif fonts
// 	gStyle->SetTextFont(132);
// 	gStyle->SetTextSize(0.08);
	gStyle->SetLabelFont(42,"x");
	gStyle->SetLabelFont(42,"y");
	gStyle->SetLabelFont(42,"z");
	gStyle->SetLabelSize(0.05,"x");
	gStyle->SetTitleSize(0.06,"x");
	gStyle->SetTitleOffset(1.0,"X");
	gStyle->SetLabelSize(0.05,"y");
	gStyle->SetTitleSize(0.06,"y");
	gStyle->SetTitleOffset(1.,"y");
	gStyle->SetLabelSize(0.05,"z");
	gStyle->SetTitleSize(0.06,"z");
	
	gStyle->SetStripDecimals(false);

	// use bold lines and markers
	gStyle->SetMarkerStyle(20);
//	gStyle->SetHistLineWidth(1.85);
	//gStyle->SetLineStyleString(2,"[12 12]"); // postscript dashes

	// get rid of X error bars and y error bar caps
	gStyle->SetErrorX(0.001);
	
	gStyle->SetHistMinimumZero(1);
	
	// do not display any of the standard histogram decorations
	gStyle->SetOptTitle(0);
	gStyle->SetOptStat(1110);
	gStyle->SetOptFit(1011);
	gStyle->SetFitFormat(".3g");
		
	gStyle->SetPalette(1);

	// put tick marks on top and RHS of plots
 	gStyle->SetPadTickX(1);
 	gStyle->SetPadTickY(1);
	gStyle->SetTickLength(0.020,"xz");
	gStyle->SetTickLength(0.015,"y");
	gStyle->SetLineScalePS(3);

	//pandaStyle->SetLabelOffset(0.02,"xyz");

	//gROOT->Reset();
	//pandaStyle->SetLineWidth(2);
	//gROOT->SetStyle("PANDA");
	//gROOT->ForceStyle();
}

// --------------------------------------------------------------------
// breakup momentum of M -> m1 m2
double qbr(double M, double m1, double m2)
{
	return sqrt(( M*M - (m1+m2)*(m1+m2) ) * ( M*M - (m1-m2)*(m1-m2) ))/2./M;
}

// --------------------------------------------------------------------
// energy dependent width of X -> J/psi pi+ pi- (pi0)
double GamJnpi(double E, double f, double Gres, double mres, double mmin, int stps=300)
{	
	double M     = E + mD0 + mD0s;
	double mmax  = M - mJpsi;
	
	bool pipicase = mmin<(2*mpi+0.01);
	
	if (mmax<mmin) return 0.;
	
	int Ed =  (int)(E*LUTprec);  // discretised energy
	
	if (useLUT) 
	{
		if (pipicase) 
		{
			if (LUT_J2pi.find(Ed)!=LUT_J2pi.end()) {lookups++;return LUT_J2pi[Ed];}
		}
		else 
		{
			if (LUT_J3pi.find(Ed)!=LUT_J3pi.end()) {lookups++;return LUT_J3pi[Ed];}
		}
	}
	
	double mbin  = (mmax-mmin)/stps;
	
	double integ = 0.;
	
	for (double m=mmin; m<mmax; m+=mbin) integ += mbin*qbr(M,m,mJpsi)/((m-mres)*(m-mres) + 0.25*Gres*Gres);
	
	double res = f*Gres*integ/twopi;
	calcus++;
	
	if (useLUT) 
	{
		if (pipicase) LUT_J2pi[Ed] = res;
		else LUT_J3pi[Ed] = res;
	}
	
	return res;
}

// --------------------------------------------------------------------
// total energy dependent width
double GamJ(double E, double frho, double fomg)
{
	double res = GamJnpi(E,frho,Grho,mrho, 2.*mpi) + GamJnpi(E,fomg,Gomg,momg,2.*mpi+mpi0) + Gam0;
	return res;
}

// --------------------------------------------------------------------
// creates LUT to speed up Gamma calculation
void createLUTs(double Emin=-5, double Emax=20)
{
	LUT_J2pi.clear();
	LUT_J3pi.clear();
	
	double step = 1./LUTprec;
	
	for (double E = Emin; E<=Emax; E+=step)
	{
		double G2pi = GamJnpi(E, fr, Grho, mrho, 2*mpi       , 500);
		double G3pi = GamJnpi(E, fo, Gomg, momg, 2*mpi + mpi0, 500);
		
		LUT_J2pi[(int)(E*LUTprec)] = G2pi;
		LUT_J3pi[(int)(E*LUTprec)] = G3pi;
	}
	
	cout <<"Created LUT(J2pi) with "<<LUT_J2pi.size()<<" entries."<<endl;
	cout <<"Created LUT(J3pi) with "<<LUT_J3pi.size()<<" entries."<<endl;
	
	useLUT=true;
}


// --------------------------------------------------------------------
// auxilliary for line shape functions
double DE2(double E, double Ef, double frho, double fomg, double g)
{
	double res=0;
	
	if (E<0)
	{
		double kap1 = sqrt(-2.*mu1*E);
		double kap2 = sqrt(2.*mu2*(delta-E));	
		double GE   = GamJ(E,frho,fomg);
		
		res  = (E - Ef - 0.5*g*kap1 - 0.5*g*kap2)*(E - Ef - 0.5*g*kap1 - 0.5*g*kap2) + GE*GE/4.;
	}
	else if (E<delta)
	{
		double k1   = sqrt(2.*mu1*E);
		double kap2 = sqrt(2.*mu2*(delta-E));
		double GE   = GamJ(E,frho,fomg);
			
		res  = (E - Ef - 0.5*g*kap2)*(E - Ef - 0.5*g*kap2) + (g*k1 + GE)*(g*k1 + GE)/4.;
	}
	else
	{
		double k1   = sqrt(2.*mu1*E);
		double k2   = sqrt(2.*mu2*(E-delta));
		double GE   = GamJ(E,frho,fomg);
			
		res  = (E - Ef)*(E - Ef) + (g*k1 + g*k2 + GE)*(g*k1 + g*k2 + GE)/4.;
	}
	
	return res;
}

// --------------------------------------------------------------------
// lineshape X -> J/psi pi+ pi-
double fcn_jpsi2pi_hanhart(double* x, double* p)
{
	double E    = x[0];
	
	double Ef   = p[0];
	double frho = p[1];
	double fomg = p[2];
	double g    = p[3];
	double B    = p[4];
	 	
	return B*GamJnpi(E,frho,Grho,mrho, 2.*mpi)/(twopi*DE2(E,Ef,frho,fomg,g));
}

// --------------------------------------------------------------------
// lineshape X -> J/psi pi+ pi- pi0
double fcn_jpsi3pi_hanhart(double* x, double* p)
{
	double E    = x[0];
	
	double Ef   = p[0];
	double frho = p[1];
	double fomg = p[2];
	double g    = p[3];
	double B    = p[4];
	 	
	return B*GamJnpi(E,fomg,Gomg,momg, 2.*mpi+mpi0)/(twopi*DE2(E,Ef,frho,fomg,g));
}

// --------------------------------------------------------------------
// lineshape X -> D*0 D0
double fcn_dds_hanhart(double* x, double* p)
{
	double E    = x[0];

	if (E<0) return 0;
		
	double Ef   = p[0];
	double frho = p[1];
	double fomg = p[2];
	double g    = p[3];
	double B    = p[4];
	
	double k1   = sqrt(2.*mu1*E);
	 	
	return B*g*k1/(twopi*DE2(E,Ef,frho,fomg,g));
}

// --------------------------------------------------------------------
// the real part of the scattering length a
// see: http://dx.doi.org/10.1103/PhysRevD.76.034007, formula (37)
double fcn_Rea(double *x, double *p)
{
	double Ef   = x[0];
	
	double frho = p[0];
	double fomg = p[1];
	double g    = p[2];

	double G  = GamJ(0,frho, fomg);
	
	double aN = sqrt(2*mu2*delta) + 2*Ef/g;
	double aD = aN*aN + G*G/(g*g);
	
	return -aN/aD*197.3;
}

// --------------------------------------------------------------------
// the convoluted J/psi pi+ pi- lineshape
double fcn_conv_jpsi2pi_hanhart(double *x, double *p)
{
	double E      = x[0];
	
	double sigE   = p[0];
	double bg     = p[1];
	double shiftE = p[2];
	
	if (sigE<=0) 
	{
		x[0] += shiftE;
		double val = fcn_jpsi2pi_hanhart(x,&p[3]);
		x[0] = E;
		return val;
	}
	
	double np   = 51;
	
	double xmin  = E + shiftE - 5*sigE;
	double xmax  = E + shiftE + 5*sigE;
	
	if (abs(p[3]+8.56)<1.5 && xmin<0 && xmax>0) np=201;
	
	double xstp  = (xmax-xmin)/(np-1);
	
	double sum=0., wsum=0.;
	
	for (double Ex=xmin; Ex<=xmax; Ex+=xstp)
	{
		x[0] = Ex;
		double G = exp(-0.5*(E-shiftE-Ex)*(E-shiftE-Ex)/(sigE*sigE));
		sum  += G;
		wsum += G*fcn_jpsi2pi_hanhart(x,&p[3]); 
		
	}
	
	x[0] = E;
	
	return wsum/sum+p[1];		
}


// -------------------------------------------
//  END LINE SHAPES
// -------------------------------------------


// --------------------------------------------------------------------

int SplitString(TString s, TString delim, StrVec &toks)
{
	toks.clear();
	TObjArray *tok = s.Tokenize(delim);
	int N = tok->GetEntries();	
	for (int i=0;i<N;++i) 
	{
		TString token = (((TObjString*)tok->At(i))->String()).Strip(TString::kBoth);
		toks.push_back(token);
	}
	return toks.size();
}

// --------------------------------------------------------------------
// saves a canvas in various types
void saveCanvas(TCanvas *c, TString name, TString ftypes="gif,C,root,pdf")
{
	StrVec types;
	SplitString(ftypes,",",types);
	
	for (unsigned int i=0;i<types.size();++i)
		c->SaveAs(savepath+savepref+name+"."+types[i]);
}

// --------------------------------------------------------------------
// saves a canvas in various types in QA directory
void saveCanvasQA(TCanvas *c, TString name, TString ftypes="gif,C,root,pdf")
{
	StrVec types;
	SplitString(ftypes,",",types);
	
	for (unsigned int i=0;i<types.size();++i)
		c->SaveAs(savepathqa+savepref+name+"."+types[i]);
}

// --------------------------------------------------------------------

void confhist(TH1 *h, int col=1, TString tx="m [GeV/c^{2}]", TString ty="", double offy=1.05)
{
	double width = (h->GetXaxis()->GetXmax()-h->GetXaxis()->GetXmin())/h->GetNbinsX()*1000.;
	
	if (width<0.001) h->SetNdivisions(505);
	
	if (ty=="") ty=TString::Format("counts / %.1f MeV/c^{2}",width);
	
	h->SetXTitle(tx);
	h->SetYTitle(ty);
	h->SetLineColor(col);
	h->SetStats(0);
	h->SetLineWidth(2);
}

// --------------------------------------------------------------------

void confgraph(TGraph *g, int col=1, int sty=20, TString tit="", TString tity="", TString titx="#Gamma [keV]")
{
	g->SetMarkerSize(1.4);
	g->SetMarkerStyle(sty);
	g->SetMarkerColor(col);
	g->SetLineColor(col);
	g->SetLineWidth(2);
	
	g->GetHistogram()->SetTitle(tit);
	g->GetHistogram()->SetXTitle(titx);
	g->GetHistogram()->SetYTitle(tity);
}

// --------------------------------------------------------------------

void divCanvas(TCanvas *c, int n)
{
	int nx = sqrt((double)n), ny = nx;
	
	if (nx*ny < n) nx += 1;
	if (nx*ny < n) ny += 1;
	
	if (n==20) {nx=4; ny=5;}
	if (n==30) {nx=5; ny=6;}
	if (n==40) {nx=5; ny=8;}
 		
	c->Divide(nx,ny);
}

// --------------------------------------------------------------------

double writeLatex(double x, double y, TString txt, int col=0, double size=0.05, int font=42)
{
	TLatex lat(x,y,txt);
	lat.SetNDC();
	
	lat.SetTextColor(col);
	lat.SetTextFont(font);
	lat.SetTextSize(size);
	
	lat.DrawLatex(x,y,txt);
	
	return lat.GetYsize();
}

// --------------------------------------------------------------------

void drawLegend(TGraph *g1=0, TGraph *g2=0, TGraph *g3=0, double x1=0.6, double y1=0.78, double x2=0.95, double y2=0.95)
{
	TLegend *leg1 = new TLegend(x1,y1,x2,y2);
	if (g1!=0 && g1->GetN()>0) leg1->AddEntry(g1,g1->GetTitle(),"lp");
	if (g2!=0 && g2->GetN()>0) leg1->AddEntry(g2,g2->GetTitle(),"lp");
	if (g3!=0 && g3->GetN()>0) leg1->AddEntry(g3,g3->GetTitle(),"lp");
	leg1->Draw();
}

// --------------------------------------------------------------------

double getGraphMax(TGraph *g)
{
	if (g==0) return -1e15;
	return TMath::MaxElement(g->GetN(),g->GetY()); 
}

// --------------------------------------------------------------------

double getGraphMin(TGraph *g)
{
	if (g==0) return 1e15;
	return TMath::MinElement(g->GetN(),g->GetY()); 
}

// --------------------------------------------------------------------

double getGraphMaxX(TGraph *g)
{
	if (g==0) return 0;
	// to set good start parameters find E with highest observed number of events
	int idx = TMath::LocMax(g->GetN(),g->GetY());
	
	double maxX,dummy;
	g->GetPoint(idx,maxX,dummy);
	
	return maxX;
}

// --------------------------------------------------------------------

double getGraphMinX(TGraph *g)
{
	if (g==0) return 0;
	// to set good start parameters find E with highest observed number of events
	int idx = TMath::LocMin(g->GetN(),g->GetY());
	
	double minX,dummy;
	g->GetPoint(idx,minX,dummy);
	
	return minX;
}

// --------------------------------------------------------------------

void timePattern(double min, double max, int steps, std::vector<double> &E)
{
	E.clear();
	
	double step = (max-min)/(steps-1);
	for (int i=0;i<steps;++i) E.push_back(min+i*step);
}

// --------------------------------------------------------------------
// creates a graph with scanned number of events
/*
void getMLfitS(double S, double B, double &Smeas, double &dSmeasl, double &dSmeash, bool savePanel=false, TString tit="")
{
	// generate signal data
	roosigfrac.setVal(1.0);
	RooDataSet *sdata = 0;
	if (S>0) sdata = roomodel.generate(roox,S) ;
	
	// generate background data
	roosigfrac.setVal(0);
	RooDataSet *bdata = 0;
	if (B>0) bdata =roomodel.generate(roox,B) ;
	
	RooDataSet data("data","data", roox);
	if (sdata) data.append(*sdata);
	if (bdata) data.append(*bdata);
	
	// Fit model to data
	roosigfrac.setVal((S+0.001)/(S+B));
	roomodel.fitTo(data, Minos(roosigfrac), PrintLevel(-1000)) ;
	
	Smeas   = roosigfrac.getVal()*(S+B);
	//dSmeas = roosigfrac.getError()*(S+B);
	dSmeasl = fabs(roosigfrac.getAsymErrorLo()*(S+B)); 
	dSmeash = fabs(roosigfrac.getAsymErrorHi()*(S+B));
	
	//cout <<Smeas<<" + "<<dSmeash<<" - "<<dSmeasl<<endl;
	 
	if (savePanel)
	{
		RooPlot* xframe = roox.frame(Title(tit),Bins(20)) ;
		data.plotOn(xframe, MarkerSize(0.5),XErrorSize(0)) ;
		roomodel.plotOn(xframe);

		//// Overlay the background component of model with a dashed line
		roomodel.plotOn(xframe,Components(roobkg),LineStyle(kDashed),LineColor(kRed)) ;
		xframe->Draw();
	}
	
	if (sdata) delete sdata;
	if (bdata) delete bdata;
}
*/
// --------------------------------------------------------------------
// creates a graph with scanned number of events
void getMLfitS(double S, double Bhd, double Bnr, double &Sm, double &dSmeasl, double &dSmeash, bool savePanel=false, TString tit="")
{
	// total number of events to be generated
	double N = S + Bhd + Bnr;

	RooDataSet *sdata = 0, *b1data = 0, *b2data = 0;
	RooPlot    *xframe = 0, *yframe = 0;
	
	if (savePanel)
	{
		xframe = rx.frame(Title(tit),Bins(20)) ;
		yframe = ry.frame(Title(tit),Bins(20)) ;
	}
	
	double inifr = S/(S+Bhd+Bnr)*gRandom->Gaus(1,0.1);
	sigfrac.setVal(inifr);
	
	switch (AnaMode) 
	{
		// ******* ML analysis in m(ll)
		// RooAddPdf  model1dJ("model1dJ", "1D model in m(ll)", RooArgSet(dgs,  p2b2),     sigfrac);  // pdf in m(ll) from d-gaus (S) and poly2 (B); NR bkg generated with signal pdf
		case 1:
		{
			// generate signal and background data
			sdata  = dgs.generate(ry,   S+Bnr);  // both signal and NR bkg
			b1data = p2b2.generate(ry,  Bhd);    // hadronic bkg
			
			// merge datasets
			RooDataSet data("data","data", ry);
			data.append(*sdata); data.append(*b1data);
			
			// and fit model to it
			model1dJ.fitTo(data, Minos(sigfrac), PrintLevel(-1000));
			
			// if we shall print, plot data and model on frame
			if (yframe) 
			{
				data.plotOn(yframe, MarkerSize(0.5), XErrorSize(0)) ;
				model1dJ.plotOn(yframe);
				model1dJ.plotOn(yframe, Components(p2b2), LineStyle(kDashed),LineColor(kGray+2));
				yframe->Draw();
			}
		}	
		break;
			
		// ******* ML analysis in m(ll)+m(pipi)
		// RooAddPdf  model1dS("model1dS", "1D model in m(ll)+m(pipi)",         , RooArgSet(nvs,  pdfbg1dS), sigfrac);  // pdf in m(ll)+m(pipi) from novosib. (S) and sum of 2 argus (B)
		case 2:
		{
			// generate signal and background data
			sdata  = nvs.generate(rx,    S);    // both signal and NR bkg
			b1data = argb1.generate(rx,  Bnr);  // NR bkg
			b2data = argb2.generate(rx,  Bhd);  // hadronic bkg

			// merge datasets
			RooDataSet data("data","data", rx);
			data.append(*sdata); data.append(*b1data); data.append(*b2data);
			
			// and fit model to it
			model1dS.fitTo(data, Minos(sigfrac), PrintLevel(-1000));
			
			// if we shall print, plot data and model on frame
			if (xframe) 
			{
				data.plotOn(xframe, MarkerSize(0.5), XErrorSize(0)) ;
				
				model1dS.plotOn(xframe);
				model1dS.plotOn(xframe,Components(argb1),LineStyle(kDashed), LineColor(kMagenta));
				model1dS.plotOn(xframe,Components(argb2),LineStyle(kDashed), LineColor(kRed));
				model1dS.plotOn(xframe,Components(RooArgSet(argb1,argb2)), LineStyle(kDashed), LineColor(kGray+2));
				
				xframe->Draw();
			}
		}	
		break;
			
		// ******* 2D ML analysis in m(ll) vs. m(ll)+m(pipi)
		// RooAddPdf  model2d ("model2d",  "2D model in m(ll) and m(ll)+m(pipi)", RooArgSet(pdfs, pdfbg2d),  sigfrac);  // pdf from full 2d models
		case 3:
		{
			// generate signal and background data
			sdata  = pdfs.generate(RooArgSet(rx,ry), S);
			b1data = pdfb1.generate(RooArgSet(rx,ry),Bnr) ;
			b2data = pdfb2.generate(RooArgSet(rx,ry),Bhd) ;
			
			// merge datasets
			RooDataSet data("data","data", RooArgSet(rx,ry));
			data.append(*sdata); data.append(*b1data); data.append(*b2data);
			
			// and fit model to it
			model2d.fitTo(data, Minos(sigfrac), PrintLevel(-1000));
		
			// if we shall print, plot data and model on frame
			if (xframe) 
			{
				data.plotOn(xframe, MarkerSize(0.5), XErrorSize(0)) ;
				
				model2d.plotOn(xframe);
				model2d.plotOn(xframe,Components(pdfb1),LineStyle(kDashed), LineColor(kMagenta));
				model2d.plotOn(xframe,Components(pdfb2),LineStyle(kDashed), LineColor(kRed));
				model2d.plotOn(xframe,Components(RooArgSet(pdfb1,pdfb2)), LineStyle(kDashed), LineColor(kGray+2));
				
				xframe->Draw();
			}
		}
		break;
			
		default: break;
	}
	
	// compute the 
	Sm      = sigfrac.getVal()*N;
	dSmeasl = fabs(sigfrac.getAsymErrorLo()*N); 
	dSmeash = fabs(sigfrac.getAsymErrorHi()*N);
	
	// delete datasets
	if (sdata) delete sdata;
	if (b1data) delete b1data;
	if (b2data) delete b2data;
}
 
// --------------------------------------------------------------------
// creates a graph with scanned number of events

TGraphAsymmErrors* doScan(TF1 *fs, double intlumi, std::vector<double> nrgs, bool verbose=false, TCanvas *csave=0)
{	
	bool savePanel = (csave==0) ? false : true;

	int N = nrgs.size();
	
	// background levels are expected to be constant across energy range
	double BhdExp = intlumi * BkgXs   * effBhd;  // hadronic DPM background
	double BnrExp = intlumi * BkgXsNR * effBnr;  // NR J/psi pi pi background; eff including BRe/BRm
	
	// create a new TGraph
	TGraphAsymmErrors *g=new TGraphAsymmErrors(nrgs.size());

	// this is the energy shift related to the absolute positioning 
	// for the first energy point
	double Eshift = 0;
	if (dEBeamAbs>0) Eshift = gRandom->Gaus(0,dEBeamAbs);
	
	//cout <<Eshift<<" "<<flush;
	currEshift=Eshift;
	
	// energy error for relative step
	double errE   = sqrt(dEBeamRel*dEBeamRel); 
	
	// position of function maximum in scan region
	double fsMaxX = fs->GetMaximumX(nrgs[0], nrgs[N-1]);
	
	// our real starting scan point
	double E_real = nrgs[0] + Eshift + gRandom->Gaus(0,dEBeamRel);
	
	for (unsigned int i=0; i<nrgs.size(); ++i)
	{
		// energy we believe to set (we need later for inserting the point in graph)
		double E_ideal = nrgs[i]; 
		
		// we go one step ( = difference between intended beam energies) forward and make a relative positioning error 
		if (i>0) E_real += (nrgs[i]-nrgs[i-1]) + gRandom->Gaus(0,dEBeamRel);
		
		// signal cross section at E_real; effS includes BRe/BRm
		double SigExp = intlumi * fs->Eval(E_real) * BRX * effS;   
		
		double Nm=0, Bm=0., Sm=0, errNl=0., errNh=0., SigRnd, BhdRnd, BnrRnd;
		
		
		// which analysis mode?
		switch (AnaMode) 
		{
		// cut and count analysis
		case 0 : 
			
			//Bm = gRandom->Poisson(BhdExp) + gRandom->Poisson(BnrExp);
			//Sm = gRandom->Poisson(SigExp);
			
			
			
			Nm = gRandom->Poisson(SigExp+BhdExp+BnrExp);//Sm + Bm;
			
			errNh =      0.5*TMath::ChisquareQuantile(1-0.158655,2*(Nm+1)) - Nm;
			errNl = Nm - 0.5*TMath::ChisquareQuantile(  0.158655,2*Nm);
			
		    break;
			
		// ML analysis 
		case 1:
		case 2:
		case 3:
			// do we want to create and store an panel with all fits? 
			if (savePanel) csave->cd(i+1);
			
			// generate 
			SigRnd = gRandom->Poisson(SigExp);
			BhdRnd = gRandom->Poisson(BhdExp);
			BnrRnd = gRandom->Poisson(BnrExp);
			
			// extract Sm via maximum likelihood fit with RooFit
			getMLfitS(SigRnd, BhdRnd, BnrRnd, Nm, errNl, errNh, savePanel, Form("E = %.2f MeV",E_ideal+3872));
			
			// if we plot the panel, we update it as well
			if (savePanel) csave->Update();
			
			break;
			
		default:
			break;			
		}
				
		if (verbose) printf("Scan point %2d/%2d : E_id = %4.1f MeV  E_re = %4.1f MeV ; S_exp = %5.1f  S_meas = %5.1f ; B_exp = %5.1f  B_meas = %5.1f\n",
		                    i, N, E_ideal, E_real, SigExp, Sm, BhdExp + BnrExp, Bm);
				
		// put extracted signal and background (Bm should be 0 for ML modes) in scan graph
		g->SetPoint(i,E_ideal, Nm);
		
		// and the according errors
		g->SetPointError(i, sqrt(i)*errE, sqrt(i)*errE, errNl, errNh);
	}
	
	if (verbose) cout <<"Current energy displacement  = "<<Eshift<<endl;
	if (verbose) cout <<"Maximum signal cross section = "<<fs->GetMaximum(nrgs[0],nrgs[N])<<endl;
	
	return g;
}

// --------------------------------------------------------------------

double doVoigtStudy(int Ntoy, TGraphErrors *gp, TGraphErrors *gb, std::vector<double> VgtWid, int steps, double intlumi, double Eres, bool ctlplot=false, int nplotsdump=0)
{
	TCanvas *cc1 = 0;
	TCanvas *cc2 = 0;
	TCanvas *cc3 = 0;
	
	// these canvas are for control plots (current fit and cumulated fit results)
	if (ctlplot) 
	{
		cc1=new TCanvas("cc1","cc1",10,10,500,500);
		cc2=new TCanvas("cc2","cc2",520,10,500,500);
		cc3=new TCanvas("cc3","cc3",1020,10,500,500);
	}
	
	// number of parameter settings
	int Nwid = VgtWid.size();
	
	// the scan graph fit functions
	TF1 *fVoigt=new TF1("fVoigt","[0]*TMath::Voigt(x-[1],[2],[3])+[4]");
	fVoigt->SetParNames("A","m_{0} [MeV]","#sigma_{beam} [MeV]","#Gamma [MeV]","a");
	fVoigt->SetNpx(300);

	// control distribution of fit results	
	TH1F *hFitG=new TH1F("hFitG","",200,-150,150);
	TF1 fGaus("fGaus","gaus(0)");
	fGaus.SetParNames("A","bias [keV]","#sigma [keV]");
	
	TH1F *hFitProb=new TH1F("hFitProb","Probabilty",100,0,1.01);
	
	for (int i=0;i<Nwid;++i)
	{
		double Gam = VgtWid[i];

		// set scan window size
		double Nsig = 2.5;  // factor for Gamma for fit window width; is 5 x FWHM
		double FWHM = 0.5346*Gam+sqrt(0.2166*Gam*Gam + 5.55*Eres*Eres);
		double xmin = -Nsig*FWHM; //max(-1.0,-Nsig*(Gam/2.36+Eres));
		double xmax =  Nsig*FWHM; //min(1.0,Nsig*(Gam/2.36+Eres));
		
		// use the window given by the freq-sweep acceptance of dp/p = 1e-3 => dE = 1.678 MeV Window width @ 3.872 GeV
		if (freqFixWin)
		{
			xmin = -0.839;
			xmax = 0.839;
		}	
			
		// we scale the convoluted function in the way that
		// the uncovoluted function has max cross section of SigXs
		fVoigt->SetRange(xmin, xmax);
		
		fVoigt->SetParameters(1,0,0,Gam,0);
		double A = SigXs/fVoigt->GetMaximum(xmin,xmax);
		
		// needed to determine sigma and mean of fitted Gammas
		double sumG=0, sumG2=0, sigG=0, meanG=0;
		
		hFitG->Reset();
		hFitProb->Reset();
		
		int Ngood=0, Nloop=0;
		
		// create the scan point array
		std::vector<double> range;
		timePattern(xmin,xmax, steps, range);
		
		// we fix the beam resolution for the fit
		fVoigt->FixParameter(2,Eres);  

		int ndump = nplotsdump;

		// wait for Ntoy good fits 
		while (Ngood<Ntoy && Nloop++<Ntoy*6)
		{
			// save the scan panel
			TCanvas *cscan = 0;
			if (ndump>0 && AnaMode>0)
			{
				cscan = new TCanvas("cscan","cscan",500,5,1000-(steps-20)*10,1000);
				divCanvas(cscan,steps);
			}
			
			fVoigt->SetParameters(A,0,Eres,Gam,0);
			double currmax = fVoigt->GetMaximum(xmin,xmax);
			
			// ********* The actual scan experiment ************
			TGraphAsymmErrors *g = doScan(fVoigt, intlumi, range, false, cscan);
	
			if (ctlplot) cc1->cd();	

			// to set good start parameters find E with highest observed number of events
			double maxX = getGraphMaxX(g), maxY = getGraphMax(g);

			// control plot of scan graph
			if (getGraphMin(g)>0) g->SetMinimum(0);
			g->GetHistogram()->GetXaxis()->SetRangeUser(xmin,xmax);
			confgraph(g,1,20,Form("Scan %d x %.0fnb^{-1} / #Gamma = %.0fkeV",steps,intlumi,Gam*1000),"events","E - E_{0} [MeV]");
			g->SetMarkerSize(1.0);

			if (ctlplot) g->Draw("AP");
			
			double bglevel = 0.5*(g->GetY()[0]+g->GetY()[steps-1]);
			// initial parameter settings for fit 
			//                        amp           ,  m0 , sigE, Gamma,              bg level 
			fVoigt->SetParameters( A/currmax*(maxY-bglevel), maxX, Eres, Gam  , bglevel );
			
			//fVoigt->Draw("same");
			//if (ctlplot) cc1->Update();
			//sleep(1);
			
			// and fit  
			g->Fit(fVoigt,"ql");
			g->Fit(fVoigt,"qlm");
			
			
			// check for good fit
			double Afit = fVoigt->GetParameter(0);
			double Fmax = fVoigt->GetMaximum(xmin,xmax);
			double Gfit = fVoigt->GetParameter(3);
			double Mfit = fVoigt->GetParameter(1);
			double PFit = fVoigt->GetProb();
			
			// signal amp <=0 or maximum of fitted function deviates too much from observed maximum or Gamma<0 -> reject
			if (Afit<=0 /*|| Fmax>1.5*maxY || Fmax<0.7*maxY*/ || Gfit<=0.001 || Gfit>Gam*15. || fVoigt->GetProb()<0.001) {delete g; continue;}

			if (ctlplot) cc1->Update();

			if ((++Ngood%50)==0) cout <<Ngood<<" "<<flush;
			
			if (cc1!=0 && ndump-->0)
			{
				TString name = Form("ScanMode%d_Voigt_%dx%.0finvnb_G%03dkeV_dE%.0fkeV_%d", AnaMode, steps, intlumi, (int)(Gam*1000), Eres*1000, ndump);
				
				saveCanvasQA(cc1,name,"gif,C,pdf");
				
				if (cscan) 
				{
					saveCanvasQA(cscan,name+"_panel","gif,C,pdf,root");
					delete cscan;
				}
									
				cc1->cd();
			}
						
			delete g;

			// control histo showing the distribution of fit values
			hFitG->Fill((Gfit-Gam)*1000);
			hFitProb->Fill(PFit);
			
			sumG += Gfit; sumG2 += Gfit*Gfit; 
		}
		
		cout <<" - "<<Ngood<<"/"<<Nloop<<" fits accepted"<< endl;
				
		// determine the sensitivity and bias
		meanG = sumG/Ngood;
		sigG  = sqrt((sumG2 - sumG*sumG/Ngood)/(Ngood-1));
		
		//gp->SetPoint(i,Gam*1000,sigG/meanG*100);
		//gb->SetPoint(i,Gam*1000,(meanG-Gam)/Gam*100);
		
		// take values from histogram
		double meanGh = hFitG->GetMean()/1000+Gam;
		double sigGh  = hFitG->GetRMS()/1000;
		
		// compute error estimates of mean and RMS
		// dmean       = RMS/sqrt(N)
		// dRMS        = RMS/sqrt(2*(N-1));
		
		double prec  = sigGh/meanGh;
		// d(RMS/mean) = RMS/mean * sqrt(1/(4(N-1)) + RMS^2/(N*mean^2))
		double dprec = prec*sqrt(0.25/(Ngood-1) + sigGh*sigGh/(Ngood*meanGh*meanGh)); 

		gp->SetPoint     (i, Gam*1000,  prec*100);
		gp->SetPointError(i,        0, dprec*100);
		
		double bias  = meanGh/Gam - 1.;
		// d(mean/G0-1) = dmean/G0 = RMS/(sqrt(N)*G0)
		double dbias = sigGh/(sqrt(Ngood)*Gam);
		
		gb->SetPoint     (i,Gam*1000,  bias*100);
		gb->SetPointError(i,       0, dbias*100);

		//gp->SetPoint(i,Gam*1000,sigGh/meanGh*100);
		//gb->SetPoint(i,Gam*1000,(meanGh-Gam)/Gam*100);
		
		
				
		cout <<i<<" : Input Gamma = "<<VgtWid[i]<<"  Output Gamma = "<<meanG<<" +- "<<sigG<<"  histo: Gamma = "<<meanGh<<" +- "<<sigGh<<endl;
		
		if (ctlplot)
		{
			cc2->cd();
			hFitG->SetTitle(Form("%dx%.0f/nb - #Gamma=%dkeV - dE=%.0fkeV",steps,intlumi,(int)(Gam*1000),Eres*1000));
			hFitG->SetXTitle("#Gamma_{fit}#minus#Gamma_{0} [keV]");
			fGaus.SetParameters(hFitG->GetMaximum(),(meanG-Gam)*1000,sigG*1000);
			hFitG->Fit(&fGaus,"q","",max(-Gam*1000,(double)-500),Gam*1000);
			hFitG->Fit(&fGaus,"qm","",max(-Gam*1000,(double)-500),500);
			hFitG->DrawCopy();
			cc2->Update();

			cc3->cd();
			hFitProb->Draw();
			hFitProb->Fit("pol1","q");
			cc3->Update();
			
			if (nplotsdump>0) saveCanvasQA(cc2,Form("QA_Mode%d_Voigt_%dx%.0finvnb_G%03dkeV_dE%.0fkeV_%dtoy",AnaMode,steps,intlumi,(int)(Gam*1000),Eres*1000,Ntoy),"gif,C,pdf");
			
		}

	}	
	
	if (cc1) delete cc1;
	if (cc2) delete cc2;
	if (cc3) delete cc3;

	delete fVoigt;
	delete hFitG;
	delete hFitProb;
}

// --------------------------------------------------------------------

void doHHStudy(int Ntoy, TGraphErrors *gp, std::vector<double> vEf, int steps, double intlumi, double Eres, bool ctlplot=false, int nplotsdump=0)
{
	
	TCanvas *cc1 = 0;
	TCanvas *cc2 = 0;
	TCanvas *cc3 = 0;
	
	// these canvas are for control plots (current fit and cumulated fit results)
	if (ctlplot) 
	{
		cc1=new TCanvas("cc1","cc1",10,10,500,500);
		cc2=new TCanvas("cc2","cc2",520,10,500,500);
		cc3=new TCanvas("cc3","cc3",1020,10,500,500);
	}

	// number of parameter settings
	int NEf = vEf.size();
	
	// the scan graph fit functions
	TF1 *fH2 = new TF1("fH2",fcn_conv_jpsi2pi_hanhart,-3,10,8);
	fH2->SetParNames("#sigma_{beam} [MeV]","a_{bg}","dE_{pos} [MeV]","E_{f} [MeV]","f_{#rho}","f_{#omega}","g","B");
	fH2->SetNpx(200);
	
	// control distribution of fit results	
	TH1F *hFitG=new TH1F("hFitG","",200,-1.5,1.5);
	TF1 fGaus("fGaus","gaus(0)");
	fGaus.SetParNames("A","bias [MeV]","#sigma [MeV]");
	
	TH1F *hFitProb=new TH1F("hFitProb","Probabilty",100,0,1.01);

	for (int i=0;i<NEf;++i)
	{
		if (ctlplot) cc1->cd();
		
		double Ef = vEf[i];
		
		double xmin = -Eres*6-0.6;
		double xmax =  Eres*6+0.4;
		if (Ef>-8.56) xmax += abs(Ef+8.56)*2.5;
		if (Ef>-8.56) xmin -= abs(Ef+8.56);

		// use the window given by the freq-sweep acceptance of dp/p = 1e-3 => dE = 1.678 MeV Window width @ 3.872 GeV
		if (freqFixWin)
		{
			xmin = -0.839;
			xmax = 0.839;
		}	
			
		//double mean = (Ef+8.56)/6.4; // empircal shift to take into account position of maximum and tail
		//if (mean>0) mean*=2;
		
		//// set scan window size
		//double xmin = -1+mean;
		//double xmax =  1+mean;
		//if (Ef>-8.56) xmax += abs(Ef+8.56)*2;
		
		fH2->SetRange(xmin,xmax);
			
		// we scale the convoluted function in the way that
		// the uncovoluted function has max cross section of SigXs
		fH2->SetParameters(0,0,0,Ef,fr,fo,g12,1);
		double A = SigXs/fH2->GetMaximum(xmin,xmax);
		fH2->SetParameters(0, 0, 0, Ef, fr, fo, g12 ,A);
		
		// needed to determine sigma and mean of fitted Ef
		double sumG=0, sumG2=0, sigG=0, meanG=0;
		
		hFitG->Reset();
		hFitProb->Reset();
		
		int Ngood=0, Nloop=0;
		
		// create the scan point array
		std::vector<double> range;
		timePattern(xmin,xmax, steps, range);
		
		// we fix some params for the fit
		fH2->FixParameter(0,Eres);  
		fH2->FixParameter(4,fr);  
		fH2->FixParameter(5,fo);  
		fH2->FixParameter(6,g12);  
		//fH2->FixParameter(2,0.); // no eshift

		int ndump = nplotsdump;

		// wait for Ntoy good fits 
		while (Ngood<Ntoy && Nloop++<Ntoy*3)
		{
			// save the scan panel
			TCanvas *cscan = 0;
			if (ndump>0 && AnaMode>0)
			{
				cscan = new TCanvas("cscan","cscan",500,5,1000-(steps-20)*10,1000);
				divCanvas(cscan,steps);
			}
			
			//                  0    1  2   3   4    5   6   7 
			fH2->SetParameters(Eres, 0, 0, Ef, fr, fo, g12 ,A);
			double currmax = fH2->GetMaximum(xmin,xmax);
			
			// ********* The actual scan experiment ************
			//double tmp = dEBeamAbs; //switch off shift
			//dEBeamAbs = 0; 
			TGraphAsymmErrors *g = doScan(fH2, intlumi, range, false, cscan);
			//dEBeamAbs = tmp;

			if (ctlplot) cc1->cd();	

			// control plot of scan graph
			if (getGraphMin(g)>0) g->SetMinimum(0);
			g->GetHistogram()->GetXaxis()->SetRangeUser(xmin,xmax);
			confgraph(g,1,20,Form("Scan %d x %.0fnb^{-1} / E_{f} = %.2fMeV",steps,intlumi,Ef),"events","E - E_{0} [MeV]");
			g->SetMarkerSize(1.0);

			if (ctlplot) g->Draw("AP");
			
			// to set good start parameters find E with highest observed number of events
			//fH2->FixParameter(2,currEshift); // no eshift
			double maxX = currEshift*gRandom->Gaus(1,0.05);//getGraphMaxX(g);
			double maxY = getGraphMax(g);
	
			//double bglevel = 0.5*(g->GetY()[0]+g->GetY()[steps-1]);
			double bglevel = 0.25*(g->GetY()[0]+g->GetY()[1]+g->GetY()[steps-2]+g->GetY()[steps-1]);
			// initial parameter settings for fit 
			//                 sig  bg level, Eshift ,    Ef, (fix parms),      amp
			fH2->SetParameters(Eres, bglevel, -maxX   ,    Ef, fr, fo, g12, A/currmax*(maxY-bglevel));

			//if (RooFitMode) 
			//{
				//fH2->SetParameter(1,intlumi*BkgXsNR*(effbeNR*BRe+effbmNR*BRm));
				////fH2->SetParLimits(1,0, intlumi*BkgXs*effb*0.2);
			//}
			
			// set some parameter limits for more stable convergence
			//fH2->SetParLimits(2,-0.5,0.5);  
			fH2->SetParLimits(7, 0.5*A/currmax*(maxY-bglevel), 2.0*A/currmax*(maxY-bglevel));
			
			// and fit  
			g->Fit(fH2,"q","",xmin,xmax);
			g->Fit(fH2,"qm","",xmin,xmax);

			if (ctlplot) cc1->Update();

			// check for good fit
			double Afit  = fH2->GetParameter(7);
			double Fmax  = fH2->GetMaximum(xmin,xmax);
			double Effit = fH2->GetParameter(3);
			double PFit  = fH2->GetProb();
			
			// signal amp <=0 or maximum of fitted function deviates too much from observed maximum or Gamma<0 -> reject
			if (Afit<=0 /*|| Fmax>1.5*maxY || Fmax<0.7*maxY */|| Effit<-50 || Effit>0 || fH2->GetProb()<0.01) {delete g; continue;}
			
			//if ((++Ngood%10)==0) 
			cout <<Ngood++<<" "<<flush;
			
			if (cc1!=0 && (ndump--)>0)
			{
				TString name = Form("ScanMode%d_Mol_%dx%.0finvnb_Ef%.1fdMeV_dE%.0fkeV_%d", AnaMode, steps, intlumi, fabs(Ef), Eres*1000, ndump);
				
				saveCanvasQA(cc1,name,"gif,C,pdf");
				
				if (cscan) 
				{
					saveCanvasQA(cscan,name+"_panel","gif,C,pdf");
					delete cscan;
				}
			}
			
			delete g;
			
			// control histo showing the distribution of fit values
			hFitG->Fill(Effit-Ef);
			hFitProb->Fill(PFit);
			
			sumG += Effit; sumG2 += Effit*Effit; 
		}
		
		
		cout <<" - "<<Ntoy<<"/"<<Nloop<<" fits accepted"<< endl;
		
		// determine the sensitivity and bias
		meanG = sumG/Ngood;
		sigG  = sqrt((sumG2 - sumG*sumG/Ngood)/(Ngood-1));
				
		gp->SetPoint(i,Ef,meanG);
		gp->SetPointError(i,Ef,sigG);
				
		cout <<i<<" : Input Gamma = "<<vEf[i]<<"  Output Gamma = "<<meanG<<" +- "<<sigG<<endl;
		
		if (ctlplot)
		{
			cc2->cd();
			hFitG->SetTitle(Form("%dx%.0f/nb - Ef=%dMeV - dE=%.0fkeV",steps,intlumi,(int)Ef,Eres*1000));
			hFitG->SetXTitle("Ef_{fit}-Ef [MeV]");
			fGaus.SetParameters(hFitG->GetMaximum(),meanG-Ef,sigG);
			hFitG->Fit(&fGaus,"q");
			hFitG->DrawCopy();
			cc2->Update();
			
			cc3->cd();
			hFitProb->Draw();
			cc3->Update();
			
			if (nplotsdump>0) saveCanvasQA(cc2,Form("QA_Mode%d_Mol_%dx%.0finvnb_Ef%.1fMeV_dE%.0fkeV_%dtoy",AnaMode, steps,intlumi,fabs(Ef),Eres*1000,Ntoy),"gif,C");
		}
		
	}	
	
	if (cc1) delete cc1;
	if (cc2) delete cc2;
	if (cc3) delete cc3;
	delete fH2;
	delete hFitG;
	delete hFitProb;
}

// --------------------------------------------------------------------

TGraph* convertSens(TGraphErrors *ge, int col=1, int sty=20, TString htit="", TString tit="", TString tity="", TString titx="#Gamma [keV]")
{
	if (ge==0) return 0x0;
	
	int n=ge->GetN();
	
	if (n==0) return 0x0;
	
	double *x=ge->GetX();
	double *y=ge->GetY();
	double *ey=ge->GetEY();
	
	TGraph *g=new TGraph(n);
	
	TF1 f1("f1","TMath::Gaus(x,0,1,1)");
	
	for (int i=0;i<n;++i)
	{
		double dInSigma = fabs((fabs(y[i])-8.56))/ey[i];
		g->SetPoint(i,x[i],f1.Integral(dInSigma,10)*100);
	}
	
	g->SetTitle(tit);
	confgraph(g, col, sty, htit, tity, titx);
	
	g->SetMinimum(-1);
	g->SetMaximum(35);
	
	return g;
}
// --------------------------------------------------------------------

void printGraph(TGraph *g, TString pref="")
{
	if (g==0) return;
	int n=g->GetN();
	if (n==0) return;
	
	double *x = g->GetX();
	double *y = g->GetY();
	
	cout <<"int N_"<<pref<<" = "<<n<<";"<<endl;
	
	cout <<"double x_"<<pref<<"["<<n<<"] = { ";
	for (int i=0;i<n;++i) cout <<x[i]<<(i<(n-1)?", ":" };\n");
	
	cout <<"double y_"<<pref<<"["<<n<<"] = { ";
	for (int i=0;i<n;++i) cout <<y[i]<<(i<(n-1)?", ":" };\n");
	
	cout <<"TGraph *g_"<<pref<<" = new TGraph("<<n<<", x_"<<pref<<", y_"<<pref<<");"<<endl;
}

// --------------------------------------------------------------------

void printGraphErr(TGraphErrors *g, TString pref="")
{
	if (g==0) return;
	int n=g->GetN();
	if (n==0) return;
	
	double *x = g->GetX();
	double *y = g->GetY();
	double *ex = g->GetEX();
	double *ey = g->GetEY();
	
	cout <<"int N_"<<pref<<" = "<<n<<";"<<endl;
	
	cout <<"double x_"<<pref<<"["<<n<<"] = { ";
	for (int i=0;i<n;++i) cout <<x[i]<<(i<(n-1)?", ":" };\n");
	
	cout <<"double y_"<<pref<<"["<<n<<"] = { ";
	for (int i=0;i<n;++i) cout <<y[i]<<(i<(n-1)?", ":" };\n");
	
	cout <<"double ex_"<<pref<<"["<<n<<"] = { ";
	for (int i=0;i<n;++i) cout <<ex[i]<<(i<(n-1)?", ":" };\n");
	
	cout <<"double ey_"<<pref<<"["<<n<<"] = { ";
	for (int i=0;i<n;++i) cout <<ey[i]<<(i<(n-1)?", ":" };\n");
	
	cout <<"TGraphErrors *g_"<<pref<<" = new TGraphErrors("<<n<<", x_"<<pref<<", y_"<<pref<<", ex_"<<pref<<", ey_"<<pref<<");"<<endl;
}


// --------------------------------------------------------------------

void xscan_final2(double Xs=100, double dEHR=83.9, int steps=40, double days=2, int Ntoy=100, int npdump=1, int fitmode=1, TString opt="BW Mol", TString pref="TestScan_")
{
	double xmin = -3, xmax = 3;
	
	SigXs = Xs;
	dEcmHR = dEHR/1000.;
	
	gROOT->ProcessLine( "gErrorIgnoreLevel = 4001;");	
	gROOT->ProcessLine( "gPrintViaErrorHandler = kTRUE;");	
	
	setPandaStyle();
	
	bool BWmode  = opt.Contains("BW");
	bool Molmode = opt.Contains("Mol");
	
	bool HLmode = !opt.Contains("!HL");
	bool HRmode = !opt.Contains("!HR");
	bool HESRrmode = !opt.Contains("!HESRr");
	
	freqFixWin = opt.Contains("fixwin");
	
	savepref = pref;
	
	AnaMode = fitmode;
	
	if (AnaMode>3) AnaMode=3;
	
	//Adapt the J/psi window for the m(ll) ML fit
	if (AnaMode==1) ry.setRange(3.097-0.16, 3.097+0.16);
	
	// Set the efficiencies according to mode
	effS   = effSige[AnaMode]*BRe + effSigm[AnaMode]*BRm;
	effBhd = effBhde[AnaMode]     + effBhdm[AnaMode];
	effBnr = effBnre[AnaMode]*BRe + effBnrm[AnaMode]*BRm;
	
	if (RooMsgService::instance().numStreams()>0)
	{
		RooMsgService::instance().deleteStream(0);
		RooMsgService::instance().deleteStream(1);
	}
	
	gStyle->SetStatX(0.95);
	gStyle->SetStatY(0.95);
	gStyle->SetStatW(0.17);
	gStyle->SetStatH(0.17);

	useLUT=true;
	
	//createLUTs();
	
	//gROOT->LoadMacro("setPandaStyle.C");
	gRandom->SetSeed();
	//setPandaStyle();
	
	// ********************************
	// define the line shape functions 
	// ********************************
	
	// conventional voigt
	TF1 *fVoigt=new TF1("fVoigt","[0]*TMath::Voigt(x-[1],[2],[3])+[4]");
	fVoigt->SetRange(-3,3);
	fVoigt->SetNpx(300);
	fVoigt->SetLineColor(1);
	fVoigt->SetParNames("A","m_{0}","#sigma_{beam} [MeV]","#Gamma [MeV]","a");
	
	// the set of widths to be studied
	std::vector<double> VgtWid;
	VgtWid.push_back(0.05);
	VgtWid.push_back(0.07);
	VgtWid.push_back(0.1);
	VgtWid.push_back(0.13);
	VgtWid.push_back(0.18);
	VgtWid.push_back(0.25);
	VgtWid.push_back(0.5);
	//VgtWid.push_back(0.8);
	int Nwid = VgtWid.size();
	
	//VgtWid.push_back(0.8);
	
	// Hanhart J/psi pi+ pi- unconvoluted
	// Hanhart J/psi pi+ pi- for HL mode 
	TF1 *fH2 = new TF1("fH2",fcn_conv_jpsi2pi_hanhart,xmin,xmax,8);
	fH2->SetParNames("#sigma_{beam} [MeV]","a_{bg}","dE_{pos} [MeV]","E_{f} [MeV]","f_{#rho}","f_{#omega}","g","B");
	fH2->SetNpx(300);
	fH2->SetLineColor(2);
	
	// the set of Ef to be studies
	std::vector<double> TestEf;
	//TestEf.push_back(-25);
	//TestEf.push_back(-22);
	//TestEf.push_back(-20.5);	
	//TestEf.push_back(-19.25);
	//TestEf.push_back(-18.25);
	//TestEf.push_back(-17.5);
	//TestEf.push_back(-16);
	//TestEf.push_back(-12);
	
	//TestEf.push_back(-14);
	//TestEf.push_back(-12);
	//TestEf.push_back(-10);	
	//TestEf.push_back(-9);
	//TestEf.push_back(-8);
	//TestEf.push_back(-7);
	//TestEf.push_back(-6);
	//TestEf.push_back(-5);
	
	//TestEf.push_back(-14);
	TestEf.push_back(-10);	
	TestEf.push_back(-9);
	TestEf.push_back(-8.8);
	TestEf.push_back(-8.3);
	TestEf.push_back(-8);
	TestEf.push_back(-7.5);
	TestEf.push_back(-7);
	//TestEf.push_back(-5);
	int NEf = TestEf.size();
	
	double Efmin=TestEf[0], Efmax = TestEf[NEf-1];
	cout <<Efmin<<" ... "<<Efmax<<endl;
	
	// ********************************
	// canvas for plotting
	// ********************************
	
	TCanvas *c1=0, *c2=0, *c3=0, *c4=0, *c5=0, *c6=0, *c7=0;
	
	//c1=new TCanvas("c1","c1",10,10,1200,800);
	//c1->Divide(3,2);

	c2 =new TCanvas("c2","c2",600,600);
	
	if (BWmode)
	{
		c3 =new TCanvas("c3","c3",1200,600);
		c3->Divide(2,1);
	}
	
	if (Molmode)
	{
		c4 =new TCanvas("c4","c4",10,50,600,600);
		c5 =new TCanvas("c5","c5",620,50,600,600);
		c6 =new TCanvas("c6","c6",1230,50,600,600);
		c7 =new TCanvas("c7","c7",620,600,600,600);
	}
	

	//int cnt=1;

	//// ********************************
	//// draw Voigtian lineshapes
	//// ********************************
	
	//// background histogram for line shape
	//TH1F *hv=new TH1F("hv","",100,-1.5,1.5);
	//hv->SetStats(0);
	//hv->SetXTitle("E-E_{0} [MeV]");
	//hv->SetYTitle("#sigma_{X} [nb]");
	//hv->SetMaximum(SigXs*1.05);
	
	//TGraph *gMaxSigmaUnFoldBW=new TGraph(Nwid);
	//TGraph *gMaxSigmaFoldBWHL=new TGraph(Nwid);
	//TGraph *gMaxSigmaFoldBWHR=new TGraph(Nwid);
	//confgraph(gMaxSigmaUnFoldBW, 1, 22, "Maximum cross section Voigtian","#sigma_{max} [nb]" );
	//confgraph(gMaxSigmaFoldBWHR, 4, 20, "Maximum cross section Voigtian","#sigma_{max} [nb]" );
	//confgraph(gMaxSigmaFoldBWHL, 2, 21, "Maximum cross section Voigtian","#sigma_{max} [nb]" );
	//gMaxSigmaUnFoldBW->SetTitle("Unfolded");
	//gMaxSigmaFoldBWHR->SetTitle("HR mode");
	//gMaxSigmaFoldBWHL->SetTitle("HL mode");

	//for (int i=0;i<Nwid;++i)
	//{
		//double Wid = VgtWid[i];
		//c1->cd(cnt);
		
		//// configure histogram
		//hv->SetTitle(Form("#Gamma = %.0f keV",Wid*1000));
		//hv->DrawCopy();
		
		//// scale lineshape, so that maximum of UNFOLDED FUNCTION = signal cross section		
		//fVoigt->SetParameters(1,0,0,Wid,0);
		//double A = SigXs/fVoigt->GetMaximum(xmin,xmax);
		
		//fVoigt->SetParameter(0,A);
		//fVoigt->SetLineColor(1);
		//if (i==0 || i==Nwid-1) fVoigt->DrawCopy("same");
		//gMaxSigmaUnFoldBW->SetPoint(i,Wid*1000,fVoigt->GetMaximum(xmin,xmax));
		
		////cout <<"int = "<<fVoigt->Integral(xmin,xmax);
		
		//fVoigt->SetParameters(A,0,dEcmHL,Wid,0);
		//fVoigt->SetLineColor(2);
		//if (i==0 || i==Nwid-1) fVoigt->DrawCopy("same");
		//gMaxSigmaFoldBWHL->SetPoint(i,Wid*1000,fVoigt->GetMaximum(xmin,xmax));
		
		////cout <<" "<<fVoigt->Integral(xmin,xmax);

		//fVoigt->SetParameters(A,0,dEcmHR,Wid,0);
		//fVoigt->SetLineColor(4);
		//if (i==0 || i==Nwid-1) 
		//{
			//fVoigt->DrawCopy("same"); 
			//cnt++;
			//writeLatex(0.16,0.85,Form("#Gamma_{X}=%.0f keV",Wid*1000.),1,0.055); 
		//}
			
		//gMaxSigmaFoldBWHR->SetPoint(i,Wid*1000,fVoigt->GetMaximum(xmin,xmax));
		////cout <<" "<<fVoigt->Integral(xmin,xmax)<<endl;

		//// add legend
		//if (i==0)
		//{
			//double yt = 0.85, xr = 0.60;
				
			//double hgt = writeLatex(xr,yt,"Unfolded",1,0.045)*1.5;
			//writeLatex(xr,yt-1*hgt,"HR (#sigma=84 keV)",4,0.045);
			//writeLatex(xr,yt-2*hgt,"HL (#sigma=168 keV)",2,0.045);
		//}

		
		//c1->Update();
	//}
	
	//c1->cd(3);
	//gMaxSigmaUnFoldBW->GetHistogram()->SetMinimum(0);
	//gMaxSigmaUnFoldBW->GetHistogram()->SetMaximum(105);
	//gMaxSigmaUnFoldBW->GetHistogram()->SetXTitle("#Gamma [keV]");
	//gMaxSigmaUnFoldBW->GetHistogram()->SetYTitle("#sigma_{max} [nb]");
	//gMaxSigmaUnFoldBW->GetHistogram()->GetXaxis()->SetRangeUser(0,550);
	//gMaxSigmaUnFoldBW->Draw("APL");
	//gMaxSigmaFoldBWHR->Draw("PL same");
	//gMaxSigmaFoldBWHL->Draw("PL same");
	//drawLegend(gMaxSigmaUnFoldBW,gMaxSigmaFoldBWHR,gMaxSigmaFoldBWHL, 0.6,0.14,0.95,0.34);
	//cnt++;
	
	//// ********************************
	//// draw Hanhart lineshapes
	//// ********************************
	//// Ref: http://dx.doi.org/10.1103/PhysRevD.76.034007
	
	//// set range for parameter E_f
	
	//// the real part of the scattering length a
	// see Ref. formula (37)	
	/*	c2->cd(); gPad->SetGridy(); gPad->SetGridx();
	TF1 *fRea = new TF1("fRea",fcn_Rea,Efmin, Efmax,3);
	fRea->SetRange(-12,-5);
	fRea->SetParameters(fr,fo,g12);
	fRea->SetLineWidth(2);
	fRea->GetHistogram()->SetXTitle("E_{f} [MeV]");
	fRea->GetHistogram()->SetYTitle("Re(a) [fm]");
	
	fRea->Draw();
	cout <<"Thresh:"<<fRea->GetX(0,Efmin,Efmax)<<endl;
	
	return;
	*/
	////c2->SaveAs("fig/X_lineshape_Re_a2.gif");
	////c2->SaveAs("fig/X_lineshape_Re_a2.C");
	////c2->SaveAs("fig/X_lineshape_Re_a2.root");
		
	
	//// background histogram for line shape
	//TH1F *h=new TH1F("h","",100,-1.5,1.5);
	//confhist(h, 1, "E-E_{0} [MeV]", "#sigma_{X} [nb]");
	//h->SetMaximum(SigXs*1.05);

	//TGraph *gMaxSigmaUnFoldMol=new TGraph(NEf);
	//TGraph *gMaxSigmaFoldMolHL=new TGraph(NEf);
	//TGraph *gMaxSigmaFoldMolHR=new TGraph(NEf);
	//confgraph(gMaxSigmaUnFoldMol, 1, 22, "Maximum cross section molecule","#sigma_{max} [nb]" );
	//confgraph(gMaxSigmaFoldMolHR, 4, 20, "Maximum cross section molecule","#sigma_{max} [nb]" );
	//confgraph(gMaxSigmaFoldMolHL, 2, 21, "Maximum cross section molecule","#sigma_{max} [nb]" );
	//gMaxSigmaUnFoldMol->SetTitle("Unfolded");
	//gMaxSigmaFoldMolHR->SetTitle("HR mode");
	//gMaxSigmaFoldMolHL->SetTitle("HL mode");
		

	//// loop to draw line shapes
	//for (int i=0; i<NEf;++i)
	//{
		//double TEf =TestEf[i];
		
		//c1->cd(cnt);
		////gPad->SetLogy();

		//// configure histogram
		//h->SetTitle(Form("E_{f} = %.2f MeV",TEf));
		//h->DrawCopy();
				
		//// scale lineshape, so that maximum of UNFOLDED FUNCTION = signal cross section		
		//fH2->SetParameters(0,0,0,TEf,fr,fo,g12,B);
		//double fmax = fH2->GetMaximum(xmin,xmax);
		//double A = B/fmax*SigXs;
		
		//double shiftE=0.0;
		
		//// unfolded
		//fH2->SetParameters(0,0,shiftE,TEf,fr,fo,g12,A);             
		//fH2->SetLineColor(1); 
		//if (i==0 || i==3) fH2->DrawCopy("same");
		//gMaxSigmaUnFoldMol->SetPoint(i,TEf,fH2->GetMaximum(xmin,xmax));
		
		//// HL mode
		//fH2->SetParameters( dEcmHL, 0.,shiftE , TEf,fr,fo,g12,A); 
		//fH2->SetLineColor(2); 
		//if (i==0 || i==3) fH2->DrawCopy("same");
		//gMaxSigmaFoldMolHL->SetPoint(i,TEf,fH2->GetMaximum(xmin,xmax));
		
		//// HR mode
		//fH2->SetParameters( dEcmHR, 0., shiftE, TEf,fr,fo,g12,A); 
		//fH2->SetLineColor(4); 
		//if (i==0 || i==3) 
		//{
			//fH2->DrawCopy("same");
			//cnt++;
			//writeLatex(0.16,0.85,Form("E_{f}=%.1f MeV",TEf),1,0.055); 
		//}
		
		//gMaxSigmaFoldMolHR->SetPoint(i,TEf,fH2->GetMaximum(xmin,xmax));
		
		
		////TLegend *leg1 = new TLegend(0.6,0.7,0.96,0.94);
		////leg1->AddEntry(fH2,"unfolded","l");
		////leg1->AddEntry(fH2cHR,"dE = 67 keV","l");
		////leg1->AddEntry(fH2cHL,"dE = 167 keV","l");
		////leg1->Draw();
		
		//c1->Update();
	//}	
	//c1->cd(6);
	//gMaxSigmaUnFoldMol->GetHistogram()->SetMinimum(0);
	//gMaxSigmaUnFoldMol->GetHistogram()->SetMaximum(105);
	//gMaxSigmaUnFoldMol->GetHistogram()->SetXTitle("E_{f} [MeV]");
	//gMaxSigmaUnFoldMol->GetHistogram()->SetYTitle("#sigma_{max} [nb]");
	//gMaxSigmaUnFoldMol->Draw("APL");
	//gMaxSigmaFoldMolHR->Draw("PL same");
	//gMaxSigmaFoldMolHL->Draw("PL same");
	//drawLegend(gMaxSigmaUnFoldMol,gMaxSigmaFoldMolHR,gMaxSigmaFoldMolHL, 0.6,0.68,0.95,0.88);
	//cnt++;
	

	//c1->SaveAs(savepath+"X_lineshapes_BW_Molecule2.gif");
	//c1->SaveAs(savepath+"X_lineshapes_BW_Molecule2.C");
	
	//saveCanvas(c1,"X_lineshapes","pdf");
	
	//return;
	
	// ********************************
	// Molecule sensitivity study
	// ********************************
	
	if (Molmode)
	{
		// ************ preparations ************	
		double molthresh = -8.56516;//-18.75;
		
		double bakdEBeamAbs = dEBeamAbs; // switch off position error
		//dEBeamAbs = 0.;
		
		TLine l1,l2;
		l1.SetLineColor(4);
		l1.SetLineWidth(2);
		l1.SetLineStyle(3);
		
		l2.SetLineColor(6);
		l2.SetLineWidth(2);
		l2.SetLineStyle(7);
		
		double Xmin, Xmax, Ymin, Ymax;

		TGraphErrors       *gHLHprec    = new TGraphErrors(NEf); gHLHprec->SetTitle("HL mode");		
		TGraphErrors       *gHRHprec    = new TGraphErrors(NEf); gHRHprec->SetTitle("HR mode");			
		TGraphErrors       *gHESRrHprec = new TGraphErrors(NEf); gHESRrHprec->SetTitle("HESRr mode");		
			
		
	
		// *********************************
		// ************ HL mode ************	
	
		if (HLmode)
		{
			cout <<"Lineshape Study HL: "<<steps<<" x "<<days<<" days data taking"<<endl<<endl;
			
			doHHStudy(Ntoy, gHLHprec, TestEf, steps, days*IntLHL, dEcmHL, true, npdump); 
			confgraph(gHLHprec,    1, 20, Form("Sensitivity Ef (%d x %.0fd in HL)",    steps,days),"E_{f} [MeV]" ,"E_{f,0} [MeV]" );
			
			//for (int i=0;i<TestEf.size();++i)
			//{
				//gHLHprec->SetPoint(i,TestEf[i],TestEf[i]+gRandom->Gaus(0,0.5));
				//gHLHprec->SetPointError(i,TestEf[i],sqrt(fabs(TestEf[i])));
			//}
			
			cout <<"HL Graph info:"<<endl<<"-----------------------"<<endl;
			gHLHprec->Print();
			cout<<"-----------------------"<<endl;
			printGraphErr(gHLHprec,"hlmol");
			cout<<"-----------------------"<<endl;
			
			//  ************ draw ************
			
			c4->cd(); gPad->SetGridy(); gPad->SetGridx();
			gHLHprec->Draw();
			
			Xmin = gHLHprec->GetHistogram()->GetXaxis()->GetXmin();
			Xmax = gHLHprec->GetHistogram()->GetXaxis()->GetXmax();
			Ymin = gHLHprec->GetHistogram()->GetYaxis()->GetXmin();
			Ymax = gHLHprec->GetHistogram()->GetYaxis()->GetXmax();
			
			l2.DrawLine(molthresh,Ymin,molthresh,Ymax);
			l2.DrawLine(Xmin,molthresh,Xmax,molthresh);
			l1.DrawLine(Xmin,Xmin,Xmax,Xmax);
			
			saveCanvas(c4,Form("Sensitivity_Mol_Mode%d_%dx%.0fd_toy%d_HL", AnaMode,steps,days,Ntoy));
		}
		
		// *********************************
		// ************ HR mode ************	
		if (HRmode)
		{
			cout <<"Lineshape Study HR: "<<steps<<" x "<<days<<" days data taking"<<endl<<endl;
			
			if (HRmode) doHHStudy(Ntoy,gHRHprec,    TestEf, steps, days*IntLHR,    dEcmHR, true, npdump); 
			confgraph(gHRHprec,    1, 20, Form("Sensitivity Ef (%d x %.0fd in HR)",    steps,days),"E_{f} [MeV]" ,"E_{f,0} [MeV]" );
			
			cout <<"HR Graph info:"<<endl<<"-----------------------"<<endl;
			gHRHprec->Print();
			cout<<"-----------------------"<<endl;
			printGraphErr(gHRHprec,"hrmol");
			cout<<"-----------------------"<<endl;
			
			// ************
			
			c5->cd(); gPad->SetGridy(); gPad->SetGridx();
			gHRHprec->Draw();
			
			Xmin = gHRHprec->GetHistogram()->GetXaxis()->GetXmin();
			Xmax = gHRHprec->GetHistogram()->GetXaxis()->GetXmax();
			Ymin = gHRHprec->GetHistogram()->GetYaxis()->GetXmin();
			Ymax = gHRHprec->GetHistogram()->GetYaxis()->GetXmax();
			
			l2.DrawLine(molthresh,Ymin,molthresh,Ymax);
			l2.DrawLine(Xmin,molthresh,Xmax,molthresh);
			l1.DrawLine(Xmin,Xmin,Xmax,Xmax);

			saveCanvas(c5,Form("Sensitivity_Mol_Mode%d_%dx%.0fd_toy%d_HR", AnaMode, steps,days,Ntoy));
		}

	
		// *********************************
		// ************ HESRr mode ************	
		if (HESRrmode)
		{
			
			cout <<"Lineshape Study HESRr: "<<steps<<" x "<<days<<" days data taking"<<endl<<endl;
			
			if (HESRrmode) doHHStudy(Ntoy,gHESRrHprec, TestEf, steps, days*IntLHESRr, dEcmHR,  true, npdump); 
			confgraph(gHESRrHprec, 1, 20, Form("Sensitivity Ef (%d x %.0fd in HESRr)", steps,days),"E_{f} [MeV]" ,"E_{f,0} [MeV]" );
		
			cout <<"HESRr Graph info:"<<endl<<"-----------------------"<<endl;
			gHESRrHprec->Print();
			cout<<"-----------------------"<<endl;
			printGraphErr(gHESRrHprec,"hesrrmol");
			cout<<"-----------------------"<<endl;
			
			// ************
			
			c6->cd(); gPad->SetGridy(); gPad->SetGridx();
			gHESRrHprec->Draw();
			
			Xmin = gHESRrHprec->GetHistogram()->GetXaxis()->GetXmin();
			Xmax = gHESRrHprec->GetHistogram()->GetXaxis()->GetXmax();
			Ymin = gHESRrHprec->GetHistogram()->GetYaxis()->GetXmin();
			Ymax = gHESRrHprec->GetHistogram()->GetYaxis()->GetXmax();
			
			l2.DrawLine(molthresh,Ymin,molthresh,Ymax);
			l2.DrawLine(Xmin,molthresh,Xmax,molthresh);
			l1.DrawLine(Xmin,Xmin,Xmax,Xmax);

			saveCanvas(c6,Form("Sensitivity_Mol_Mode%d_%dx%.0fd_toy%d_HESRr", AnaMode, steps,days,Ntoy));
		}
		
		// ********************************
		// draw prob stuff
		// ********************************
		
		TGraph *ghl=0, *ghr=0, *ghesr=0; 
		
		if (HLmode)    ghl   = convertSens(gHLHprec,    2, 21, Form("Prob. mismatch (%d x %dd)",steps,(int)days), "HL mode",    "P [%]","E_{f,0} [MeV]");
		if (HRmode)    ghr   = convertSens(gHRHprec,    4, 22, Form("Prob. mismatch (%d x %dd)",steps,(int)days), "HR mode",    "P [%]","E_{f,0} [MeV]");
		if (HESRrmode) ghesr = convertSens(gHESRrHprec, 1, 20, Form("Prob. mismatch (%d x %dd)",steps,(int)days), "HESRr mode", "P [%]","E_{f,0} [MeV]");

		c7->cd(); gPad->SetGridy();
		TString opt="APL";
		
		if (ghl!=0 && ghl->GetN()>0)     { ghl->Draw(opt);   opt="PL same";}
		if (ghr!=0 && ghr->GetN()>0)     { ghr->Draw(opt);   opt="PL same";}
		if (ghesr!=0 && ghesr->GetN()>0) { ghesr->Draw(opt); opt="PL same";}
		drawLegend(ghesr, ghr, ghl);

		l2.DrawLine(molthresh,-1,molthresh,35);

		
		saveCanvas(c7,Form("Prob_Mol_Mode%d_mis_%dx%.0f_toy%d", AnaMode, steps,days,Ntoy));

		dEBeamAbs = bakdEBeamAbs; // restore beam positioning error
	}
	
	if (BWmode)
	{
				
		// ********************************
		// HL BW sensitivity study
		// ********************************
		cout <<endl<<"HL Study: "<<steps<<" x "<<days<<" days data taking"<<endl<<endl;
		
		TGraphErrors  *gHLVprec = new TGraphErrors();	gHLVprec->SetTitle("HL mode");	
		TGraphErrors  *gHLVbias = new TGraphErrors();	gHLVbias->SetTitle("HL mode");	
		
		if (HLmode) doVoigtStudy(Ntoy, gHLVprec, gHLVbias, VgtWid, steps, days*IntLHL, dEcmHL, true, npdump); 

		cout <<"HL Graph info:"<<endl<<"-----------------------"<<endl;
		cout <<"Precision [%]"<<endl;
		gHLVprec->Print();
		cout <<"Bias [%]"<<endl;
		gHLVbias->Print();
		
		cout<<"-----------------------"<<endl;
		printGraphErr(gHLVprec,"hlprec");
		printGraphErr(gHLVbias,"hlbias");
		cout<<"-----------------------"<<endl;
		
		confgraph(gHLVprec, 2, 21, Form("Sensitivity BW #Gamma (%d x %.0fd)", steps,days),"#Delta#Gamma_{meas}/#bar{#Gamma}_{meas} [%]" );
		confgraph(gHLVbias, 2, 21, Form("Bias BW #Gamma (%d x %.0fd)", steps,days)       ,"(#Gamma_{meas}-#Gamma_{0})/#Gamma_{0} [%]");
		
		// ********************************
		// HR BW sensitivity study
		// ********************************
		cout <<endl<<"HR Study: "<<steps<<" x "<<days<<" days data taking"<<endl<<endl;
		
		TGraphErrors  *gHRVprec = new TGraphErrors();	gHRVprec->SetTitle("HR mode");		
		TGraphErrors  *gHRVbias = new TGraphErrors();	gHRVbias->SetTitle("HR mode");		
		
		if (HRmode) doVoigtStudy(Ntoy, gHRVprec, gHRVbias, VgtWid, steps, days*IntLHR, dEcmHR, true, npdump); 

		cout <<"HR Graph info:"<<endl<<"-----------------------"<<endl;
		cout <<"Precision [%]"<<endl;
		gHRVprec->Print();
		cout <<"Bias [%]"<<endl;
		gHRVbias->Print();
		
		cout<<"-----------------------"<<endl;
		printGraphErr(gHRVprec,"hrprec");
		printGraphErr(gHRVbias,"hrbias");
		cout<<"-----------------------"<<endl;
		
		confgraph(gHRVprec, 4, 22, Form("Sensitivity BW #Gamma (%d x %.0fd)", steps,days), "#Delta#Gamma/#bar{#Gamma} [%]" );
		confgraph(gHRVbias, 4, 22, Form("Bias BW #Gamma (%d x %.0fd)", steps,days)       , "(#Gamma-#Gamma_{0})/#Gamma_{0} [%]");
		
		// ********************************
		// HESRr BW sensitivity study 
		// ********************************
		
		cout <<"HESRr Study: "<<steps<<" x "<<days<<" days data taking"<<endl<<endl;
		
		TGraphErrors       *gHESRrVprec = new TGraphErrors();	gHESRrVprec->SetTitle("HESRr mode");		
		TGraphErrors       *gHESRrVbias = new TGraphErrors();	gHESRrVbias->SetTitle("HESRr mode");		
		
		if (HESRrmode) doVoigtStudy(Ntoy, gHESRrVprec, gHESRrVbias, VgtWid, steps, days*IntLHESRr, dEcmHR,  true, npdump); 

		cout <<"HESRr Graph info:"<<endl<<"-----------------------"<<endl;
		cout <<"Precision [%]"<<endl;
		gHESRrVprec->Print();
		cout <<"Bias [%]"<<endl;
		gHESRrVbias->Print();
		
		cout<<"-----------------------"<<endl;
		printGraphErr(gHESRrVprec,"hesrrprec");
		printGraphErr(gHESRrVbias,"hesrrprec");
		cout<<"-----------------------"<<endl;
		
		confgraph(gHESRrVprec, 1, 20, Form("Sensitivity #Gamma (%d x %.0fd)", steps,days),"#Delta#Gamma / #bar{#Gamma} [%]" );
		confgraph(gHESRrVbias, 1, 20, Form("Bias #Gamma (%d x %.0fd)", steps,days)       ,"(#bar{#Gamma}#minus#Gamma_{0}) / #Gamma_{0} [%]");
		
		// ********************************
		// draw stuff
		// ********************************
		
		c3->cd(1); gPad->SetGridy();
		
		int gymin=0, gymax=80;
		
		gHLVprec->SetMinimum(gymin);    gHLVprec->SetMaximum(gymax);
		gHESRrVprec->SetMinimum(gymin); gHESRrVprec->SetMaximum(gymax);
		gHRVprec->SetMinimum(gymin);    gHRVprec->SetMaximum(gymax);
		
		opt="APL";
		if (gHLVprec->GetN()>0)    {gHLVprec->Draw(opt); opt="PL same";}
		if (gHESRrVprec->GetN()>0) {gHESRrVprec->Draw(opt); opt="PL same";}
		if (gHRVprec->GetN()>0)    {gHRVprec->Draw(opt); opt="PL same";}
		
		
		drawLegend(gHESRrVprec,gHRVprec , gHLVprec, 0.6);		
		c3->cd(2);gPad->SetGridy();
		
		
		gHLVbias->SetMinimum(-10);    gHLVbias->SetMaximum(30);
		gHRVbias->SetMinimum(-10);    gHRVbias->SetMaximum(30);
		gHESRrVbias->SetMinimum(-10); gHESRrVbias->SetMaximum(30);
		
		TString opt="APL";
		if (gHLVbias->GetN()>0)    {gHLVbias->Draw(opt); opt="PL same";}
		if (gHESRrVbias->GetN()>0) {gHESRrVbias->Draw(opt); opt="PL same";}
		if (gHRVbias->GetN()>0)    {gHRVbias->Draw(opt); opt="PL same";}
		
		drawLegend(gHESRrVbias,gHRVbias , gHLVbias, 0.6);
		
		saveCanvas(c3, Form("Sensitivity_BW_Mode%d_%dx%.0fd_toy%d", AnaMode, steps,days,Ntoy));
		
	}
}