ActionPlotLimits.cxx

#include "QFramework/TQFolder.h"
#include "QFramework/TQIterator.h"
#include "QFramework/TQUtils.h"
#include "QFramework/TQPathManager.h"

#include "SFramework/TSStatisticsManager.h"
#include "SFramework/TSStatisticsPlotter.h"
#include "SFramework/TSUtils.h"

#include "THistPainter.h"
#include "TFormula.h"

#include "TSystem.h"
#include <limits>

// #define _DEBUG_
#include "QFramework/TQLibrary.h"

#define inf std::numeric_limits<double>::infinity()

/*<cafdoc name=PlotLimits>
  PlotLimits
  ===========================
  
  Plot a limit scan.
  
  Usage:
  ---------------

  The following is an example for a single expected limit.
  ```
  +PlotLimits {
    +allMasses {
      <outputFile="plot.pdf">
      +expected{
        +mh600 {
          # identify the masses and nominal cross sections of each mass point
          <x=600>
          <yScale=1.3207e-02>
          <source="/mh600/Limits/asimov">
        }
        +mh800 {
          # identify the masses and nominal cross sections of each mass point
          <x=800>
          <yScale=9.1080e-03>
          <source="/mh800/Limits/asimov">
        }
        +mh1000 {
          # identify the masses and nominal cross sections of each mass point
          <x=1000>
          <yScale=1.5939e-03>
          <source="/mh1000/Limits/asimov">
        }
        # identify the tags to be used from each data source
        <y="exp_upper_med", yUp.1="exp_upper_p1s", yDn.1="exp_upper_m1s", yUp.2="exp_upper_p2s", yDn.2="exp_upper_m2s"> @ ?;
        # plot cosmetics
        <style.title="expected">
        <style.lineStyle=2>
        <xLabel="m_{H}">
        <yLabel="#sigma_{H}">
      }
    }
  }
  ```

In some cases, you may want to set limits on a model parameter that is not directly proportional to the cross section.
You can use TFormula expressions to translate the y value used for plotting into the model parameter, e.g.
```
        # identify the tags to be used from each data source
        <y="exp_upper_med", yUp.1="exp_upper_p1s", yDn.1="exp_upper_m1s", yUp.2="exp_upper_p2s", yDn.2="exp_upper_m2s", yFormula="sqrt(y)"> @ ?;
```
All tags are expanded automatically.
However, "y" should not be used as a tag in the formula.
If written as a simple letter, it will be treated as a parameter of the TFormula and the uncertainty bands are scaled correctly.
if `yScale` is specified, the result of the formula will be scaled by that factor.
</cafdoc> */

namespace TSBaseActions {
  
  class PlotLimits : public TSStatisticsManager::Action {


    struct HEPdataEntry{
      struct YError {
	double minus;
	double plus;
	std::string label;

	YError(const std::string& lbl, double p, double m): minus(m),plus(p),label(lbl) {};
	YError(const YError& other) = default;
	~YError() = default;
	YError() = default;	
	
      };
      std::string xname;
      std::string yname;
      std::string xunit;
      std::string yunit;            
      std::vector<double> x;
      std::vector<double> y;
      std::vector<std::vector<YError> > yerr;
    };
      
    bool writeHEPdataYML(TQFolder* /*config*/, const HEPdataEntry& data, const TString& outfilename) const {
      
      std::ofstream out(outfilename.Data());
      
      TQUtils::writeHEPdataHeader(out,"independent",data.xname,data.xunit);
      for(size_t i=0; i<data.x.size(); ++i){
	TQUtils::writeHEPdataValue(out,data.x.at(i));
      }
      
      TQUtils::writeHEPdataHeader(out,"dependent",data.yname,data.yunit);
      for(size_t i=0; i<data.y.size(); ++i){
	if(data.yerr.size() > i && data.yerr.at(i).size() > 0){
	  TQUtils::writeHEPdataValueWithErrors(out,data.y.at(i));
	  for(const auto& e:data.yerr.at(i)){
	    TQUtils::writeHEPdataError(out,e.plus,e.minus,e.label);
	  }
	} else {
	  TQUtils::writeHEPdataValue(out,data.y.at(i));
	}
      }
      return true;
    }

    
    bool execute(TQFolder * config) const override {
      
      TQFolder * results = manager->getResults();
      if (!config) {
	return 0;
      }
      
      bool verbose = config->getTagBoolDefault("verbose",false);

      TString outputFile = config->getTagStringDefault("outputFile", "plot.pdf");
      TSStatisticsPlotter pl;
      TString tqpath = TQLibrary::getTQPATH();
      
      TString path = config->getTagStringDefault("arrangement",tqpath+"/../SFramework/share/templates/PlotLimits.txt");
      gSystem->ExpandPathName(path);
      
      TQFolder* plot = TQFolder::newFolder("plot");

      
      int nOK = 0;
      TQFolderIterator itr(config->getListOfFolders("?"));
      while(itr.hasNext()){
        TQFolder* plotconfig = itr.readNext();
        if(!plotconfig) continue;
        
        int nPoints=0;
        bool showBands = plotconfig->getTagBoolDefault("showBands",false);
        bool showN = plotconfig->getTagBoolDefault("showN",false);	
        bool isReference = plotconfig->getTagBoolDefault("isReference",false);
        bool isObserved = plotconfig->getTagBoolDefault("isObserved",!showBands);

	double lastx = -inf;

	HEPdataEntry curve;
	TString xtitle = plotconfig->getTagStringDefault("style.title.xAxis","x");
	TString ytitle = plotconfig->getTagStringDefault("style.title.yAxis","y");		
	curve.xname = TQStringUtils::getWithoutUnit(xtitle);
	curve.xunit = TQStringUtils::getUnit(xtitle);
	curve.yname = TQStringUtils::getWithoutUnit(ytitle) + " (" + plotconfig->getTagStringDefault("style.title") + ")";
	curve.yunit = TQStringUtils::getUnit(ytitle);	
	
        TQFolderIterator points(plotconfig->getListOfFolders("?"));
        while(points.hasNext()){
          TQFolder* point = points.readNext();
          if(!point) continue;

          double x;
	  if(!point->getTagDouble("x",x)) continue;

	  if(x < lastx) continue;
	  lastx = x;

	  double yval = 0.;
	  double yScale = point->getTagDoubleDefault("~yScale",1.);
	  bool useYFormula = false;
	  if(!point->getTagDouble("y",yval)){
	    TString sourcePath = point->getTagStringDefault("~source","Limits/asimov");
	    TQFolder* result = results->getFolder(config->getName()+"+");
	    if(!result){
	      manager->error(TString::Format("unable to access result '%s'",config->GetName()));
	      continue;
	    }
	    TQFolder* source = result->getFolder(sourcePath);
	    if(!source){
	      manager->error(TString::Format("unable to access limit point '%s' from path '%s' from result '%s'",point->getPath().Data(),sourcePath.Data(),config->GetName()));
	      continue;
	    }
	    TString yname = "";
	    if(isObserved) {
	      yname = point->getTagStringDefault("~y","upper");
	    } else {
	      yname = point->getTagStringDefault("~y","exp_upper_med");
	    }

	    TString yFormula = "";
	    if(isReference) { 
	      yval = 1.0;
	    } else {
	      if(!source->getTagDouble(yname,yval)) continue;
	    }
	    
	    TFormula yCalculation;
	    if (point->getTagString("~yFormula", yFormula)) {
	      useYFormula = true;
	      std::string yvalstring = std::to_string(yval);
	      yFormula.ReplaceAll(yname, yvalstring.c_str());
	      yFormula = source->replaceInText(yFormula);
	      yFormula = point->replaceInText(yFormula);
	      DEBUG("creating TFormula yCalculation: %s", yFormula.Data());
	      yCalculation = TFormula("yCalculation", yFormula);
	      // no parameter named x should be in formula
	    }

	    if (useYFormula) {
	      yval = yCalculation.Eval(0, yval);
	    }
	    
	    if(showBands){
	      int i=1;
	      TString ynameUp, ynameDn;
	      std::vector<TString> ynamesUp, ynamesDn;            
	      while(point->getTagString(TString::Format("~yUp.%d",i),ynameUp) && point->getTagString(TString::Format("~yDn.%d",i),ynameDn)){
		++i;
		ynamesUp.push_back(ynameUp);
		ynamesDn.push_back(ynameDn);
	      }
	      if(ynamesUp.size() == 0){
		for(size_t i=0; i<2; ++i){
		  ynamesUp.push_back(TString::Format("exp_upper_p%ds",(int)i+1));
		  ynamesDn.push_back(TString::Format("exp_upper_m%ds",(int)i+1));
		}
	      }
	      std::vector<HEPdataEntry::YError> hepdata_errors;
	      for(size_t i=ynamesUp.size(); i>0; --i){
		double yUp = source->getTagDoubleDefault(ynamesUp[i-1],yval);
		double yDn = source->getTagDoubleDefault(ynamesDn[i-1],yval);
		double yBandUp = 0;
		double yBandDn = 0;
		TQFolder* pt = plot->getFolder(TString::Format("Overlay.%s.band.%ds/p.%d+",plotconfig->GetName(),(int)i,nPoints));
		pt->setTagDouble("x",x);
		pt->setTagDouble("y",yval*yScale);
		if (useYFormula) {
		  yBandUp = yScale*yCalculation.Eval(0, yval + std::abs(yUp-yval));
		  yBandDn = yScale*yCalculation.Eval(0, yval - std::abs(yDn-yval));
		} else {
		  yBandUp = yScale*(yval + std::abs(yUp-yval));
		  yBandDn = yScale*(yval - std::abs(yDn - yval));
		}
		DEBUG("setting uncertainty bands. up: %g, down: %g", yBandUp, yBandDn);
		pt->setTagDouble("yp", yBandUp);
		pt->setTagDouble("yn", yBandDn);
		
		HEPdataEntry::YError hepdata_error(TString::Format("%zusigma",i).Data(),yBandUp - yval * yScale,yBandDn - yval * yScale);		
		hepdata_errors.emplace_back(hepdata_error);
	      }
	      curve.yerr.emplace_back(hepdata_errors);	      
	    }

	    int nToys;
	    if(showN && source->getTagInteger("n",nToys)){
	      TQFolder* text = plot->getFolder(TString::Format("Text.showN.%s.p.%d+",plotconfig->GetName(),nPoints));
	      text->setTagDouble("y",yval*yScale);
	      text->setTagDouble("x",x);
	      text->setTagDouble("style.textSize",0.02);
	      text->setTagString("text",TString::Format("%d",nToys));
	    }
	  }

	  double y = yval * yScale;
	  if(!TQUtils::isNum(x) || !TQUtils::isNum(yval)) continue;
	  DEBUG("adding point x=%g, y=%g", x, y);
	    
          TQFolder* pt = plot->getFolder(TString::Format("Overlay.%s.med/p.%d+",plotconfig->GetName(),nPoints));
          pt->setTagDouble("x",x);
          pt->setTagDouble("y",y);

	  curve.x.push_back(x);
	  curve.y.push_back(y);	  
          
          plotconfig->exportTags(plot->getFolder(TString::Format("Overlay.%s.med",plotconfig->GetName())),"","style.*");
          
	  if(verbose) manager->info(TString::Format("  adding point %g / %g",x,y));
  
          nPoints++;
        }
        if(nPoints>0){
	  manager->info(TString::Format("created plot '%s' with %d points",plotconfig->GetName(),nPoints));
	  nOK++;
	  
	  TString yml;
	  if(plotconfig->getTagString("writeHEPdata",yml)){
	    manager->info(TString::Format("writing '%s'",yml.Data()));
	    writeHEPdataYML(config,curve,yml);
	  }
	}
      }

      plot->setGlobalOverwrite(false);
      plot->importFromTextFile(TQPathManager::getPathManager()->getTargetPath(path).c_str());
      TString addElements;
      if(config->getTagString("addElements",addElements)){
        TString errMsg;
        if(!plot->importFromText(addElements,errMsg)){
          manager->error(errMsg);
        }
      }
      plot->setGlobalOverwrite(true);
      plot->importTags(config);
      TCanvas * c = pl.plot(plot);
      if (c) {
        c->SaveAs(TQPathManager::getPathManager()->getTargetPath(outputFile).c_str());
        delete c;
      } 
      delete plot;
      return c;
    };
        
  };  
    
    
  namespace {
    bool available = TSStatisticsManager::registerAction(new PlotLimits(),"PlotLimits");
  }
}