//Macro to try and get a handle on making a likelihood ratio test statistic
//More realisitic idea now though
//Use totalHadrons and the xi distribition and try to combine into a single test statisitic ::squeeee::
// L = Pois(totalHadrons)xProd[pdf(xi)]_0^totalHadrons
//Ian Connelly
//24 April 2012

//Note that xi and projP are technically correlated.
//They are being treated as not correlated for now
//Need to investigate if I can generate a 2D pdf constrained between xi = 0.1->0.9

#include "TFile.h"
#include "RooRealVar.h"
#include "RooWorkspace.h"
#include "RooArgSet.h"
#include "RooPoisson.h"
#include "RooDataSet.h"
#include "RooKeysPdf.h"
//#include "RooFit.h"
#include "RooPlot.h"
#include "TMath.h"
#include "TH1.h"
#include "TString.h"
#include "RooHistPdf.h"

using namespace RooFit;
using std::cout;
using std::endl;
//using namespace RooStats;

int main(){
    
    //*******************//*******************//*******************//
    // When editing for different background comparison
    // change - root file name
    //        - workspace which is not w0
    //        - name of pdfs being Get'ted
    // NOT the poisson mean or pdf variable
    //        - name of root file mv'd in clusterRun.sh
    //*******************//*******************//*******************//
    
    
    
    //Remove the numeric integration stream
    //Note that RooMsgService::instance().Print() will show the active streams
    RooMsgService::instance().getStream(1).removeTopic(NumIntegration);
    
    TString path = "/home/connelly/scratch3/testarea/analysis/kernalestimation/output";
    TFile* f0 = new TFile(path+"/workspace0.root","OPEN");
    RooWorkspace* w0 = (RooWorkspace*)f0->Get("w"); // sig_pyth
    TFile* f1 = new TFile(path+"/workspace1.root","OPEN");
    RooWorkspace* w1 = (RooWorkspace*)f1->Get("w"); // sig_herw
    TFile* f2 = new TFile(path+"/workspace2.root","OPEN");
    RooWorkspace* w2 = (RooWorkspace*)f2->Get("w"); // bkg_zjj0
    TFile* f3 = new TFile(path+"/workspace3.root","OPEN");
    RooWorkspace* w3 = (RooWorkspace*)f3->Get("w"); // bkg_zjj1
    TFile* f4 = new TFile(path+"/workspace4.root","OPEN");
    RooWorkspace* w4 = (RooWorkspace*)f4->Get("w"); // bkg_tt
    
    cout << "=========>Workspaces loaded." << endl;
    
    //File to hold the teststat results
    TFile* f_teststat = new TFile("TEST.root","RECREATE");
    
    //Try for now just using a single signal and single background before looking into combining the backgrounds (and maybe signals)
    //Variables
    RooRealVar* totalHadrons = w0->var("totalHadrons");
    RooRealVar* xi = w0->var("xi");
    RooRealVar* nHad = w0->var("nHadrons");
    RooRealVar* projP = w0->var("hadronPlaneMomentum");
    RooRealVar* xi_s = new RooRealVar("xi_s","xi_s",0.1,0.9); //RooRealVar for generation in xi range
    RooRealVar* xis = new RooRealVar("xis","xis",0.1,0.9); // New variable after tests with the 2D pdf - not chaning the above one as it works for now
    //RooRealVar* xi_s = xi;
    //xi_s->setRange("xirange",0.1,0.9);
    //xi->setRange("xirange2",0.1,0.9);
    
    TH1D* h_teststatSignal = new TH1D("h_teststatSignal","Likelihood ratio under signal toy MC data",500,0,20);
    TH1D* h_teststatBkg    = new TH1D("h_teststatBkg","Likelihood ratio under background toy MC data",500,0,20);
    
    TH1D* h_teststatSignal2D = new TH1D("h_teststatSignal2D","Likelihood ratio under signal toy MC data with 2D pdfs",500,0,20);
    TH1D* h_teststatBkg2D    = new TH1D("h_teststatBkg2D","Likelihood ratio under background toy MC data with 2D pdfs",500,0,20);
    
    //Mean number of total hadrons per process to give an overview
    double mean0 = w0->data("data_sig_pyth_totalHadrons")->mean(*totalHadrons);
    double mean1 = w1->data("data_sig_herw_totalHadrons")->mean(*totalHadrons);
    double mean2 = w2->data("data_bkg_zjj0_totalHadrons")->mean(*totalHadrons);
    double mean3 = w3->data("data_bkg_zjj1_totalHadrons")->mean(*totalHadrons);
    double mean4 = w4->data("data_bkg_tt_totalHadrons")->mean(*totalHadrons);
    cout << mean0 << " " << mean1 << " " << mean2 << " " << mean3 << " " << mean4 << endl;
    //These are the mean values actually needed in the Poissons
    double mean0P1 = w0->data("data_sig_pyth_varExcP1")->mean(*nHad);
    double mean0P2 = w0->data("data_sig_pyth_varExcP2")->mean(*nHad);
    double mean0P3 = w0->data("data_sig_pyth_varExcP3")->mean(*nHad);
    //Will consider zjj0 as single background to begin with
    double mean2P1 = w4->data("data_bkg_tt_varExcP1")->mean(*nHad);
    double mean2P2 = w4->data("data_bkg_tt_varExcP2")->mean(*nHad);
    double mean2P3 = w4->data("data_bkg_tt_varExcP3")->mean(*nHad);
    
    RooRealVar poisMean0P1 ("poisMean0P1","Mean of sig Poisson - pyth - P1",mean0P1,"");
    RooRealVar poisMean0P2 ("poisMean0P2","Mean of sig Poisson - pyth - P2",mean0P2,"");
    RooRealVar poisMean0P3 ("poisMean0P3","Mean of sig Poisson - pyth - P3",mean0P3,"");
    RooRealVar poisMean2P1 ("poisMean2P1","Mean of bkg Poisson - tt - P1",mean2P1,"");
    RooRealVar poisMean2P2 ("poisMean2P2","Mean of bkg Poisson - tt - P2",mean2P2,"");
    RooRealVar poisMean2P3 ("poisMean2P3","Mean of bkg Poisson - tt - P3",mean2P2,"");
    
    //Components - Poisson Distributions
    RooPoisson* poissonSigPythP1 = new RooPoisson("poissonSigPythP1","Poisson distribution of P1 hadrons - pyth", *nHad, poisMean0P1);
    RooPoisson* poissonSigPythP2 = new RooPoisson("poissonSigPythP2","Poisson distribution of P2 hadrons - pyth", *nHad, poisMean0P2);
    RooPoisson* poissonSigPythP3 = new RooPoisson("poissonSigPythP3","Poisson distribution of P3 hadrons - pyth", *nHad, poisMean0P3);
    RooPoisson* poissonBkgZjj0P1 = new RooPoisson("poissonBkgZjj0P1","Poisson distribution of P1 hadrons - zjj1", *nHad, poisMean2P1);
    RooPoisson* poissonBkgZjj0P2 = new RooPoisson("poissonBkgZjj0P2","Poisson distribution of P2 hadrons - zjj1", *nHad, poisMean2P2);
    RooPoisson* poissonBkgZjj0P3 = new RooPoisson("poissonBkgZjj0P3","Poisson distribution of P3 hadrons - zjj1", *nHad, poisMean2P3);
    //Components - xi pdf
    RooKeysPdf* pdf0xiP1 = (RooKeysPdf*)w0->pdf("pdfkeys_sig_pyth_xi_P1");
    RooKeysPdf* pdf0xiP2 = (RooKeysPdf*)w0->pdf("pdfkeys_sig_pyth_xi_P2");
    RooKeysPdf* pdf0xiP3 = (RooKeysPdf*)w0->pdf("pdfkeys_sig_pyth_xi_P3");
    RooKeysPdf* pdf2xiP1 = (RooKeysPdf*)w4->pdf("pdfkeys_bkg_tt_xi_P1");
    RooKeysPdf* pdf2xiP2 = (RooKeysPdf*)w4->pdf("pdfkeys_bkg_tt_xi_P2");
    RooKeysPdf* pdf2xiP3 = (RooKeysPdf*)w4->pdf("pdfkeys_bkg_tt_xi_P3");
    //Components - projP pdf
    RooKeysPdf* pdf0projP1 = (RooKeysPdf*)w0->pdf("pdfkeys_sig_pyth_hpm_P1");
    RooKeysPdf* pdf0projP2 = (RooKeysPdf*)w0->pdf("pdfkeys_sig_pyth_hpm_P2");
    RooKeysPdf* pdf0projP3 = (RooKeysPdf*)w0->pdf("pdfkeys_sig_pyth_hpm_P3");
    RooKeysPdf* pdf2projP1 = (RooKeysPdf*)w4->pdf("pdfkeys_bkg_tt_hpm_P1");
    RooKeysPdf* pdf2projP2 = (RooKeysPdf*)w4->pdf("pdfkeys_bkg_tt_hpm_P2");
    RooKeysPdf* pdf2projP3 = (RooKeysPdf*)w4->pdf("pdfkeys_bkg_tt_hpm_P3");
    //*****NEW*****
    //2D RooHistPdf xi - 0.1->0.9, hpm - normal range
    // - Unfortunate syntax - two numbers next to each other 0 = signal, other 2 = 2Dimensional
    RooHistPdf* pdf02DP1 = (RooHistPdf*)w0->pdf("histpdf_sig_pyth_varIncP1");
    RooHistPdf* pdf02DP2 = (RooHistPdf*)w0->pdf("histpdf_sig_pyth_varIncP2");
    RooHistPdf* pdf02DP3 = (RooHistPdf*)w0->pdf("histpdf_sig_pyth_varIncP3");
    RooHistPdf* pdf22DP1 = (RooHistPdf*)w4->pdf("histpdf_bkg_tt_varIncP1");
    RooHistPdf* pdf22DP2 = (RooHistPdf*)w4->pdf("histpdf_bkg_tt_varIncP2");
    RooHistPdf* pdf22DP3 = (RooHistPdf*)w4->pdf("histpdf_bkg_tt_varIncP3");

    pdf02DP1->Print();
    
    cout << "=========>Signal and background pdfs are initalised" << endl;
    
    cout << "=========>Generation start" << endl;
    
    // Datasets to hold the generated toy data to calculate L/L with
    RooDataSet* toyHadP1  = new RooDataSet("toyHadP1","toyHadP1",RooArgSet(*nHad));
    RooDataSet* toyHadP2  = new RooDataSet("toyHadP2","toyHadP2",RooArgSet(*nHad));
    RooDataSet* toyHadP3  = new RooDataSet("toyHadP3","toyHadP3",RooArgSet(*nHad));
    RooDataSet* toyXiP1   = new RooDataSet("toyXiP1","toyXiP1",RooArgSet(*xi_s));
    RooDataSet* toyXiP2   = new RooDataSet("toyXiP2","toyXiP2",RooArgSet(*xi_s));
    RooDataSet* toyXiP3   = new RooDataSet("toyXiP3","toyXiP3",RooArgSet(*xi_s));
    RooDataSet* toyProjP1 = new RooDataSet("toyProjP1","toyProjP1",RooArgSet(*projP));
    RooDataSet* toyProjP2 = new RooDataSet("toyProjP1","toyProjP1",RooArgSet(*projP));
    RooDataSet* toyProjP3 = new RooDataSet("toyProjP1","toyProjP1",RooArgSet(*projP));
    //***NEW***
    RooDataSet* toy2DP1 = new RooDataSet("toy2DP1","toy2DP1",RooArgSet(*xis,*projP));
    RooDataSet* toy2DP2 = new RooDataSet("toy2DP2","toy2DP2",RooArgSet(*xis,*projP));
    RooDataSet* toy2DP3 = new RooDataSet("toy2DP3","toy2DP3",RooArgSet(*xis,*projP));

    
    // Generation loop
    
    //Generating signal - code 0
    //I will then have to repeat everything here but generating for background
    //Gen signal - assess the probabilities for sig versus bkg and ratio Likelihood
    //Gen bkg - do the same
    //Plot both
    
    int toyGen = 5000; // Number of toys to be generated
    
    //Impose a loop with if conditions for signal and background generation
    
    for(int i = 0; i < 2; i++){
        
        //I still want my ratio to have the form of t_s = L(s|s)/L(b|s) and t_b = L(s|b)/L(b|b)
        //Therefore the code which is in place will function correctly
        //It is just the generate of - nHadrons, xi, projP that need to be handled with care
        // ie with an if(i==) statement
        if(i == 0){
            cout << "========>Generation code 0" << endl;
            //Signal Generation - code 0
            toyHadP1 = poissonSigPythP1->generate(*nHad,toyGen);
            toyHadP2 = poissonSigPythP2->generate(*nHad,toyGen);
            toyHadP3 = poissonSigPythP3->generate(*nHad,toyGen);
        }
        if(i == 1){
            cout << "========>Generation code 1" << endl;
            //Background Generation - code 1
            //Note - is Poisson for which ever bkg event is selected and the mean is passed to the pdf
            toyHadP1 = poissonBkgZjj0P1->generate(*nHad,toyGen);
            toyHadP2 = poissonBkgZjj0P2->generate(*nHad,toyGen);
            toyHadP3 = poissonBkgZjj0P3->generate(*nHad,toyGen);
        }
        //Sets to hold variables etc
        //Sets to access the dataset
        RooArgSet* setHadP1;
        RooArgSet* setHadP2;
        RooArgSet* setHadP3;
        //****
        RooArgSet* setXiP1;
        RooArgSet* setXiP2;
        RooArgSet* setXiP3;
        RooArgSet* setProjP1;
        RooArgSet* setProjP2;
        RooArgSet* setProjP3;
        //Set for the pdf values
        RooArgSet* setPdfHadP1;
        RooArgSet* setPdfHadP2;
        RooArgSet* setPdfHadP3;
        //****
        RooArgSet* setPdfXiP1;
        RooArgSet* setPdfXiP2;
        RooArgSet* setPdfXiP3;
        RooArgSet* setPdfProjP1;
        RooArgSet* setPdfProjP2;
        RooArgSet* setPdfProjP3;
        //Set for the bkg pdf values
        RooArgSet* setBkgPdfHadP1;
        RooArgSet* setBkgPdfHadP2;
        RooArgSet* setBkgPdfHadP3;
        //***
        RooArgSet* setBkgPdfXiP1;
        RooArgSet* setBkgPdfXiP2;
        RooArgSet* setBkgPdfXiP3;
        RooArgSet* setBkgPdfProjP1;
        RooArgSet* setBkgPdfProjP2;
        RooArgSet* setBkgPdfProjP3;
        
        //***NEW***
        //Set to hold generated data
        RooArgSet* set2DP1;
        RooArgSet* set2DP2;
        RooArgSet* set2DP3;
        //Set for the 2D pdf values
        RooArgSet* setPdf2DP1;
        RooArgSet* setPdf2DP2;
        RooArgSet* setPdf2DP3;
        //Set for the bkg 2D
        RooArgSet* setBkgPdf2DP1;
        RooArgSet* setBkgPdf2DP2;
        RooArgSet* setBkgPdf2DP3;

        
        
        int genHadP1 = 0;
        int genHadP2 = 0;
        int genHadP3 = 0;
        double genXiP1 = 0;
        double genXiP2 = 0;
        double genXiP3 = 0;
        double genProjP1 = 0;
        double genProjP2 = 0;
        double genProjP3 = 0;
        
        double genProbHadP1 = 0;
        double genProbHadP2 = 0;
        double genProbHadP3 = 0;
        double genProbXiP1 = 0;
        double genProbXiP2 = 0;
        double genProbXiP3 = 0;
        double genProbProjP1 = 0;
        double genProbProjP2 = 0;
        double genProbProjP3 = 0;
        
        //---bkg---
        double genBkgProbHadP1 = 0;
        double genBkgProbHadP2 = 0;
        double genBkgProbHadP3 = 0;
        double genBkgProbXiP1 = 0;
        double genBkgProbXiP2 = 0;
        double genBkgProbXiP3 = 0;
        double genBkgProbProjP1 = 0;
        double genBkgProbProjP2 = 0;
        double genBkgProbProjP3 = 0;
        
        
        setPdfHadP1 = (RooArgSet*)poissonSigPythP1->getObservables(toyHadP1);
        setPdfHadP2 = (RooArgSet*)poissonSigPythP2->getObservables(toyHadP2);
        setPdfHadP3 = (RooArgSet*)poissonSigPythP3->getObservables(toyHadP3);
        setPdfHadP1->Print();
        //---bkg---
        setBkgPdfHadP1 = (RooArgSet*)poissonBkgZjj0P1->getObservables(toyHadP1);
        setBkgPdfHadP2 = (RooArgSet*)poissonBkgZjj0P2->getObservables(toyHadP2);
        setBkgPdfHadP3 = (RooArgSet*)poissonBkgZjj0P3->getObservables(toyHadP3);
        
        for(int n = 0; n < toyGen; n++){
            cout <<"=========>Event " << n << endl;
            //Need to get the get then number of hadrons generated
            setHadP1 = (RooArgSet*)toyHadP1->get(n);
            nHad = (RooRealVar*)setHadP1->find(nHad->GetName());
            genHadP1 = static_cast<int>(nHad->getVal());
            
            setHadP2 = (RooArgSet*)toyHadP2->get(n);
            nHad = (RooRealVar*)setHadP2->find(nHad->GetName());
            genHadP2 = static_cast<int>(nHad->getVal());
            
            setHadP3 = (RooArgSet*)toyHadP3->get(n);
            nHad = (RooRealVar*)setHadP3->find(nHad->GetName());
            genHadP3 = static_cast<int>(nHad->getVal());
            
            cout << genHadP1 << "  " << genHadP2 << "  " << genHadP3 << endl;
            
            //Generation of xi, projP
            //These take the number of hadrons value generated from the Poisson
            //The poisson number generated is dependent on signal or background
            //Hence accordingly the number is related to i == 0/1
            if(i == 0){
                cout << "=========>Code 0" << endl;
                toyXiP1 = pdf0xiP1->generate(*xi_s,genHadP1);
                toyProjP1 = pdf0projP1->generate(*projP,genHadP1);
                
                toyXiP2 = pdf0xiP2->generate(*xi_s,genHadP2);
                toyProjP2 = pdf0projP2->generate(*projP,genHadP2);
                
                toyXiP3 = pdf0xiP3->generate(*xi_s,genHadP3);
                toyProjP3 = pdf0projP3->generate(*projP,genHadP3);
                
                //***NEW***
                toy2DP1 = pdf02DP1->generate(RooArgSet(*xis,*projP),genHadP1);
                toy2DP1->Print();
                toy2DP2 = pdf02DP2->generate(RooArgSet(*xis,*projP),genHadP2);
                toy2DP3 = pdf02DP3->generate(RooArgSet(*xis,*projP),genHadP3);
            }
            if(i == 1){
                cout << "=========>Code 1" << endl;
                toyXiP1 = pdf2xiP1->generate(*xi_s,genHadP1);
                toyProjP1 = pdf2projP1->generate(*projP,genHadP1);
                
                toyXiP2 = pdf2xiP2->generate(*xi_s,genHadP2);
                toyProjP2 = pdf2projP2->generate(*projP,genHadP2);
                
                toyXiP3 = pdf2xiP3->generate(*xi_s,genHadP3);
                toyProjP3 = pdf2projP3->generate(*projP,genHadP3);
                
                //***NEW***
                toy2DP1 = pdf22DP1->generate(RooArgSet(*xis,*projP),genHadP1);
                toy2DP2 = pdf22DP2->generate(RooArgSet(*xis,*projP),genHadP2);
                toy2DP3 = pdf22DP3->generate(RooArgSet(*xis,*projP),genHadP3);
            }
            cout << "=========>Hadron xi,projP generated according to nHadrons generated in each plane" << endl;
            
            //It may be redundant creating three of each RooArgSets for each pdf used
            setPdfXiP1 = pdf0xiP1->getObservables(*xi);
            setPdfXiP2 = pdf0xiP2->getObservables(*xi);
            setPdfXiP3 = pdf0xiP3->getObservables(*xi);
            setPdfProjP1 = pdf0projP1->getObservables(*projP);
            setPdfProjP2 = pdf0projP2->getObservables(*projP);
            setPdfProjP3 = pdf0projP3->getObservables(*projP);
            //---bkg---
            setBkgPdfXiP1 = pdf2xiP1->getObservables(*xi);
            setBkgPdfXiP2 = pdf2xiP2->getObservables(*xi);
            setBkgPdfXiP3 = pdf2xiP3->getObservables(*xi);
            setBkgPdfProjP1 = pdf2projP1->getObservables(*projP);
            setBkgPdfProjP2 = pdf2projP2->getObservables(*projP);
            setBkgPdfProjP3 = pdf2projP3->getObservables(*projP);
            
            //***NEW***
            setPdf2DP1 = pdf02DP1->getObservables(RooArgSet(*xis,*projP));
            setPdf2DP2 = pdf02DP2->getObservables(RooArgSet(*xis,*projP));
            setPdf2DP3 = pdf02DP3->getObservables(RooArgSet(*xis,*projP));
            setBkgPdf2DP1 = pdf22DP1->getObservables(RooArgSet(*xis,*projP));
            setBkgPdf2DP2 = pdf22DP2->getObservables(RooArgSet(*xis,*projP));
            setBkgPdf2DP3 = pdf22DP3->getObservables(RooArgSet(*xis,*projP));
            

            //Check the dependancies
            //toyXiP1->Print();
            //setPdfXiP1->Print();
            //toyProjP1->Print();
            //setPdfProjP1->Print();
            
            //setBkgPdfXiP1->Print();
            //setBkgPdfProjP1->Print();
            //pdf2xiP1->Print();
            
            //Poisson probabilities for the number of hadrons - use the TMath Poisson function (as this is used to evaluate the RooPoisson anyway)
            //These probabilities are evaluated the same way whether the dataset is signal generated of background
            
            genProbHadP1 = TMath::Poisson(genHadP1,mean0P1);
            genProbHadP2 = TMath::Poisson(genHadP2,mean0P2);
            genProbHadP3 = TMath::Poisson(genHadP3,mean0P3);
            //Poisson probs for bkg poissons
            genBkgProbHadP1 = TMath::Poisson(genHadP1,mean2P1);
            genBkgProbHadP2 = TMath::Poisson(genHadP2,mean2P2);
            genBkgProbHadP3 = TMath::Poisson(genHadP3,mean2P3);
            
            
            cout << "=========>Poisson probabilities (signal)" << endl;
            cout << genProbHadP1 << "  -  " << genProbHadP2 << "  -  " << genProbHadP3 << endl;
            long double totPoisProb = genProbHadP1*genProbHadP2*genProbHadP3;
            cout << totPoisProb << endl;
            cout << "=========>Poisson probabilities (bkg)" << endl;
            cout << genBkgProbHadP1 << "  -  " << genBkgProbHadP2 << "  -  " << genBkgProbHadP3 << endl;
            long double totBkgPoisProb = genBkgProbHadP1*genBkgProbHadP2*genBkgProbHadP3;
            cout << totBkgPoisProb << endl;
            
            //Loop within toy datasets 
            //Plane1
            cout << "=========>Plane 1" << endl;
            long double totXiProbP1 = 1;
            long double totProjProbP1 = 1;
            long double totBkgXiProbP1 = 1;
            long double totBkgProjProbP1 = 1;
            //***NEW
            long double tot2DProbP1 = 1;
            long double tot2DBkgProbP1 =1;
            long double tot2DProbP2 = 1;
            long double tot2DBkgProbP2 =1; 
            long double tot2DProbP3 = 1;
            long double tot2DBkgProbP3 =1;
            
            for(int m = 0; m < genHadP1; m++){
                //---xi---(remember to set xi with the xi_s value to eval pdf)
                setXiP1 = (RooArgSet*)toyXiP1->get(m);
                xi_s = (RooRealVar*)setXiP1->find(xi_s->GetName());
                genXiP1 = xi_s->getVal();
                xi->setVal(genXiP1);
                //---projP--- (no need to set the value as I am not limiting my projP range)
                setProjP1 = (RooArgSet*)toyProjP1->get(m);
                projP = (RooRealVar*)setProjP1->find(projP->GetName());
                genProjP1 = projP->getVal();
                //---Probabilities---
                *setPdfXiP1 = *toyXiP1->get(m);
                genProbXiP1 =  pdf0xiP1->getVal(setPdfXiP1);
                *setPdfProjP1 = *toyProjP1->get(m);
                genProbProjP1 = pdf0projP1->getVal(setPdfProjP1);
                //pdf0xiP1->Print();
                cout << genXiP1 << "  -  " << genProjP1 << "  |  " << genProbXiP1 << "  -  " << genProbProjP1 << "  |  " ;
                //---Product of probabilities---
                totXiProbP1 *= genProbXiP1;
                totProjProbP1 *= genProbProjP1;
                //---Bkg stuff---
                setBkgPdfXiP1->setRealValue("xi",genXiP1);
                genBkgProbXiP1 = pdf2xiP1->getVal(setBkgPdfXiP1);
                setBkgPdfProjP1->setRealValue(projP->GetName(),genProjP1); //Use projP->GetName as it is actually called hadronPlaneMomentum  
                genBkgProbProjP1 = pdf2projP1->getVal(setBkgPdfProjP1);
                cout << genBkgProbXiP1 << "  -  " << genBkgProbProjP1 << endl;
                //---Product of bkg probabilities--
                totBkgXiProbP1 *= genBkgProbXiP1;
                totBkgProjProbP1 *= genBkgProbProjP1;
		//Delete pointers

            }
            
            //***NEW***
            for(int m = 0; m < genHadP1; m++){
                set2DP1 = (RooArgSet*)toy2DP1->get(m);
                //set2DP1->Print();
                xis = (RooRealVar*)set2DP1->find(xis->GetName());
                //xis->Print();
                projP = (RooRealVar*)set2DP1->find(projP->GetName());
                //projP->Print();
                genXiP1 = xis->getVal();
                xi->setVal(genXiP1);
                genProjP1 = projP->getVal();
                cout << "Values from 2D - " << genXiP1 << " -- " << genProjP1 << endl;
                //---Probabilities
                *setPdf2DP1 = *toy2DP1->get(m);
                setPdf2DP1->setRealValue("xi",genXiP1);
                genProbXiP1 = pdf02DP1->getVal(setPdf2DP1); // Ensures the xi being evaluated has correct 0.1-0.9 value
                cout << "Combined probability : " << genProbXiP1 << endl;
                //---Bkg
                *setBkgPdf2DP1 = *toy2DP1->get(m);
                setBkgPdf2DP1->setRealValue("xi",genXiP1);
                genBkgProbXiP1 = pdf22DP1->getVal(setBkgPdf2DP1);
                cout << "Combined prob (bkg)  : " << genBkgProbXiP1 << endl;
                //---Total Proabailities
                tot2DProbP1 *= genProbXiP1;
                tot2DBkgProbP1 *= genBkgProbXiP1;
                
            }
            cout << "2D Probs : " << tot2DProbP1 << " -- " << tot2DBkgProbP1 << endl;
                
                cout << totXiProbP1 << "  -  " << totProjProbP1 << "  |  " << totBkgXiProbP1 << "  -  " << totBkgProjProbP1 << endl;
            //Plane2
            cout << "=========>Plane 2" << endl;
            long double totXiProbP2 = 1;
            long double totProjProbP2 = 1;
            long double totBkgXiProbP2 = 1;
            long double totBkgProjProbP2 = 1;
            for(int m = 0; m < genHadP2; m++){
                //---xi---
                setXiP2 = (RooArgSet*)toyXiP2->get(m);
                xi_s = (RooRealVar*)setXiP2->find(xi_s->GetName());
                genXiP2 = xi_s->getVal();
                xi->setVal(genXiP2);
                //---projP---
                setProjP2 = (RooArgSet*)toyProjP2->get(m);
                projP = (RooRealVar*)setProjP2->find(projP->GetName());
                genProjP2 = projP->getVal();
                //---Probabilities---
                *setPdfXiP2 = *toyXiP2->get(m);
                genProbXiP2 = pdf0xiP2->getVal(setPdfXiP2);
                *setPdfProjP2 = *toyProjP2->get(m);
                genProbProjP2 = pdf0projP2->getVal(setPdfProjP2);
                cout << genXiP2 << "  -  " << genProjP2 << "  |  " << genProbXiP2 << "  -  " << genProbProjP2 << "  |  ";  
                //---Product of probabilities---
                totXiProbP2 *= genProbXiP2;
                totProjProbP2 *= genProbProjP2;
                //---Bkg stuff---
                setBkgPdfXiP2->setRealValue("xi",genXiP2);
                genBkgProbXiP2 = pdf2xiP2->getVal(setBkgPdfXiP2);
                setBkgPdfProjP2->setRealValue(projP->GetName(),genProjP2); //Use projP->GetName as it is actually called hadronPlaneMomentum  
                genBkgProbProjP2 = pdf2projP2->getVal(setBkgPdfProjP2);
                cout << genBkgProbXiP2 << "  -  " << genBkgProbProjP2 << endl;
                //---Product of bkg probabilities--
                totBkgXiProbP2 *= genBkgProbXiP2;
                totBkgProjProbP2 *= genBkgProbProjP2;
                //Delete pointers

            }
            
            //***NEW***
            for(int m = 0; m < genHadP2; m++){
                set2DP2 = (RooArgSet*)toy2DP2->get(m);
                //set2DP1->Print();
                xis = (RooRealVar*)set2DP2->find(xis->GetName());
                //xis->Print();
                projP = (RooRealVar*)set2DP2->find(projP->GetName());
                //projP->Print();
                genXiP2 = xis->getVal();
                xi->setVal(genXiP2);
                genProjP2 = projP->getVal();
                cout << "Values from 2D - " << genXiP2 << " -- " << genProjP2 << endl;
                //---Probabilities
                *setPdf2DP2 = *toy2DP2->get(m);
                setPdf2DP2->setRealValue("xi",genXiP2);
                genProbXiP2 = pdf02DP2->getVal(setPdf2DP2); // Ensures the xi being evaluated has correct 0.1-0.9 value
                cout << "Combined probability : " << genProbXiP2 << endl;
                //---Bkg
                *setBkgPdf2DP2 = *toy2DP2->get(m);
                setBkgPdf2DP2->setRealValue("xi",genXiP2);
                genBkgProbXiP2 = pdf22DP2->getVal(setBkgPdf2DP2);
                cout << "Combined prob (bkg)  : " << genBkgProbXiP2 << endl;
                //---Total Proabailities
                tot2DProbP2 *= genProbXiP2;
                tot2DBkgProbP2 *= genBkgProbXiP2;
                
            }
            cout << "2D Probs : " << tot2DProbP2 << " -- " << tot2DBkgProbP2 << endl;
            
            cout << totXiProbP2 << "  -  " << totProjProbP2 << "  |  " << totBkgXiProbP2 << "  -  " << totBkgProjProbP2 << endl;
            //Plane3
            cout << "=========>Plane 3" << endl;
            long double totXiProbP3 = 1;
            long double totProjProbP3 = 1;
            long double totBkgXiProbP3 = 1;
            long double totBkgProjProbP3 = 1;
            for(int m = 0; m < genHadP3; m++){
                //---xi---
                setXiP3 = (RooArgSet*)toyXiP3->get(m);
                xi_s = (RooRealVar*)setXiP3->find(xi_s->GetName());
                genXiP3 = xi_s->getVal();
                xi->setVal(genXiP3);
                //---projP---
                setProjP3 = (RooArgSet*)toyProjP3->get(m);
                projP = (RooRealVar*)setProjP3->find(projP->GetName());
                genProjP3 = projP->getVal();
                //---Probabilities---
                *setPdfXiP3 = *toyXiP3->get(m);
                genProbXiP3 = pdf0xiP3->getVal(setPdfXiP3);
                *setPdfProjP3 = *toyProjP3->get(m);
                genProbProjP3 = pdf0projP3->getVal(setPdfProjP3);
                cout << genXiP3 << "  -  " << genProjP3 << "  |  " << genProbXiP3 << "  -  " << genProbProjP3 << "  |  ";
                //---Product of probabilities---
                totXiProbP3 *= genProbXiP3;
                totProjProbP3 *= genProbProjP3;
                //---Bkg stuff---
                setBkgPdfXiP3->setRealValue("xi",genXiP3);
                genBkgProbXiP3 = pdf2xiP3->getVal(setBkgPdfXiP3);
                setBkgPdfProjP3->setRealValue(projP->GetName(),genProjP3); //Use projP->GetName as it is actually called hadronPlaneMomentum  
                genBkgProbProjP3 = pdf2projP3->getVal(setBkgPdfProjP3);
                cout << genBkgProbXiP3 << "  -  " << genBkgProbProjP3 << endl;
                //---Product of bkg probabilities--
                totBkgXiProbP3 *= genBkgProbXiP3;
                totBkgProjProbP3 *= genBkgProbProjP3;
		//Delete pointers

            }
            
            //***NEW***
            for(int m = 0; m < genHadP3; m++){
                set2DP3 = (RooArgSet*)toy2DP3->get(m);
                //set2DP1->Print();
                xis = (RooRealVar*)set2DP3->find(xis->GetName());
                //xis->Print();
                projP = (RooRealVar*)set2DP3->find(projP->GetName());
                //projP->Print();
                genXiP3 = xis->getVal();
                xi->setVal(genXiP3);
                genProjP3 = projP->getVal();
                cout << "Values from 2D - " << genXiP3 << " -- " << genProjP3 << endl;
                //---Probabilities
                *setPdf2DP3 = *toy2DP3->get(m);
                setPdf2DP3->setRealValue("xi",genXiP3);
                genProbXiP3 = pdf02DP3->getVal(setPdf2DP3); // Ensures the xi being evaluated has correct 0.1-0.9 value
                cout << "Combined probability : " << genProbXiP3 << endl;
                //---Bkg
                *setBkgPdf2DP3 = *toy2DP3->get(m);
                setBkgPdf2DP3->setRealValue("xi",genXiP3);
                genBkgProbXiP3 = pdf22DP3->getVal(setBkgPdf2DP3);
                cout << "Combined prob (bkg)  : " << genBkgProbXiP3 << endl;
                //---Total Proabailities
                tot2DProbP3 *= genProbXiP3;
                tot2DBkgProbP3 *= genBkgProbXiP3;
                
            }
            cout << "2D Probs : " << tot2DProbP3 << " -- " << tot2DBkgProbP3 << endl;
            
            
            cout << totXiProbP3 << "  -  " << totProjProbP3 << "  |  " << totBkgXiProbP3 << "  -  " << totBkgProjProbP3 << endl;
            long double totProb = totPoisProb*totXiProbP1*totXiProbP2*totXiProbP3;//*totProjProbP1*totProjProbP2*totProjProbP3;
            long double totBkgProb = totBkgPoisProb*totBkgXiProbP1*totBkgXiProbP2*totBkgXiProbP3;//*totBkgProjProbP1*totBkgProjProbP2*totBkgProjProbP3;
            cout << "====>Signal L under signal     : " <<totProb << endl;
            cout << "====>Background L under signal : " << totBkgProb << endl;
            cout << "====>L(s|s)/L(b|s)             : " << totProb/totBkgProb << endl;
            
            long double tot2DProb = totPoisProb*tot2DProbP1*tot2DProbP2*tot2DProbP3;
            long double tot2DBkgProb = totBkgPoisProb*tot2DBkgProbP1*tot2DBkgProbP2*tot2DBkgProbP3;
            cout << "====>Signal L under signal 2D    : " <<tot2DProb << endl;
            cout << "====>Background L under signal 2D: " << tot2DBkgProb << endl;
            cout << "====>L(s|s)/L(b|s) 2D            : " << tot2DProb/tot2DBkgProb << endl;
            
            //Fill relevant histogram according to whether the dataset is signal or bkg generated
            if(i == 0){
                h_teststatSignal->Fill(totProb/totBkgProb); 
                h_teststatSignal2D->Fill(tot2DProb/tot2DBkgProb);               

            }
            if(i == 1){
                h_teststatBkg->Fill(totProb/totBkgProb);
                h_teststatBkg2D->Fill(tot2DProb/tot2DBkgProb);

            }
        }
    }
    
    
    f_teststat->Write();
    f_teststat->Close();
    cout << "~~~Fin.~~~" << endl;
    return 0;
    
}