I have been trying desperately to create a likelihood ratio test statistic but keep coming up against problems evaluating pdfs. I have one code which is working using the createNLL function but I want to also evaluate it by hand with the pdfs I am using to verify that createNLL is working as I think it is, as previously I had different answers.
I have decided to rewrite my code from scratch making better use of functions in order to control the flow of the program but I seem to have a problem. I have generated 1D RooKeysPdfs using RooDataSets created from TTree outputs I have generated using my analysis program over some MC backgrounds. They appear to be working okay, ie: I can plot them and they look as I would expect them to.
The issue I have now is that I am using the advice Wouter gave me [url=https://root-forum.cern.ch/t/accessing-pdf-values-for-different-x-values-rookeyspdf/14194/1 ignoring the other stuff I have written there. The issue now seems to be that the values I am getting returned to me are, more often than not, greater than 1! Seeing as these should be probabilities, I don’t understand where this has come from.
I cannot post a working model of this without uploading my workspaces, but I will attach my code. The problem lies in the EvaluateEvent() method, but as I am not getting any sort of errors, I don’t know what to do. I think I am doing it correctly, but I do not have anyone I can ask to look over my code who knows more about RooFit than I do unfortunately.
As far as I can tell, I have either implemented the function wrong, or my pdf is incorrect. I notice that when I naively do:
pdf->createIntegral(*variable)
the value I get out is less than 1; 0.968955 to be precise, but as I might not be applying the range correctly I do not know if this is sensible to take into consideration. If this should return 1 though, it implies that RooKeysPdf did not generate correctly. If it is of any consequence, the x-range of this pdf is 0.1->0.9.
If someone could have a look at my code, I would very much appreciate it. The EvaluateEvent() method is reproduced here:
[code]void EvaluateEvent(int poisP1, int poisP2, int poisP3, RooDataSet* xiP1, RooDataSet* xiP2, RooDataSet* xiP3){
//Easy part - Poisson probabilities
double probnHadGenP1 = TMath::Poisson(poisP1, meanGenP1);
double probnHadGenP2 = TMath::Poisson(poisP2, meanGenP2);
double probnHadGenP3 = TMath::Poisson(poisP3, meanGenP3);
double probnHadComP1 = TMath::Poisson(poisP1, meanComP1);
double probnHadComP2 = TMath::Poisson(poisP2, meanComP2);
double probnHadComP3 = TMath::Poisson(poisP3, meanComP3);
//Harder part - Product of probs
double probXiGenP1;
//Plane 1
RooArgSet* setPdfXiGenP1 = pdfXiGenP1->getObservables(xiP1);
RooArgSet* setPdfXiComP1 = pdfXiComP1->getObservables(xiP1);
RooArgSet* setXiP1;
cout << "Int" << pdfXiGenP1->createIntegral(*xi)->getVal() << endl;
for (int j = 0; j < poisP1; j++){
//setXiP1 = (RooArgSet*)xiP1->get(j);
//xi = (RooRealVar*)setXiP1->find(xi->GetName());
*setPdfXiGenP1 = *xiP1->get(j);
probXiGenP1 = pdfXiGenP1->getVal(setPdfXiGenP1);
cout << xi->getVal() << " " << probXiGenP1 << endl;
}
[/code]
A typical output of value and probability from pdf is :
0.838481 1.44625
0.158841 2.10157
0.629804 0.944437
0.280225 1.39589
0.239268 1.52438
0.317163 1.28884
0.155004 2.13711
0.245083 1.49627
I really cannot work out what could be causing this.
Thanks,
Ian
EDIT: I have noticed that ->createIntegral(*xi,*xi)->getVal() does indeed return 1 which is what should happen when I specify the integration as noted [url=https://root-forum.cern.ch/t/getting-integral-of-fitted-pdf/13779/4 but this doesn’t solve my issue.
code copy.txt (7.86 KB)