Intro:

This section will serve as a documentation for SideBandSubtractionUnitTest.C as well as another example of how SideBandSubtraction.C can be used. More information on SideBandSubtraction.C can be found here. What the unit test does is, makes a 2D "signal" PDF, a 2D "background" PDF, and then adds them together and feeds the sum PDF into SideBandSubtraction.C. Since the ultimate goal of SidebandSubtraction.C is to get a signal PDF from a Signal + Background PDF, the output of SideBandSubtraction.C should match the signal PDF that we had at the beginning. The purpose of this test is so that a user can run this test to get a feel for how SideBandSubtraction.C works and test to make sure that their own PDF's are working properly. Note that at this time this page is in a state of flux, so things may change abruptly.

Customization:

At this time the signal PDF consists of a Gaussian along a "mass" axis and another gaussian along a "pt" axis and a background PDF consisting of a flat line PDF along the mass axis and another Gaussian along the pt axis. There were plans to allow users to specify their own functions and then run the unit_test() function on them, but since this would require so much input from the user, it was impractical. Instead, the user should just edit the PDF's that get put into the RooProdPdf in the SideBandSubractionUnitTest.C function. The user can also switch the bool verbose from false to true to have graphs output. Further details on this are found in the documentation as well as in comments in the code when applicable. In the code, the areas a user can edit are marked "TO USERS". These include more detailed explanations than can be provided here. Be sure to check all areas marked "TO USERS" to ensure that there are no problems.

Obtaining and Running:

At this point, the unit test is included with SideBandSubtract. Following the steps in that page to download SideBandSubtract will also download the unit test. To run the unit test, simply run the following commands:

  cd UserCode/DavidBjergaard/test/
  root -l -b -q run_unit_test.C

To edit the PDF's yourself, simply navigate the UserCode/DavidBjergaard/src/ and edit SideBandSubtractionUnitTest.C.

After the ROOT code runs the script will then output a status message. This tells why (or if) the test failed. At this point the possible messages are:

1. No Problems
2. Kolmogorov test gives unsatisfactory results
3. bg_pdf < 2-dimensions
4. bg_pdf: #events before cut < # events after cut
5. signal_pdf < 2-dimensions
6. bg_pdf: #events before cut < # events after cut
7. Defined signal/sideband regions overlap
8. Fails Chi-Square test
9. Pulls [(fitted_value-original_value)/fitting_error] too large
10. Number of interest variables does not equal (total number of variables)-1

After this message the run_unit_test.C script will do a little clean up. This involves making a folder for any output called plots. This is done automatically using mkdir, so if a folder named plots already exists it just doesn't do anything. After that it moves any .gifs and .eps files to this plots folder. If the verbose option is turned of (see the "customization" section for more info on this) then there are no files output for it to move. These shouldn't cause any problems, but users should be aware that these are in place in case they have .gif or .eps files in this folder. Of course, a user could just edit run_unit_test.C to their needs.

NOTES:

If your rootlogon.C file has the following lines:

   gSystem->Load("libRooFit");
   using namespace RooFit;

The program won't run properly the first time you run it. Simply running it again should fix the problem.

Overview:

The unit test is essentially a 4 step process.

1. Make background and signal PDF's
2. Add them together
3. Make an SBS object a process it, using the plots that were generated in (1) and (2).
4. Make sure that the output of the SideBandSubtraction.C (3) matches the input signal PDF(1)

Step by Step:

Step 1:

First we define our variables:

  const float m1 = 3.2;//mean for bg_pdf over pt
  const float m2 = 3.1;//mean for signal_pdf over mass
  const float m3 = 2.5;//mean for signal_pdf over pt
  const float s1 = .2;//sigmas, numbered the same as the means
  const float s2 = .1;
  const float s3 = .2;
  RooRealVar mass("mass","mass",2,4,"GeV");//separation variable
  RooRealVar pt("pt","pt",2,4,"GeV");
  RooRealVar mean1("mean1","mean_of_gaussian_1",m1,2,4,"GeV");
  RooRealVar sigma1("sigma1","width_of_gaussian_1",s1,"GeV");
  RooRealVar mean2("mean2","mean_of_gaussian_2",m2,2,4,"GeV");
  RooRealVar sigma2("sigma2","width_of_gaussian_2",s2,"GeV");
  RooRealVar mean3("mean3","mean_of_gaussian_3",m3,2,4,"GeV");
  RooRealVar sigma3("sigma3","width_of_gaussian_3",s3,"GeV");

Then make 4 pdf's (background over mass, background over pt, signal over mass, signal over pt):

  RooPolynomial bg("bg","Constant Background",mass,RooArgList());
  RooGaussian gauss1("gauss1","Gaussian of PDF 1",pt,mean1,sigma1);
  RooGaussian gauss2("gauss2","Gaussian_of_PDF_2",mass,mean2,sigma2);
  RooGaussian gauss3("gauss3","Gaussian_of_PDF_3",pt,mean3,sigma3);

Visually, these are:

Then, to finish make 2 RooProdPdf's, bg and signal:

  RooProdPdf bg_pdf("signal pdf","bg*mass",RooArgSet(bg,gauss1));
  RooProdPdf signal_pdf("signal_pdf","gauss2*gauss3",RooArgSet(gauss2,gauss3));

NOTES:

These PDF's can be anything, as long as they act like real data. If they do not, it may cause the SideBandSubtract object to segfault. If you decide the change the PDF's to be put through the unit test, just make sure the RooProdPdf's have the names specified above or the rest of code will not work. Any other PDF names/types can be the same.

Step 2:

Add the two to make our final 2D pdf.

  RooRealVar gfrac("gfrac","fraction of gaussian",0.5);
  RooAddPdf sum_pdf("sum_pdf","bg_pdf+signal_pdf",RooArgList(bg_pdf,signal_pdf),gfrac);

Visually, these look like:
Over mass:

Over pt:

Step 3:

Make the SBS Object and run the doGlobalFit() function:

  TH1F massHist("mass","mass",100,2,4);
  TH1F ptHist("pt","pt",100,2,4);

  vector<TH1F*> basehistos(0);
  basehistos.push_back(&massHist);
  basehistos.push_back(&ptHist);

  SideBandSubtract sbs(&sum_pdf,&bg_pdf,datasum,&mass,basehistos,true);

  sbs.addSideBandRegion(2,2.8);
  sbs.addSignalRegion(2.8,3.4);
  sbs.addSideBandRegion(3.4,4.0);
  sbs.doGlobalFit();
  sbs.printResults();
  vector<TH1F> subtractedHistos(sbs.getSBSHistos());

NOTES:

For more information about the arguments for the SideBandSubtract object, check it's main page here. The output of sbs.printResults() can also be found there under the section titled, "Pedagogical step by step explanation of the side-band subtraction method."
Of particular note is that the name of the histos in the basehistos vector need to match the names of the variables used. For example, notice the bold areas is the following: RooRealVar mass(" mass ","mass",2,4,"GeV"); and TH1F massHist(" mass ","mass",100,2,4);. These must match.

Step 4:

We check that the output matches the input:
To check visually, we can just overlay the pt projection of the output histo from doGlobalFit() with the signal pdf over pt (edit SideBandSubtractionUnitTest.C to have bool verbose=true to enable these graphs):

  RooDataHist SBSroodata("SBSroodata","SBSroodata",RooArgList(pt),&subtractedHistos[0]);
  signal_pdf.fitTo(SBSroodata);
  RooPlot* ptframe = pt.frame();
  SBSroodata.plotOn(ptframe);
  signal_pdf.plotOn(ptframe);
  TCanvas canvas;
  ptframe->Draw();
  canvas.SaveAs("overlay.gif");

And to get a numerical value of accuracy, we use a KolmogorovTest, a Chi-square test, and computing the pulls on a fitted histogram:

  RooDataSet *datasig = signal_pdf.generate(RooArgSet(mass,pt),1000);
  TH2F* modelsig = datasig->createHistogram(mass,pt,"1000","1000");
  Float_t KS_result = subtractedHistos[0].KolmogorovTest(modelsig->ProjectionX("gauss3",0,-1));
  cout<<"Test Result: "<< KS_result << endl<<endl;

  RooPlot* ptframe = pt.frame();
  SBSroodata.plotOn(ptframe);
  signal_pdf.plotOn(ptframe);
  TCanvas canvas;
  Double_t chiresult;
  chiresult = ptframe->chiSquare();

  float pull = ((fitted-m2)/error);
  if(pull > (3*error)) return 1;

NOTES:

For the sake of the test, and since there needed to be a passing condition, the constraint that the result for the kolmogorov test be less that .1 was imposed. Obviously, for the above graph you can see that the results are quite good and thus the result is <<.1.

-- AnthonyDiChiara - 26 Jun 2009