Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

PDFs for B2G-23-008 analysis #876

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Conversation

anmalara
Copy link

@anmalara anmalara commented Dec 4, 2023

No description provided.

@guitargeek
Copy link
Contributor

guitargeek commented Dec 4, 2023

Hello @anmalara! I have tome questions and suggestions about this:

why not use RooGenericPdf for simple expression pdfs like this? Adding simple custom pdfs that only implement evaluate() anyway has several disadvantages:

  • It's a big commitment for backwards compatibility! With the ClassDef(ExpGaussExp,1), you make your class serializable, meaning people can workspaces with it in a file. That puts a burden on the combine maintainers, because in principle these classes can never be changed or removed without causing possible problems because people can't open their workspace
  • By only implementing the evaluate function, and not other RooFit interfaces like for the vectorized/GPU evaluation or C++ code generation from the model, you are not benefiting from the new stuff that RooFit has to offer
  • It's just a lot of duplicated code

Maybe an alternative approach to this could be to simply define a header with free inlined functions to use in RooGenericPdfs? That would be a slimmer and more future-proof solution, and it could even be inspected by compiler plugins for automatic differentiation like Clad because it's in the headers. E.g.:

#ifndef customPdfs_h
#define customPdfs_h

inline double expGaussExp(double x, double p0, double p1, double p2, double p3)
{
  Double_t std=(x-p0)/p1;
  Double_t result=0;

  if (std<-p2){
    result=exp(p2*p2/2.+p2*std);
  } else if (-p2<=std && std<p3){
    result=exp(-0.5*pow(std, 2));
  } else {
    result=exp(p3*p3/2.-p3*std);
  }

  return result;
}

#endif // customPdfs_h

Then it could be used like this C++:

#include <TInterpreter.h>
#include <RooGenericPdf.h>
#include <RooRealVar.h>
#include <RooWorkspace.h>

void initCustomPdfs()
{
    static bool hasInit = false;
    if(hasInit) return;
    gInterpreter->ProcessLine("#include \"customPdfs.h\"");
    hasInit = true;
}

auto expGaussExpFormula = "expGaussExp(x[0], x[1], x[2], x[3], x[4])";

void example()
{
    initCustomPdfs();

    RooRealVar x{"x", "x", 1.0, 0.0, 10.0};
    RooRealVar p0{"p0", "p0", 1.0, 0.0, 10.0};
    RooRealVar p1{"p1", "p1", 1.0, 0.0, 10.0};
    RooRealVar p2{"p2", "p2", 1.0, 0.0, 10.0};
    RooRealVar p3{"p3", "p3", 1.0, 0.0, 10.0};

    RooGenericPdf expGaussExp{"expGaussExp", expGaussExpFormula, {x, p0, p1, p2, p3}};

    expGaussExp.Print();
}

That would at least be my suggestion from my role as the RooFit maintainer :) In ROOT, we put quite a bit an emphasis on the interfaces where users pass string formulas, because that just very flexible (see also RDataFrame).

Sorry for chiming in as a non-CMS member, but I think it's important to know there are superior and more maintainable alternatives. I'm also very much welcoming other ideas and suggestions on how to make the RooGenericPdf/RooFormulaVar and friends work better!

@anmalara
Copy link
Author

anmalara commented Dec 5, 2023

Hi @guitargeek,

thanks for your feedback. I agree with the idea of having simplified code. I wasn't aware of the method you suggested. I think I have a compiled version with the changes you suggested for all functions. I cannot test if it actually works because the workspaces I had used the old class definition. This will take longer to check since I would need to recreate them first.

@anmalara
Copy link
Author

hi @guitargeek,

sorry to go back to this only now.
I stand corrected and even though the compilation works, I still didn't manage to get my code running.
Actually, trying to run your code I get this error. Do you have any suggestion on how to fix this?
If it helps, I'm running with ROOT 6.14.

Thanks,
Andrea

Error in <RooFormula::Compile>:  Bad numerical expression : "expGaussExp(x[0],x[1],x[2],x[3],x[4])"
RooGenericPdf::expGaussExp[ actualVars=(x,p0,p1,p2,p3) formula="expGaussExp(x[0], x[1], x[2], x[3], x[4])" ] = [#0] ERROR:InputArguments -- RooFormula::RooFormula(expGaussExp): compile error
[#0] ERROR:Eval -- RooFormula::eval(expGaussExp): Formula doesn't compile: expGaussExp(x[0], x[1], x[2], x[3], x[4])
[#0] ERROR:Eval -- RooFormula::eval(expGaussExp): Formula doesn't compile: expGaussExp(x[0], x[1], x[2], x[3], x[4])
0

@guitargeek
Copy link
Contributor

guitargeek commented Dec 12, 2023

Hi! I can really recommend updating ROOT here. The 6.14 series is from 2018, which is quite old, and the problem is fixed since ROOT 6.20.

In the combine docs, the recommended environment is CMSSW_11_3_4, which comes with ROOT 6.22.09. What is the specific reason for using the much older environment?

By the way, also according to the combine docs, the previously recommended CMSSW version was CMSSW_10_2_X with combine v8, which used ROOT 6.12.07. So also given the fact that you are using a ROOT version that combine doesn't support explicitly, it would probably be very good to upgrade ROOT.

@anmalara
Copy link
Author

Hi @guitargeek,

Thanks for looking into it. I see the problem.
I am using the versions you suggested for combine, but the creating of the workspace is done at analysis level, which uses CMSSW 106, as per CMS recommendations.

Although I see the point of moving to the latest recommended version, this would require some extra work needed in my workflow to update to these different releases. I'd prefer to not move at this stage of my analysis (unblinded already).
I'm not sure how to proceed from here. If it's fine with you, I'd suggest to keep this PR as it is since it follows the same structure of the other class definition. Even though it's not optimal, maybe a different PR can take care of your suggestions?

Maybe @anigamova have some suggestions since I'd like to have a working setup in the next release.
Thanks :-)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants