Add together an arbitrary number of signals plus background

Newbie
1

Hello,

The title explains what I want to do. If it matters, I am fitting Brillouin spectra and don’t always know how many peaks there will be or which I want to fit.

This what I have so far:

void region_1() {

  TH1F* hist = ReadBS("Region_1.csv");
  if ( hist->GetEntries() == 0 ) // Empty histo.
    return;

  // hist->Draw("AP");
  
  // Import into RooDataHist
  double xlow, xhigh;
  xlow = hist->GetXaxis()->GetXmin();
  xhigh = hist->GetXaxis()->GetXmax();
  
  RooRealVar x("x","x",xlow,xhigh);
  RooDataHist data("data","dataset with x",x,hist);
  RooPlot* frame = x.frame() ;
  data.plotOn(frame) ;
  frame->Draw();

  std::vector<lorentzian> lorentzians;
  
  // Lorentzian
  // A/pi ((G/2)/((x-P)^2 + (G/2)^2))
  // RooRealVar A1("A1","Amplitude of Lorentzian", 0.5, 0,100);
  // RooRealVar P1("P1", "Peak of Lorentzian", -74,-77, -72);
  // RooRealVar G1("G1", "Width of lorentzian",0.1,0,2);
  // RooGenericPdf L1{"L1", "A1 * ((G1/2)/((x-P1)^2 + (G1/2)^2))", RooArgList(A1,P1,G1,x)};

  // The info used here should come from TSpectrum
  for (int i = 0; i < 2; ++i) {
    lorentzian l;
    if ( i == 0 ) {
    l.A = 0.5; l.Alow = 0.0; l.Ahigh = 100.0;
    l.P = -74; l.Plow = -77; l.Phigh = -72;
    l.G = 0.1; l.Glow = 0.0; l.Ghigh = 2.0;
    } else {// -77,-80, -75
    l.A = 0.5; l.Alow = 0.0; l.Ahigh = 100.0;
    l.P = -77; l.Plow = -80; l.Phigh = -75;
    l.G = 0.1; l.Glow = 0.0; l.Ghigh = 2.0;
    }
    lorentzians.push_back(l);
  }
  std::cout << "lorentsians: " << lorentzians.size() << "\n";

  // DEBUG
  // for (int i = 0; i < 10; ++i)
  //   std::cout << "counter: " << UNIQUE_ID(poo) << "\n";
  
  // Now construct a generic PDF for our lorentzian
  std::vector<RooGenericPdf> gpdfs;
  std::vector<RooArgList> args;
  for (int i = 0; i < lorentzians.size(); ++i) {
    RooRealVar A("A", "Amplitude of Lorentzian", lorentzians.at(i).A,lorentzians.at(i).Alow,lorentzians.at(i).Ahigh);
    RooRealVar P("P", "Peak of Lorentzian", lorentzians.at(i).P,lorentzians.at(i).Plow,lorentzians.at(i).Phigh);
    RooRealVar G("G", "FWHM of Lorentzian", lorentzians.at(i).G,lorentzians.at(i).Glow,lorentzians.at(i).Ghigh);
    RooArgList arg;
    arg.add(A);
    arg.add(P);
    arg.add(G);
    arg.add(x);
    RooGenericPdf L{"L", "A * ((G/2)/((x-P)^2 + (G/2)^2))", arg};
    // RooGenericPdf LL(L);
    gpdfs.push_back(L);
  }
  std::cout << "gpdfs: " << gpdfs.size() << "\n";
  std::cout << "args: " << args.size() << "\n";
  
  // Need to sum L1 and L2 and fit that.
  RooRealVar f("f","f",0,1);
  RooAddPdf S1("S1", "Sum of Lorentians", args,f);

  S1.fitTo(data);
  // S1.plotOn(frame);
  // frame->Draw();

  // How would one assemble an arbitrary number of peaks plus a background?
  delete hist;
}

It won’t compile because RooAddPdf does not take a std::vector of RooArgList. I have looked through some of the examples in the tutorials but I can’t find one that clearly does what I want.

I appreciate any help, and by all means point me to the tutorial that I may have missed.

sprock

Please read tips for efficient and successful posting and posting code

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided

3

Hi,

Indeed, as you say, the RooAddPdf can be constructed with a RooArgList and not a vector ROOT: RooAddPdf Class Reference .

Maybe I am missing something, apologies if that’s the case.

Cheers,
Danilo

4

Hello @sprock,

I think there is a misunderstanding in how PDFs and their coefficients are composed. You don’t need to mention their parameters again when making the RooAddPdf. See for example this tutorial:

They add two PDFs, and they do it using:

   RooGaussian sig1("sig1", "Signal component 1", x, mean, sigma1);
   RooGaussian sig2("sig2", "Signal component 2", x, mean, sigma2);
 
   // Sum the signal components into a composite signal pdf
   RooRealVar sig1frac("sig1frac", "fraction of component 1 in signal", 0.8, 0., 1.);
   RooAddPdf sig("sig", "Signal", RooArgList(sig1, sig2), sig1frac);

You see that by putting N pdfs in the RooArgList constructor, you could have summed over N PDFs. The fourth argument can be used to determine the relative fractions of the PDFs. Since the PDFs need to sum to 1, you use N-1 fraction coefficients for summing over N PDFs.