//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" using namespace RooFit; using std::cout; using std::endl; //using namespace RooStats; int main(){ //Remove the numeric integration stream //Note that RooMsgService::instance().Print() will show the active streams RooMsgService::instance().getStream(1).removeTopic(NumIntegration); TFile* f0 = new TFile("../workspace0.root","OPEN"); RooWorkspace* w0 = (RooWorkspace*)f0->Get("w"); // sig_pyth TFile* f1 = new TFile("../workspace1.root","OPEN"); RooWorkspace* w1 = (RooWorkspace*)f1->Get("w"); // sig_herw TFile* f2 = new TFile("../workspace2.root","OPEN"); RooWorkspace* w2 = (RooWorkspace*)f2->Get("w"); // bkg_zjj0 TFile* f3 = new TFile("../workspace3.root","OPEN"); RooWorkspace* w3 = (RooWorkspace*)f3->Get("w"); // bkg_zjj1 TFile* f4 = new TFile("../workspace4.root","OPEN"); RooWorkspace* w4 = (RooWorkspace*)f4->Get("w"); // bkg_tt cout << "=========>Workspaces loaded." << endl; //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 //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 = w2->data("data_bkg_zjj0_varExcP1")->mean(*nHad); double mean2P2 = w2->data("data_bkg_zjj0_varExcP2")->mean(*nHad); double mean2P3 = w2->data("data_bkg_zjj0_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 - zjj0 - P1",mean2P1,""); RooRealVar poisMean2P2 ("poisMean2P2","Mean of bkg Poisson - zjj0 - P2",mean2P2,""); RooRealVar poisMean2P3 ("poisMean2P3","Mean of bkg Poisson - zjj0 - 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 - zjj0", *nHad, poisMean2P1); RooPoisson* poissonBkgZjj0P2 = new RooPoisson("poissonBkgZjj0P2","Poisson distribution of P2 hadrons - zjj0", *nHad, poisMean2P2); RooPoisson* poissonBkgZjj0P3 = new RooPoisson("poissonBkgZjj0P3","Poisson distribution of P3 hadrons - zjj0", *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*)w2->pdf("pdfkeys_bkg_zjj0_xi_P1"); RooKeysPdf* pdf2xiP2 = (RooKeysPdf*)w2->pdf("pdfkeys_bkg_zjj0_xi_P2"); RooKeysPdf* pdf2xiP3 = (RooKeysPdf*)w2->pdf("pdfkeys_bkg_zjj0_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*)w2->pdf("pdfkeys_bkg_zjj0_hpm_P1"); RooKeysPdf* pdf2projP2 = (RooKeysPdf*)w2->pdf("pdfkeys_bkg_zjj0_hpm_P2"); RooKeysPdf* pdf2projP3 = (RooKeysPdf*)w2->pdf("pdfkeys_bkg_zjj0_hpm_P3"); cout << "=========>Signal and background pdfs are initalised" << endl; //Checking it works //RooPlot here was causing seg fault at end of program /* RooPlot* plot = nHad->frame(); poissonSigPythP1->plotOn(plot); plot->Draw(); */ //Components - Product per hadron on xi //---------------Testing-----------------// cout << "-----------------------------------Testing--------------------------------" << endl; //You have to access the variables inside the pdf AND the dataset -> They are not the same thing //So you need the other rooargset of observables from the pdf which also gets calculated //I think it loops in the same way as the data but will take the data value at that point and then return the pdf value at that point //Seems to work /* RooDataSet* d0 = (RooDataSet*)w0->data("data_sig_pyth_varIncP1"); //P1 xi,projP data RooDataSet* d0_ = (RooDataSet*)w0->data("data_sig_pyth_varExcP1"); //P1 nHad RooRealVar* nHadrons = (RooRealVar*)w0->var("nHadrons"); RooArgSet* set; //To hold data variable RooArgSet* set_; //---"-"--- RooArgSet* pdfxi = pdf0xiP1->getObservables(d0); //Needed to access the pdf observables pdfxi->Print(); double val = 0; //Init value for data double val_ = 0; int entries = d0->numEntries(); int entries_ = d0_->numEntries(); cout << "Number of events : " << entries_; cout << "Total number of tracks : " << entries << endl; int start = 0; //Counters to keep track of place through the xi dataset due to all lumped in together int end = 0; // ----"---"---- // Loop through the nHad dataset getting the number of entries to loop through the (xi,p) dataset with for(int j = 0; j < entries_ ; j++){ set_ = (RooArgSet*)d0_->get(j); //Get event nHadrons = (RooRealVar*)set_->find(nHadrons->GetName()); //Find the desired entry in dataset val_ = nHadrons->getVal(); //Get the value for the entry end += val_; long double prob = 0; //Init probability from pdf long double totprob = 1; // Not zero as going to do *= // What to do here for the case of zero hadrons -> Gets a probability of 1 // This is equivalent to it not contributing to the overall likelihood // Only alternative is setting equal to zero, but then that would set the entire likelihood to zero // So maybe leave this as it is. cout << "Start : " << start << " End : " << end << " Number of hadrons : " << val_ << endl; for(int i = start; i < end; i++){ *pdfxi = *d0->get(i); //Getting the entry with which to pass to the pdf set = (RooArgSet*)d0->get(i); //Getting xi value xi = (RooRealVar*)set->find(xi->GetName());//This line seems a bit redundant (as it seems to work without it) val = xi->getVal(); prob = pdf0xiP1->getVal(); //cout << val << " --- " << prob << endl; totprob = totprob*prob; } cout << j << " - " << "Total contribution to likelihood : " << totprob << endl; start = end; } //To check the number of events between 0.1 and 0.9 (may not work) //This does work //Need to use a new RooRealVar which has the reduced range and then generate using the pdf // and applying that RooRealVar xi->setRange("xirange",0.1,0.9); RooArgSet xiset("xi"); RooAbsReal* fracInt = pdf0xiP1->createIntegral(xiset,xiset,"xirange"); cout << fracInt->getVal() << endl; cout << d0->sumEntries(0,"xirange") << endl; RooRealVar xismall("xismall","xismall",0.1,0.9); RooDataSet* newdata = pdf0xiP1->generate(xismall,50); cout << newdata->numEntries() << endl; cout << xismall->Print() << endl; RooRealVar* newvar; //Variables and sets to hold the new variable RooArgSet* newset; for(int k = 0; k < 50; k++){ newset = (RooArgSet*)newdata->get(k); newdata->Print(); newset->Print(); newvar = (RooRealVar*)newset->find(xismall->GetName()); newvar->Print(); double value = newvar->getVal(); cout << value << endl; } //--------End of testing------------// cout << "----------------------------End of Testing--------------------------------" << endl; */ //-------Start of generation and L(s)/L(b) calculation--------// // Likelihood = Product(P1->P3)[Poisson(nHadrons)*Product(1->nHad)[pdf(xi|nHad)*pdf(projP|nHad)]] // To test - Generate signal events and then calculate Likelihood ratio between L(s)/L(b) // - Do this multiple times -> This will give how the variable varies under signal // - Then do the same for background // - All the pieces are effectively in place // Assume I have vectors? Might get issues with pdfs in vectors plus its a macro so won't do this 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)); // Generation loop //Generating signal - code 0 int toyGen = 50; // Number of toys to be generated toyHadP1 = poissonSigPythP1->generate(*nHad,toyGen); toyHadP2 = poissonSigPythP2->generate(*nHad,toyGen); toyHadP3 = poissonSigPythP3->generate(*nHad,toyGen); //Sets to hold variables etc RooArgSet* setHadP1; RooArgSet* setHadP2; RooArgSet* setHadP3; RooArgSet* setXiP1; RooArgSet* setXiP2; RooArgSet* setXiP3; RooArgSet* setProjP1; RooArgSet* setProjP2; RooArgSet* setProjP3; RooArgSet* setPdfHadP1; RooArgSet* setPdfHadP2; RooArgSet* setPdfHadP3; RooArgSet* setPdfXiP1; RooArgSet* setPdfXiP2; RooArgSet* setPdfXiP3; RooArgSet* setPdfProjP1; RooArgSet* setPdfProjP2; RooArgSet* setPdfProjP3; 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; setPdfHadP1 = (RooArgSet*)poissonSigPythP1->getObservables(toyHadP1); setPdfHadP2 = (RooArgSet*)poissonSigPythP2->getObservables(toyHadP2); setPdfHadP3 = (RooArgSet*)poissonSigPythP3->getObservables(toyHadP3); setPdfHadP1->Print(); for(int n = 0; n < toyGen; n++){ //Need to get the get then number of hadrons generated setHadP1 = (RooArgSet*)toyHadP1->get(n); nHad = (RooRealVar*)setHadP1->find(nHad->GetName()); genHadP1 = static_cast(nHad->getVal()); setHadP2 = (RooArgSet*)toyHadP2->get(n); nHad = (RooRealVar*)setHadP2->find(nHad->GetName()); genHadP2 = static_cast(nHad->getVal()); setHadP3 = (RooArgSet*)toyHadP3->get(n); nHad = (RooRealVar*)setHadP3->find(nHad->GetName()); genHadP3 = static_cast(nHad->getVal()); cout << genHadP1 << " " << genHadP2 << " " << genHadP3 << 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); cout << "=========>Hadron xi,projP generated according to nHadrons generated in each plane" << endl; setPdfXiP1 = pdf0xiP1->getObservables(toyXiP1); toyXiP1->Print(); setPdfXiP1->Print(); setPdfXiP2 = (RooArgSet*)pdf0xiP2->getObservables(toyXiP2); setPdfXiP3 = (RooArgSet*)pdf0xiP3->getObservables(toyXiP3); setPdfProjP1 = (RooArgSet*)pdf0projP1->getObservables(toyProjP1); setPdfProjP2 = (RooArgSet*)pdf0projP2->getObservables(toyProjP2); setPdfProjP3 = (RooArgSet*)pdf0projP3->getObservables(toyProjP3); //Loop within toy datasets //Plane1 cout << "=========>Plane 1" << endl; for(int m = 0; m < genHadP1; m++){ setXiP1 = (RooArgSet*)toyXiP1->get(m); xi_s = (RooRealVar*)setXiP1->find(xi_s->GetName()); genXiP1 = xi_s->getVal(); setProjP1 = (RooArgSet*)toyProjP1->get(m); projP = (RooRealVar*)setProjP1->find(projP->GetName()); genProjP1 = projP->getVal(); *setPdfXiP1 = *toyXiP1->get(m); genProbXiP1 = pdf0xiP1->getVal(setPdfXiP1);//Like this or with nothing () - still not returning the access value cout << genXiP1 << " - " << genProjP1 << " | " << genProbXiP1 <Plane 2" << endl; for(int m = 0; m < genHadP2; m++){ setXiP2 = (RooArgSet*)toyXiP2->get(m); xi_s = (RooRealVar*)setXiP2->find(xi_s->GetName()); genXiP2 = xi_s->getVal(); setProjP2 = (RooArgSet*)toyProjP2->get(m); projP = (RooRealVar*)setProjP2->find(projP->GetName()); genProjP2 = projP->getVal(); cout << genXiP2 << " - " << genProjP2 << endl; } //Plane3 cout << "=========>Plane 3" << endl; for(int m = 0; m < genHadP3; m++){ setXiP3 = (RooArgSet*)toyXiP3->get(m); xi_s = (RooRealVar*)setXiP3->find(xi_s->GetName()); genXiP3 = xi_s->getVal(); setProjP3 = (RooArgSet*)toyProjP3->get(m); projP = (RooRealVar*)setProjP3->find(projP->GetName()); genProjP3 = projP->getVal(); cout << genXiP3 << " - " << genProjP3 << endl; } cout << n << endl; } cout << "End." << endl; return 0; }