//vim: sw=4 ts=4 fdm=marker et:                                                                                                                  

// -----------------------------------------------                                                                                             
//       Author: Niladribihari SAHOO <nsahoo@cern.ch>                                                                                               
//       Created:   [2016-03-30 WED 10:11]                                                                                               
// -----------------------------------------------        

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sstream>
#include <sys/stat.h>
#include <math.h>
#include <string.h>
#include <getopt.h> // passing unix-like arguments
//#include <regex> //c++11 feature should be fine using gcc491.

#include <TSystem.h>
#include <TStyle.h>
#include "TMatrixDSym.h"
#include "RooMultiVarGaussian.h"
#include <TF1.h>
#include <TF2.h>
#include <TH1.h>
#include <TH2.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TMath.h>
#include <TMinuit.h>
#include <TRandom3.h>
#include <TFile.h>
#include <TPad.h> 
#include <TCanvas.h> 
#include <TChain.h> 
#include <TChainElement.h>
#include <TPaveText.h>
#include <TLatex.h>
#include <TLine.h>
#include <TString.h>
#include <TGraphAsymmErrors.h>
#include <TLorentzVector.h>
#include <TVirtualFitter.h>

#include <RooConstVar.h>
#include <RooRealVar.h>
#include <RooAbsPdf.h>
#include <RooAddPdf.h>
#include <RooGaussian.h>
#include <RooBifurGauss.h>
#include <RooChebychev.h> 
#include <RooGenericPdf.h> 
#include <RooExponential.h>
#include <RooBernstein.h>
#include <RooPolynomial.h>
#include <RooPoisson.h>
#include <RooExtendPdf.h>
#include <RooProdPdf.h>
#include <RooDataHist.h>
#include <RooDataSet.h>
#include <RooAbsData.h>
#include <RooFitResult.h>
#include <RooPlot.h>
#include <RooAddition.h>
#include <RooProduct.h>
#include <RooMinuit.h>
#include <RooWorkspace.h>
#include <RooRandom.h>

//#include <stdlib.h>     /* srand, rand */
//#include <time.h>       /* time */

#include <TRandom3.h>
#include "tools.h" 
#define PI 3.14159265358979

using namespace std; 
using namespace RooFit;

// Tags configration
bool is7TeVCheck = false; // Using 2011 efficiency map.
int isCDFcut = 0; // 0 for off, 1 for cdf, 2 for LHCb
TChain *ch=new TChain("tree");
TString plotpath="./plots";
//TString iwspacepath=".";
//TString owspacepath=".";
TString idatacardpath=".";
TString odatacardpath=".";
//Constants, Fit results for efficiency, etc.. //{{{
char genQ2range[11][32] = {"genQ2 < 2.00 && genQ2 > 1.00",
                           "genQ2 < 4.30 && genQ2 > 2.00",
                           "genQ2 < 8.68 && genQ2 > 4.30",
                           "genQ2 <10.09 && genQ2 > 8.68",
                           "genQ2 <12.86 && genQ2 >10.09",
                           "genQ2 <14.18 && genQ2 >12.86",
                           "genQ2 <16.00 && genQ2 >14.18",
                           "genQ2 <19.00 && genQ2 >16.00",
                           "genQ2 < 4.30 && genQ2 > 1.00",
                           "genQ2 <19.00 && genQ2 > 1.00",
                           "genQ2 < 6.00 && genQ2 > 1.00"};
char q2range[11][32] = {"Q2 < 2.00 && Q2 > 1.00",
                        "Q2 < 4.30 && Q2 > 2.00",
                        "Q2 < 8.68 && Q2 > 4.30",
                        "Q2 <10.09 && Q2 > 8.68",
                        "Q2 <12.86 && Q2 >10.09",
                        "Q2 <14.18 && Q2 >12.86",
                        "Q2 <16.00 && Q2 >14.18",
                        "Q2 <19.00 && Q2 >16.00",
                        "Q2 < 4.30 && Q2 > 1.00",
                        "Q2 <19.00 && Q2 > 1.00",
                        "Q2 < 6.00 && Q2 > 1.00"};
char q2rangeLatex[11][32] = {   " 1.00 < q^{2} < 2.00",
                                " 2.00 < q^{2} < 4.30",
                                " 4.30 < q^{2} < 8.68",
                                " 8.68 < q^{2} <10.09",
                                "10.09 < q^{2} <12.86",
                                "12.86 < q^{2} <14.18",
                                "14.18 < q^{2} <16.00",
                                "16.00 < q^{2} <19.00",
                                " 1.00 < q^{2} < 4.30",
                                " 1.00 < q^{2} <19.00",
                                " 1.00 < q^{2} < 6.00"};
double q2rangedn[11] = {1.00 , 2.00 , 4.30 , 8.68  , 10.09 , 12.86 , 14.18 , 16.00 , 1.00 ,  1.00 , 1.00};
double q2rangeup[11] = {2.00 , 4.30 , 8.68 , 10.09 , 12.86 , 14.18 , 16.00 , 19.00 , 4.30 , 19.00 , 6.00};
char mumuMassWindow[7][512] = { "Mumumass > 0",
                                "(Mumumass > 3.096916+3.5*Mumumasserr || Mumumass < 3.096916-5.5*Mumumasserr) && (Mumumass > 3.686109+3.5*Mumumasserr || Mumumass < 3.686109-3.5*Mumumasserr)",
                                "(Mumumass < 3.096916+3*Mumumasserr && Mumumass > 3.096916-5*Mumumasserr) || (Mumumass < 3.686109+3*Mumumasserr && Mumumass > 3.686109-3*Mumumasserr)",
                                //"Mumumass < 3.096916+3*Mumumasserr && Mumumass > 3.096916-5*Mumumasserr",
                                //"Mumumass < 3.686109+3*Mumumasserr && Mumumass > 3.686109-3*Mumumasserr",
                                "(Mumumass*Mumumass<8.68 && Bmass-Mumumass>2.182+0.16 && Bmass-Mumumass<2.182-0.16) || (Mumumass*Mumumass>10.09 && Mumumass*Mumumass<12.86 && Bmass-Mumumass>1.593+0.06 && Bmass-Mumumass<1.593-0.06) || (Mumumass*Mumumass>14.18 && Bmass-Mumumass>1.593+0.06 && Bmass-Mumumass<1.593-0.06)",
                                "(Mumumass*Mumumass < 8.68 && ( Bmass-Mumumass < 2.182+0.16 || Bmass-Mumumass > 2.182-0.16)) || (Mumumass*Mumumass>8.68 && Mumumass*Mumumass<10.09) || (Mumumass*Mumumass > 10.09 && Mumumass < 12.86 && ( Bmass-Mumumass < 1.593+0.06 || Bmass-Mumumass > 1.593-0.06)) || (Mumumass*Mumumass>12.86 && Mumumass*Mumumass<14.08) || (Mumumass*Mumumass > 14.08 && (Bmass-Mumumass > 1.593-0.06 || Bmass-Mumumass < 1.593+0.06))",
                                //"Bmass > 5.168 || (Bmass < 5.168 && Bmass-Mumumass-2.182 > 0.15) || (Bmass < 5.168 && Bmass-Mumumass-2.182 < -0.15) || (Bmass < 5.168 && Bmass-Mumumass-1.593 > 0.07) || (Bmass < 5.168 && Bmass-Mumumass-1.593 < -0.07)",
                                //"Bmass > 5.23 || (Bmass < 5.23 && Bmass-Mumumass-2.182 < 0.15 && Bmass-Mumumass-2.182 > -0.15) || (Bmass < 5.23 && Bmass-Mumumass-1.593 < 0.07 && Bmass-Mumumass-1.593 > -0.07)"
                                "Mumumasserr > 0",
                                "Mumumasserr > 0"
};//None, sig, bkg, #Jpsi, #Psi2S, CDF, anti-CDF, LHCb, anti-LHCb
double genAfb[9]={-0.160,-0.066,0.182,0.317,0.374,0.412,0.421,0.376, -0.098};
double genAfberr[9]={0.000434,0.000285,0.000223,0.000383,0.000277,0.000421,0.000395,0.000422,0.000333};
double genFl [9]={0.705,0.791,0.649,0.524,0.454,0.399,0.369,0.341, 0.7620};
double genFlerr[9]={0.000568,0.000409,0.000271,0.000420,0.000287,0.000415,0.000369,0.000361, 0.000240};
std::string f_accXrecoEff_ord0[9] = { // default values
    "11627.982364*((1.195422e-04*exp(-0.5*((CosThetaL-(-1.727343e-01))/2.021796e-01)**2)+1.156964e-04*exp(-0.5*((CosThetaL-(2.507083e-01))/2.478225e-01)**2)+4.629809e-05*exp(-0.5*((CosThetaL-(-5.148565e-01))/1.407258e-01)**2))*(7.165504e-05-2.621913e-05*CosThetaK+1.453609e-04*CosThetaK**2+2.274953e-05*CosThetaK**3-2.398253e-04*CosThetaK**4-4.428545e-05*CosThetaK**5+9.677067e-05*CosThetaK**6))",
    "9916.540629*((7.278431e-05*exp(-0.5*((CosThetaL-(-4.905860e-01))/1.878949e-01)**2)+7.448700e-05*exp(-0.5*((CosThetaL-(5.058518e-01))/2.003984e-01)**2)+1.425194e-04*exp(-0.5*((CosThetaL-(1.313125e-02))/2.957232e-01)**2))*(8.311598e-05-2.316101e-05*CosThetaK+1.476586e-04*CosThetaK**2-2.367362e-05*CosThetaK**3-1.683845e-04*CosThetaK**4+5.042865e-06*CosThetaK**5+1.243843e-05*CosThetaK**6))",
    "8.353802e+03*((1.344021e-04+1.980409e-05*CosThetaL+4.029664e-05*CosThetaL**2+1.560540e-05*CosThetaL**3-2.131400e-04*CosThetaL**4-3.310795e-05*CosThetaL**5+5.462426e-05*CosThetaL**6)*(1.136958e-04-3.718097e-05*CosThetaK+6.443598e-05*CosThetaK**2+5.683602e-05*CosThetaK**3-4.802073e-05*CosThetaK**4-7.413557e-05*CosThetaK**5-3.107261e-05*CosThetaK**6))",
    "4.392777e+04*((2.114489e-05+2.400662e-06*CosThetaL+2.759247e-05*CosThetaL**2+1.100568e-06*CosThetaL**3-4.538219e-05*CosThetaL**4-2.412249e-06*CosThetaL**5+5.307765e-06*CosThetaL**6)*(2.406814e-05-7.583489e-06*CosThetaK-9.968329e-06*CosThetaK**2+1.463576e-05*CosThetaK**3+3.247851e-05*CosThetaK**4-1.619795e-05*CosThetaK**5-2.949584e-05*CosThetaK**6))",
    "6.506619e+03*((1.349742e-04+1.528919e-05*CosThetaL+8.605597e-05*CosThetaL**2+1.312572e-05*CosThetaL**3-2.948919e-05*CosThetaL**4-9.566140e-06*CosThetaL**5-5.879247e-05*CosThetaL**6)*(1.581494e-04-3.384666e-05*CosThetaK-1.447583e-05*CosThetaK**2+3.758161e-05*CosThetaK**3+6.777260e-05*CosThetaK**4-5.585069e-05*CosThetaK**5-8.495213e-05*CosThetaK**6))",
    "4.625695e+04*((1.803216e-05+6.423635e-07*CosThetaL+9.704679e-06*CosThetaL**2+1.065779e-05*CosThetaL**3-1.658277e-06*CosThetaL**4-1.799046e-05*CosThetaL**5+6.089049e-06*CosThetaL**6)*(2.270524e-05-9.322913e-06*CosThetaK-1.587276e-05*CosThetaK**2+2.152708e-05*CosThetaK**3+5.615584e-05*CosThetaK**4-1.901528e-05*CosThetaK**5-4.887378e-05*CosThetaK**6))",
    "5.118383e+03*((1.668313e-04+1.911185e-05*CosThetaL+1.716389e-05*CosThetaL**2-3.192265e-05*CosThetaL**3+2.000329e-04*CosThetaL**4+1.783316e-05*CosThetaL**5-1.334724e-04*CosThetaL**6)*(2.056593e-04-4.151040e-05*CosThetaK-6.658669e-05*CosThetaK**2+3.742139e-05*CosThetaK**3+1.666491e-04*CosThetaK**4-5.072888e-05*CosThetaK**5-1.492963e-04*CosThetaK**6))",
    "3.837453e+03*((2.362599e-04-4.438020e-06*CosThetaL+3.318080e-05*CosThetaL**2+1.313482e-05*CosThetaL**3+7.878926e-05*CosThetaL**4-3.939653e-06*CosThetaL**5-2.211163e-05*CosThetaL**6)*(2.669904e-04-4.272653e-05*CosThetaK+1.487773e-05*CosThetaK**2+1.983652e-05*CosThetaK**3-9.317172e-05*CosThetaK**4-3.937610e-05*CosThetaK**5+4.831201e-05*CosThetaK**6))",
    "9916.540629*((7.278431e-05*exp(-0.5*((CosThetaL-(-4.905860e-01))/1.878949e-01)**2)+7.448700e-05*exp(-0.5*((CosThetaL-(5.058518e-01))/2.003984e-01)**2)+1.425194e-04*exp(-0.5*((CosThetaL-(1.313125e-02))/2.957232e-01)**2))*(8.311598e-05-2.316101e-05*CosThetaK+1.476586e-04*CosThetaK**2-2.367362e-05*CosThetaK**3-1.683845e-04*CosThetaK**4+5.042865e-06*CosThetaK**5+1.243843e-05*CosThetaK**6))"
};
// Lumi = Nreco/(cross section*branch factor*filter efficiency), cross section is 49.59e9 [pb] for 8TeV and 48.44e9 [pb] for 7TeV.
// BF_BuToK*MuMu = 1.07E-6, 1.12E-6(2014)
// BF_BuToK*Jpsi = 1.43E-3, 1.44E-3(2014)
// BF_BuToK*Psi2S = 6.7E-4, 6.7 E-4(2014)
// BF_JpsToMuMu  = 5.96E-2
// BF_Psi2sToMuMu= 6.70E-4
// BF_K*ToK0Pi  = 2/3  (K* decays to Kpi)
// BF_K0ToKs  = 1/2
// BF_KsToPiPi = 2/3
double datasetLumi[5] = {19.4,37378.629,295.761,218.472,9.81};//data, BuToKstarMuMu(16281.440+21097.189), BuToKstarJpsi(118.201+177.560), BuToKstarPsi2S(63.103,155.369), JpsiX
//}}}

double readParam(int iBin, const char parName[], int iColumn, double defVal=0., double forceReturn=999.)
{//{{{
    // Remark: first value is at iColumn=0.
    if (forceReturn != 999.) return forceReturn;

    std::vector<double> output;
    char lineBuff[1024];
    char *valBuff;
    memset(lineBuff,' ',1024*sizeof(char));
    FILE *fp = fopen(TString::Format("%s/fitParameters%d.txt",idatacardpath.Data(),iBin),"r");
    if (!fp){
        printf("WARNING: readParam, missing parameter files, by default return %f.\n",defVal);
        return defVal;
    }
    while(fgets(lineBuff,1024,fp) != NULL ){
        valBuff = strtok(lineBuff," ");
        if ( strcmp(valBuff,parName) == 0 ){
	  printf("INFO: readParam, matched %s!\n",valBuff);
            valBuff = strtok(NULL," ");
            while(valBuff != NULL){
                //output.push_back(stof(valBuff));//stof if c++11 function, use other function
                if (strcmp(valBuff,"nan") == 0 || strcmp(valBuff,"inf") == 0 ){
                    output.push_back(defVal);
                }else{
                    output.push_back(std::atof(valBuff));
                }
                valBuff = strtok(NULL," ");
            }
            break;
        }
        memset(lineBuff,' ',1024*sizeof(char));
    }
    fclose(fp);
    
    if (iColumn < output.size() ){
      printf("INFO: readParam, get %s[%d]=%e\n",parName,iColumn,output.at(iColumn));
        return output.at(iColumn);
    }else{
      printf("WARNING: readParam, empty column! Return %s[%d]=defVal=%f.\n",parName,iColumn,defVal);
        return defVal;
    }
}//}}}

std::string readParam(int iBin, const char parName[], string defVal="", string forceReturn="defaultForceReturn")
{//{{{
    // Remark: first value is at iColumn=0.
    if (forceReturn != "defaultForceReturn") return forceReturn;

    string output;
    char lineBuff[1024];
    char *valBuff;
    memset(lineBuff,' ',1024*sizeof(char));
    FILE *fp = fopen(TString::Format("%s/fitParameters%d.txt",idatacardpath.Data(),iBin),"r");
    if (!fp){
        printf("WARNING: readParam, missing parameter files, by default return %s.",defVal.c_str());
        return defVal;
    }
    while(fgets(lineBuff,1024,fp) != NULL ){
        valBuff = strtok(lineBuff," ");
        if ( strcmp(valBuff,parName) == 0 ){
            printf("INFO: readParam, matched %s!\n",valBuff);
            valBuff = strtok(NULL,"\n");
            output=string(valBuff);
            break;
        }
        memset(lineBuff,' ',1024*sizeof(char));
    }
    fclose(fp);
    
    if (output != ""){
        printf("INFO: readParam, get %s=%s\n",parName,output.c_str());
        return output;
    }else{
        printf("WARNING: readParam, empty item! Return %s=defVal=%s.\n",parName, defVal.c_str());
        return defVal;
    }
}//}}}
void writeParam(int iBin, const char parName[], double *val, int nVal=2, bool overwrite=true)
{//{{{
    if ( !overwrite ) return;

    struct stat fiBuff;
    FILE *fi = 0;
    if (stat(TString::Format("%s/fitParameters%d.txt",odatacardpath.Data(),iBin),&fiBuff) == 0){
        rename(TString::Format("%s/fitParameters%d.txt",odatacardpath.Data(),iBin),TString::Format("%s/fitParameters%d.txt.temp",odatacardpath.Data(),iBin));
        fi = fopen(TString::Format("%s/fitParameters%d.txt.temp",odatacardpath.Data(),iBin),"r");
    }else{
        fi = fopen(TString::Format("%s/fitParameters%d.txt.temp",odatacardpath.Data(),iBin),"w");
    }
    
    bool parExist = false;
    char lineBuff[1024];
    char *valBuff = 0;
    memset(lineBuff,' ',1024*sizeof(char));
    FILE *fp = fopen(TString::Format("%s/fitParameters%d.txt",odatacardpath.Data(),iBin),"w");
    while(fgets(lineBuff,1024,fi) != NULL ){
        valBuff = strtok(lineBuff," ");
        if ( strcmp(valBuff,parName) == 0 ){
            fprintf(fp,"%s",parName);
            int iVal = 0;
            while(iVal < nVal){
                fprintf(fp," %e",val[iVal]);
                iVal++;
            }
            fprintf(fp,"\n");
            parExist = true;
        }else{
            fprintf(fp,"%s",lineBuff);
            valBuff = strtok(NULL," ");
            while( valBuff != NULL ){
                fprintf(fp," %s",valBuff);
                valBuff = strtok(NULL," ");
            }
        }
        memset(lineBuff,' ',512*sizeof(char));
    }
    if (parExist == false){
        fprintf(fp,"%s",parName);
        int iVal = 0;
        while(iVal < nVal){
            fprintf(fp," %e",val[iVal]);
            iVal++;
        }
        fprintf(fp,"\n");
    }
    fclose(fp);
    fclose(fi);
    remove(TString::Format("%s/fitParameters%d.txt.temp",odatacardpath.Data(),iBin));
}//}}}
void writeParam(int iBin, const char parName[], string instring, bool overwrite=true)
{//{{{
    if ( !overwrite ) return;

    struct stat fiBuff;
    FILE *fi = 0;
    if (stat(TString::Format("%s/fitParameters%d.txt",odatacardpath.Data(),iBin),&fiBuff) == 0){
        rename(TString::Format("%s/fitParameters%d.txt",odatacardpath.Data(),iBin),TString::Format("%s/fitParameters%d.txt.temp",odatacardpath.Data(),iBin));
        fi = fopen(TString::Format("%s/fitParameters%d.txt.temp",odatacardpath.Data(),iBin),"r");
    }else{
        fi = fopen(TString::Format("%s/fitParameters%d.txt.temp",odatacardpath.Data(),iBin),"w");
    }
    
    bool parExist = false;
    char lineBuff[1024];
    char *valBuff = 0;
    memset(lineBuff,' ',1024*sizeof(char));
    FILE *fp = fopen(TString::Format("%s/fitParameters%d.txt",odatacardpath.Data(),iBin),"w");
    while(fgets(lineBuff,1024,fi) != NULL ){
        valBuff = strtok(lineBuff," ");
        if ( strcmp(valBuff,parName) == 0 ){
            fprintf(fp,"%s %s\n", parName, instring.c_str());
            parExist = true;
        }else{
            fprintf(fp,"%s",lineBuff);
            valBuff = strtok(NULL," ");
            while( valBuff != NULL ){
                fprintf(fp," %s",valBuff);
                valBuff = strtok(NULL," ");
            }
        }
        memset(lineBuff,' ',512*sizeof(char));
    }
    if (parExist == false){
        fprintf(fp,"%s %s\n", parName, instring.c_str());
    }
    fclose(fp);
    fclose(fi);
    remove(TString::Format("%s/fitParameters%d.txt.temp",odatacardpath.Data(),iBin));
    return;
}//}}}

// Physically allowed ranges from AN2014_129_v14, p25.
// Transformation rule comes from AN2014_129_v14, p28.
double toUnboundedFl(double fl){
    return TMath::Tan((fl-0.5)*TMath::Pi());
}
double toBoundedFl(double fl_ubd){
    return 0.5+TMath::ATan(fl_ubd)/TMath::Pi();
}
double toUnboundedAfb(double afb, double fl){
    return TMath::Tan(2./3.*afb*TMath::Pi()/(1-fl));
}
double toBoundedAfb(double afb_ubd, double fl_ubd){
    return 3./2.*(0.5-TMath::ATan(fl_ubd)/TMath::Pi())*TMath::ATan(afb_ubd)/TMath::Pi();
}
double toTransformedAs(double fs, double fl, double as)
{
    return as/(2*sqrt(3*fs*(1-fs)*fl)*0.89);
}
double toOriginAs(double as_tr, double fs, double fl_ubd)
{
    double fl=toBoundedFl(fl_ubd);
    return 2*sqrt(3*fs*(1-fs)*fl)*0.89*as_tr;
}

bool scanFlAfbPositivePdf(double afb, double fl)
{//{{{
    // Create test function to find possible domain for fl/afb. CosThetaL as x, CosThetaK as y.
    TString f2_format = "2.*[0]*y**2*(1.-x**2)+0.5*(1-[0])*(1.-y**2)*(1.+x**2)+4./3*[1]*(1-y**2)*x";
    TF2 *f2_model = new TF2("f2_model", f2_format.Data(),-1.,1.,-1.,1.);
    f2_model->SetTitle("PDF value;cos#theta_{L};cos#theta_{K}");

    // Find minimum position
    // partial_f2  / partial_x  = "-4.*[0]*y**2*x+(1-[0])*(1-y**2)*x+4./3*[1]*(1-y**2)"
    //      zero points:    x = 4*[1]*(1-y**2)/3/(4*[0]*y**2-(1-[0])*(1-y**2))
    // partial_f2  / partial_y  = "4*[0]*y*(1-x**2)-(1-[0])*y*(1.+x**2)-8./3*[1]*y*x" = "y*(4*[0]*(1-x**2)-(1-[0])*(1+x**2)-8/3*[1]*x)"
    //      zero points:    y = 0, x = (-b #pm sqrt(b**2-4*a*c))/2/a
    // Don't use f2_model->GetMinimumXY(xxMin,yyMin), the minimum cannot be found.
    f2_model->SetParameter(0,afb);
    f2_model->SetParameter(1,fl);
    bool isPositivePDF=true;
    if (f2_model->Eval(1.,0.) < 0 || f2_model->Eval(-1,0) < 0){
        isPositivePDF = false;
    }else{
        for (int i = 0; i < 100; ++i) {//cosThetaL
            for (int j = 0; j < 100; ++j) {//cosThetaK
                if (f2_model->Eval(0.02*i-1,0.02*j-1) < 0){
                    isPositivePDF = false;
                    break;
                }
            }
            if (!isPositivePDF) break;
        }
    }

    return isPositivePDF;
}//}}}

void scanFlAfbPositivePdf()
{//{{{
    // Create test function to find possible domain for fl/afb. CosThetaL as x, CosThetaK as y.
    TCanvas *canvas= new TCanvas("canvas");
    TH2F *h2_minPdfValue = new TH2F("h2_minPdfValue","",200,-1,1,200,0,1);
    h2_minPdfValue->SetStats(false);
    h2_minPdfValue->SetXTitle("A_{FB}");
    h2_minPdfValue->SetYTitle("F_{L}");

    for( int xBin = 1; xBin <= h2_minPdfValue->GetNbinsX(); xBin++){//afb
        for( int yBin = 1; yBin <= h2_minPdfValue->GetNbinsY(); yBin++){//fl
            bool isPositivePDF = scanFlAfbPositivePdf((yBin-0.5)/h2_minPdfValue->GetNbinsY(),(2.*xBin-1.)/h2_minPdfValue->GetNbinsX()-1);
            if(isPositivePDF){
                h2_minPdfValue->SetBinContent(xBin,yBin,1);
            }else{
            }
        }
    }

    // Draw contour
    h2_minPdfValue->Draw("COL");
    canvas->Update();
    canvas->SaveSource(TString::Format("%s/scanFlAfbPositivePdf.cc",plotpath.Data()));
    canvas->Print(TString::Format("%s/scanFlAfbPositivePdf.pdf",plotpath.Data()));
    return;
}//}}}



void printListOfTChainElements(TChain *chain){
    TObjArray *fileElements=chain->GetListOfFiles();
    int nFiles = fileElements->GetEntries();

    TIter next(fileElements);
    TChainElement *chEl=0;
    for( int entry=0; entry < nFiles; entry++ ) {
        chEl=(TChainElement*)next();
        printf("%s\n",chEl->GetTitle());
    }
    printf("DEBUG\t\t: %d files in the chain\n",nFiles);
}


void MultivariateGaussianTest(int iBin,  const char outfile[] = "accXrecoEff"){

  int dim = 7;
  RooArgList yVec;
  RooArgList muVec;
  int i,j;
  RooRealVar* y;
  RooRealVar* mu_y;
  for (i = 0; i < dim; i++) {
    char* name = Form("y%d", i);
    double mmean;
    mmean = readParam(iBin,"accXrecoEff2_ord0_L",i);
    double emean;
    emean = readParam(iBin,"accXrecoEff2Err_ord0_L",i);
    y = new RooRealVar(name, name, mmean-emean, mmean+emean);
    yVec.add(*y);
    char* mu_name = Form("mu_y%d",i);
    mu_y = new RooRealVar(mu_name, mu_name, readParam(iBin,"accXrecoEff2_ord0_L",i));
    muVec.add(*mu_y);
  }

  TMatrixDSym cov(dim);
  for (i = 0; i < dim; i++) {
    for (j = 0; j < dim; j++) {
      double errPar = readParam(iBin,TString::Format("covMat_L%d",i+1), j);
      if (! cov.IsValid()) continue;
      cov(i,j) = errPar;
      
    }
  }


  RooMultiVarGaussian mvg("mvg", "mvg", yVec, muVec, cov);
  RooDataSet* meee = mvg.generate(yVec, 10);
  meee->RooDataSet::write(TString::Format("eff_par_%d.txt",iBin));

}



void readNumbFromTxtFile(int iBin){

  int i;
  for (i=0; i<20; i++){
    double accXrecoEff    = readParam(iBin,"accXrecoEff2",i);
    double accXrecoEffErr = readParam(iBin,"accXrecoEff2Err",i);
    printf("acceptance efficiency par %d in bin %d = %.15f +/- %.15f \n", i, iBin, accXrecoEff, accXrecoEffErr);
    printf("----------------------------------------------------\n");
  }

}


void readEffParFromFile(int iBin){

  int i,j;
  for (i=0; i<20; i++){
    for (j=0; j<20; j++){
      //double EffPar    = readParam(iBin,"EffErr",i);
      //double EffParErr = readParam(iBin,"accXrecoEff2Err",i);
    double effPar = readParam(iBin,TString::Format("EffErr%d",i+1), j);
    printf("efficiency par %d in bin %d = %.15f  \n", i, iBin, effPar);
    printf("----------------------------------------------------\n");
    }
  }

}




int main(int argc, char** argv) {
  //      Tags                                                                                                                                         
  is7TeVCheck = false;
  //      Help message                                                                                                                 
  /*
  if (argc <= 2) {
    printf("./fit func iBin \n");
    printf("func: MultivariateGaussianTest readNumbFromTxtFile iBin \n");
    return 0;
  }
  */

  TString func    = argv[1];
  //TString infile  = argv[2];                                                                                                                      
  int iBin        = atoi(argv[2]);



  if (func == "MultivariateGaussianTest") {
    const char outfile[]="accXrecoEff";
    printf("\n===========================================================\n");
    printf("INFO:\t func used: %s \t Bin: %d \n", func.Data(), iBin);
    printf("============================================================\n");
    MultivariateGaussianTest(iBin, outfile);

  } else if (func == "readNumbFromTxtFile"){
    printf("\n============================================================\n");
    printf("INFO:\t func used: %s \t Bin: %d \n", func.Data(), iBin);
    printf("=============================================================\n");
    readNumbFromTxtFile(iBin);

  } else if (func == "readEffPar"){
    printf("\n============================================================\n");
    printf("INFO:\t func used: %s \t Bin: %d \n", func.Data(), iBin);
    printf("=============================================================\n");
    readEffParFromFile(iBin);

  } else {
    printf("No function available for: %s \n", func.Data());
    printf("./fit [func] [iBin] \n");
    printf("func: MultivariateGaussianTest readNumbFromTxtFile readEffPar \n");

  }

}