MorphFunctions.cc

Go to the documentation of this file.

#include "CombineHarvester/CombinePdfs/interface/MorphFunctions.h"
#include <iostream>
#include <set>
#include <vector>
#include <string>
#include "boost/lexical_cast.hpp"
#include "boost/format.hpp"
#include "boost/multi_array.hpp"
#include "TGraphErrors.h"
#include "RooFitResult.h"
#include "RooRealVar.h"
#include "RooDataHist.h"
#include "RooProduct.h"
#include "RooConstVar.h"
#include "CombineHarvester/CombineTools/interface/Logging.h"
 
namespace ch {
 
std::string BuildRooMorphing(RooWorkspace& ws, CombineHarvester& cb,
                      std::string const& bin, std::string const& process,
                      RooAbsReal& mass_var, std::string norm_postfix,
                      bool allow_morph, bool verbose, bool force_template_limit, TFile * file) {
  // To keep the code concise we'll make some using-declarations here
  using std::set;
  using std::vector;
  using std::string;
  using boost::lexical_cast;
  using boost::multi_array;
  using boost::extents;
 
  // RooFit can be rather noisy on stdout, unless verbose = true we'll kill
  // everything that's not an error, and restore the user's MsgLevel at the end
  RooFit::MsgLevel backup_msg_level =
      RooMsgService::instance().globalKillBelow();
  if (!verbose) RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
 
  // Everything here works on the assumption that the bin,process,mass
  // combinations we find in this CH instance refer to single processes. So bin
  // should uniquely label an event category, and mass should be defined and
  // convertible from string to float for every process entry.
  if (verbose) {
    std::cout << ">> Bin: " << bin << "  Process: " << process << "\n";
  }
  TString key = bin + "_" + process;
 
  // For the sake of efficiency, take a filtered shallow-copy of the input CH
  // instance that just contains the bin,process combination we've been asked to
  // turn into a RooMorphingPdf
  CombineHarvester cb_bp = cb.cp().bin({bin}).process({process});
 
  // Get a vector of the mass values
  vector<string> m_str_vec = Set2Vec(cb_bp.SetFromProcs(
      std::mem_fn(&ch::Process::mass)));
 
  // We now sort these based on numerical value, hence the lexical_cast to double
  // At the moment this also serves as our check that all mass values are
  // float-convertible, as an exception will be thrown by lexical_cast if not.
  std::sort(m_str_vec.begin(), m_str_vec.end(),
    [](string const& s1, string const& s2) {
      return lexical_cast<double>(s1) < lexical_cast<double>(s2);
    });
 
  // Convert the sorted vector of mass strings to an actual vector of doubles
  vector<double> m_vec;
  for (auto const& s : m_str_vec) {
    if (verbose) std::cout << ">>>> Mass point: " << s << "\n";
    m_vec.push_back(lexical_cast<double>(s));
  }
  // So, we have m mass points to consider
  unsigned m = m_vec.size();
 
  // ss = "shape systematic"
  // Make a list of the names of shape systematics affecting this process
  vector<string> ss_vec =
      Set2Vec(cb_bp.cp().syst_type({"shape"}).syst_name_set());
  // Now check if all shape systematics are present for all mass points
  for (auto const& s : m_str_vec) {
    if (cb_bp.cp().syst_type({"shape"}).mass({s}).syst_name_set().size() !=
        ss_vec.size()) {
      throw std::runtime_error(FNERROR(
          "Some mass points do not have the full set of shape systematics, "
          "this is currently unsupported"));
    }
  }
 
  unsigned ss = ss_vec.size();  // number of shape systematics
 
  // ls = "lnN systematic"
  // Make a list of the regular lnN normalisation systematics affecting this
  // process
  vector<string> ls_vec =
      Set2Vec(cb_bp.cp().syst_type({"lnN"}).syst_name_set());
  unsigned ls = ls_vec.size();  // number of lnN systematics
 
  // Create a bunch of empty arrays to store the information we need in a more
  // convenient format. We use the boost multi_array class because it's a nice
  // multi-dimensional array implementation, and safer to use than standard
  // C-arrays
 
  // Store pointers to each ch::Process (one per mass point) in the CH instance
  multi_array<ch::Process *,  1> pr_arr(extents[m]);
  // Store pointers pointers to each ch::Systematic for each mass point (hence
  // an ss * m 2D array)
  multi_array<ch::Systematic *, 2> ss_arr(extents[ss][m]);
  // With shape systematics we have to support cases where the value in the
  // datacard is != 1.0, i.e. we are scaling the parameter that goes into the
  // Gaussian constraint PDF - we have to pass this factor on when we build the
  // normal vertical-interpolation PDF for each mass point. We will store these
  // scale factors in this array, one per shape systematic.
  multi_array<double, 1> ss_scale_arr(extents[ss]);
  // Really just for book-keeping, we'll set this flag to true when the shape
  // systematic scale factor != 1
  multi_array<bool, 1> ss_must_scale_arr(extents[ss]);
  // Similar to the ss_arr above, store the lnN ch::Systematic objects. Note
  // implicit in all this is the assumption that the processes at each mass
  // point have exactly the same list of systematic uncertainties.
  multi_array<ch::Systematic *, 2> ls_arr(extents[ls][m]);
 
  // This array holds pointers to the RooRealVar objects that will become our
  // shape nuisance parameters, e.g. "CMS_scale_t_mutau_8TeV"
  multi_array<std::shared_ptr<RooRealVar>, 1> ss_scale_var_arr(extents[ss]);
  // And this array holds the constant scale factors that we'll build from
  // ss_scale_arr above
  multi_array<std::shared_ptr<RooConstVar>, 1> ss_scale_fac_arr(extents[ss]);
  // Finally this array will contain the scale_var * scale_fac product where we
  // need to provide a scaled value of the nuisance parameter instead of the
  // parameter itself in building the vertical-interp. PDF
  multi_array<std::shared_ptr<RooProduct>, 1> ss_scale_prod_arr(extents[ss]);
 
 
  // Now let's fill some of these arrays...
  for (unsigned mi = 0; mi < m; ++mi) {
    // The ch::Process pointers
    cb_bp.cp().mass({m_str_vec[mi]}).ForEachProc([&](ch::Process *p) {
      pr_arr[mi] = p;
    });
    for (unsigned ssi = 0; ssi < ss; ++ssi) {
      // The ch::Systematic pointers for shape systematics
      cb_bp.cp().mass({m_str_vec[mi]}).syst_name({ss_vec[ssi]})
        .ForEachSyst([&](ch::Systematic *n) {
            ss_arr[ssi][mi] = n;
        });
    }
    for (unsigned lsi = 0; lsi < ls; ++lsi) {
      // The ch::Systematic pointers for lnN systematics
      cb_bp.cp().mass({m_str_vec[mi]}).syst_name({ls_vec[lsi]})
        .ForEachSyst([&](ch::Systematic *n) {
            ls_arr[lsi][mi] = n;
        });
    }
  }
 
  // We need to build a RooArgList of the vertical morphing parameters for the
  // vertical-interpolation pdf - this will be the same for each mass point so
  // we only build it once
  RooArgList ss_list;
  for (unsigned ssi = 0; ssi < ss; ++ssi) {
    // Make the nuisance parameter var. We use shared_ptrs here as a convenience
    // because they will take care of automatically cleaning up the memory at
    // the end
    ss_scale_var_arr[ssi] =
        std::make_shared<RooRealVar>(ss_vec[ssi].c_str(), "", 0);
    // We'll make a quick check that the scale factor for this systematic is the
    // same for all mass points. We could do a separate scaling at each mass
    // point but this would create a lot of complications
    set<double> scales;
    // Insert the scale from each mass point into the set, if it has a size
    // larger than one at the end then we have a problem!
    for (unsigned mi = 0; mi < m; ++mi) {
      scales.insert(ss_arr[ssi][mi]->scale());
    }
    if (scales.size() > 1) {
      // Don't let the user proceed, we can't build the model they want
      std::runtime_error(FNERROR(
          "Shape morphing parameters that vary with mass are not allowed"));
    } else {
      // Everything ok, set the scale value in its array
      ss_scale_arr[ssi] = *(scales.begin());
      // Handle the case where the scale factor is != 1
      if (std::fabs(ss_scale_arr[ssi] - 1.0) > 1E-6) {
        ss_must_scale_arr[ssi] = true;
        // Build the RooConstVar with the value of the scale factor
        ss_scale_fac_arr[ssi] = std::make_shared<RooConstVar>(
            TString::Format("%g", ss_scale_arr[ssi]), "",
            ss_scale_arr[ssi]);
        // Create the product of the floating nuisance parameter and the
        // constant scale factor
        ss_scale_prod_arr[ssi] = std::make_shared<RooProduct>(
            ss_vec[ssi] + "_scaled_" + key, "",
            RooArgList(*(ss_scale_var_arr[ssi]), *(ss_scale_fac_arr[ssi])));
        // Add this to the list
        ss_list.add(*(ss_scale_prod_arr[ssi]));
      } else {
        // If the scale factor is 1.0 then we just add the nuisance parameter
        // directly to our list
        ss_list.add(*(ss_scale_var_arr[ssi]));
      }
    }
  }
 
  // Summarise the info on the shape systematics and scale factors
  if (verbose) {
    std::cout << ">> Shape systematics: " << ss << "\n";
    for (unsigned ssi = 0; ssi < ss; ++ssi) {
      std::cout << boost::format("%-50s %-5i %-8.3g\n")
        % ss_vec[ssi] % ss_must_scale_arr[ssi] % ss_scale_arr[ssi];
    }
  }
 
 
  // lms = "lnN morphing systematic"
  // Now we have some work to do with the lnN systematics. We can consider two cases:
  //  a) The uncertainty is the same for each mass point => we can leave it in
  //     the datacard as is and let text2workspace do its normal thing
  //  b) The uncertainty varies between mass points => we can't capture this
  //     information in the text datacard in the usual way, so we'll build a RooFit
  //     object that effectively makes the lnN uncertainty a function of the mass
  //     variable
  // We'll use "lms" to refer to case b), which we'll try to figure out now...
  vector<string> lms_vec;  // vec of systematic names
  set<string> lms_set;     // set of systematic names
  // index positions in our full ls_arr array for the lms systematics
  vector<unsigned > lms_vec_idx;  
  for (unsigned lsi = 0; lsi < ls; ++lsi) {
    // Extra complication is that the user might have been evil and mixed
    // symmetric and asymmetric lnN values, we'll try and detect changes in
    // either
    set<double> k_hi;
    set<double> k_lo;
    // Go through each mass point for this systematic and add the uncertainty
    // values (so-called "kappa" values)
    for (unsigned mi = 0; mi < m; ++mi) {
      Systematic *n = ls_arr[lsi][mi];
      k_hi.insert(n->value_u());
      if (n->asymm()) {
        k_lo.insert(n->value_d());
      }
    }
    // If either of these sets has more than one entry then this is a lms case
    if (k_hi.size() > 1 || k_lo.size() > 1) {
      lms_vec.push_back(ls_vec[lsi]);
      lms_set.insert(ls_vec[lsi]);
      lms_vec_idx.push_back(lsi);
    }
  }
  unsigned lms = lms_vec.size();
  // New array for the pointers to the lms Systematic objects
  multi_array<ch::Systematic *, 2> lms_arr(extents[lms][m]);
  for (unsigned lmsi = 0; lmsi < lms; ++lmsi) {
    for (unsigned mi = 0; mi < m; ++mi) {
      lms_arr[lmsi][mi] = ls_arr[lms_vec_idx[lmsi]][mi];
    }
  }
  
  // We will need to create the nuisance parameters for these now
  multi_array<std::shared_ptr<RooRealVar>, 1> lms_var_arr(extents[lms]);
  for (unsigned lmsi = 0; lmsi < lms; ++lmsi) {
    lms_var_arr[lmsi] =
        std::make_shared<RooRealVar>(lms_vec[lmsi].c_str(), "", 0);
  }
 
  // Give a summary of the lms systematics to the user
  if (verbose) {
    std::cout << ">> lnN morphing systematics: " << lms << "\n";
    for (unsigned lmsi = 0; lmsi < lms; ++lmsi) {
      std::cout << ">>>> " << lms_vec[lmsi] << "\n";
    }
  }
 
  // Now we move into a phase of building all the objects we need:
 
  // 2D array of all input histograms, size is (mass points * (nominal +
  // 2*shape-systs)). The factor of 2 needed for the Up and Down shapes
  multi_array<std::shared_ptr<TH1F>, 2> hist_arr(extents[m][1+ss*2]);
  // We also need the array of process yields vs mass, because this will have to
  // be interpolated too
  multi_array<double, 1> rate_arr(extents[m]);
  // Also store the yield uncertainty - we don't actually need this for the signal
  // model, but will include it in the debug TGraph
  multi_array<double, 1> rate_unc_arr(extents[m]);
  // The vertical-interpolation PDF needs the TH1 inputs in the format of a TList
  multi_array<std::shared_ptr<TList>, 1> list_arr(extents[m]);
  // Combine always treats the normalisation part of shape systematics as
  // distinct from the actual shape morphing. Essentially the norm part is
  // treated as an asymmetric lnN. We have to make the kappa_hi and kappa_lo a
  // function of the mass parameter too, so we need two more arrays in (ss * m)
  multi_array<double, 2> ss_k_hi_arr(extents[ss][m]);
  multi_array<double, 2> ss_k_lo_arr(extents[ss][m]);
  // For each shape systematic we will build a RooSpline1D, configured to
  // interpolate linearly between the kappa values
  multi_array<std::shared_ptr<RooSpline1D>, 1> ss_spl_hi_arr(extents[ss]);
  multi_array<std::shared_ptr<RooSpline1D>, 1> ss_spl_lo_arr(extents[ss]);
  // To define the actually process scaling as a function of these kappa values,
  // we need an AsymPow object per mass point
  multi_array<std::shared_ptr<AsymPow>, 1> ss_asy_arr(extents[ss]);
 
  // Similar set of objects needed for the lms normalisation systematics
  multi_array<double, 2> lms_k_hi_arr(extents[lms][m]);
  multi_array<double, 2> lms_k_lo_arr(extents[lms][m]);
  multi_array<std::shared_ptr<RooSpline1D>, 1> lms_spl_hi_arr(extents[lms]);
  multi_array<std::shared_ptr<RooSpline1D>, 1> lms_spl_lo_arr(extents[lms]);
  multi_array<std::shared_ptr<AsymPow>, 1> lms_asy_arr(extents[lms]);
 
  // Now we'll fill these objects..
 
  for (unsigned mi = 0; mi < m; ++mi) {
    // Grab the nominal process histograms. We also have to convert every
    // histogram to a uniform integer binning, because this is what
    // text2workspace will do for all the non-morphed processes in our datacard,
    // and we need the binning of these to be in sync.
    hist_arr[mi][0] =
        std::make_shared<TH1F>(RebinHist(AsTH1F(pr_arr[mi]->shape())));
    // If the user supplied a TFile pointer we'll dump a bunch of info into it
    // for debugging
    if (file) {
      file->WriteTObject(pr_arr[mi]->shape(), key + "_" + m_str_vec[mi]);
    }
    // Store the process rate
    rate_arr[mi] = pr_arr[mi]->rate();
    auto proc_hist = pr_arr[mi]->ClonedScaledShape();
    proc_hist->IntegralAndError(1, proc_hist->GetNbinsX(), rate_unc_arr[mi]);
    // Do the same for the Up and Down shapes
    for (unsigned ssi = 0; ssi < ss; ++ssi) {
      hist_arr[mi][1 + 2 * ssi] =
          std::make_shared<TH1F>(RebinHist(AsTH1F(ss_arr[ssi][mi]->shape_u())));
      hist_arr[mi][2 + 2 * ssi] =
          std::make_shared<TH1F>(RebinHist(AsTH1F(ss_arr[ssi][mi]->shape_d())));
      if (file) {
        file->WriteTObject(ss_arr[ssi][mi]->shape_u(),
                           key + "_" + m_str_vec[mi] + "_" + ss_vec[ssi] + "Up");
        file->WriteTObject(ss_arr[ssi][mi]->shape_d(),
                           key + "_" + m_str_vec[mi] + "_" + ss_vec[ssi] + "Down");
      }
      // Store the uncertainty ("kappa") values for the shape systematics
      ss_k_hi_arr[ssi][mi] = ss_arr[ssi][mi]->value_u();
      ss_k_lo_arr[ssi][mi] = ss_arr[ssi][mi]->value_d();
      // For the normalisation we scale the kappa instead of putting the scaling
      // parameter as the variable
      if (std::fabs(ss_scale_arr[ssi] - 1.0) > 1E-6) {
        ss_k_hi_arr[ssi][mi] = std::pow(ss_arr[ssi][mi]->value_u(), ss_scale_arr[ssi]);
        ss_k_lo_arr[ssi][mi] = std::pow(ss_arr[ssi][mi]->value_d(), ss_scale_arr[ssi]);
      }
    }
    // And now the uncertainty values for the lnN systematics that vary with mass
    // We'll force these to be asymmetric even if they're not
    for (unsigned lmsi = 0; lmsi < lms; ++lmsi) {
      lms_k_hi_arr[lmsi][mi] = lms_arr[lmsi][mi]->value_u();
      if (lms_arr[lmsi][mi]->asymm()) {
        lms_k_lo_arr[lmsi][mi] = lms_arr[lmsi][mi]->value_d();
      } else {
        lms_k_lo_arr[lmsi][mi] = 1. / lms_arr[lmsi][mi]->value_u();
      }
    }
  }
 
  // Now we've made all our histograms, we'll put them in the TList format
  // of [nominal, syst_1_Up, syst_1_Down, ... , syst_N_Up, syst_N_Down]
  for (unsigned mi = 0; mi < m; ++mi) {
    list_arr[mi] = std::make_shared<TList>();
    for (unsigned xi = 0; xi < (1 + ss * 2); ++xi) {
      list_arr[mi]->Add(hist_arr[mi][xi].get());
    }
  }
 
  // Print the values of the yields and kappa factors that will be inputs
  // to our spline interpolation objects
  if (verbose) {
    for (unsigned mi = 0; mi < m; ++mi) {
      std::cout << boost::format("%-10s") % m_str_vec[mi];
    }
    std::cout << "\n";
    for (unsigned mi = 0; mi < m; ++mi) {
      std::cout << boost::format("%-10.5g") % rate_arr[mi];
    }
    std::cout << "\n";
    for (unsigned ssi = 0; ssi < ss; ++ssi) {
      std::cout << ss_vec[ssi] << " Up" << std::endl;
      for (unsigned mi = 0; mi < m; ++mi) {
        std::cout << boost::format("%-10.5g") % ss_k_hi_arr[ssi][mi];
      }
      std::cout << "\n";
      std::cout << ss_vec[ssi] << " Down" << std::endl;
      for (unsigned mi = 0; mi < m; ++mi) {
        std::cout << boost::format("%-10.5g") % ss_k_lo_arr[ssi][mi];
      }
      std::cout << "\n";
    }
    for (unsigned lmsi = 0; lmsi < lms; ++lmsi) {
      std::cout << lms_vec[lmsi] << " Up" << std::endl;
      for (unsigned mi = 0; mi < m; ++mi) {
        std::cout << boost::format("%-10.5g") % lms_k_hi_arr[lmsi][mi];
      }
      std::cout << "\n";
      std::cout << lms_vec[lmsi] << " Down" << std::endl;
      for (unsigned mi = 0; mi < m; ++mi) {
        std::cout << boost::format("%-10.5g") % lms_k_lo_arr[lmsi][mi];
      }
      std::cout << "\n";
    }
  }
 
  // Can do more sophistical spline interpolation if we want, but let's use
  // simple LINEAR interpolation for now
  TString interp = "LINEAR";
  
  // Here when the force_template_limit option is requested 
  // we have to add some extra terms to the vector of masses to ensure that 
  // the signal pdf goes to 0 outside of the MC template range.
 
  vector<double> new_m_vec(m_vec);
  // Insert an entry at either end of the vector for a mass just slightly
  // outside of the range
  new_m_vec.insert(new_m_vec.begin(),m_vec[0]-1E-6);
  new_m_vec.push_back(m_vec[m-1]+1E-6);
  // Create a corresponding rate array with 0 entries for these new masses
  multi_array<double, 1> new_rate_arr(extents[m+2]);
  new_rate_arr[0] = 0.0;
  for(unsigned i = 0; i < m; ++i) new_rate_arr[i+1] = rate_arr[i] ;
  new_rate_arr[m+1] = 0.0;
 
  if (verbose && force_template_limit) {
     std::cout << ">>>> Forcing rate to 0 outside of template range:" << "\n";
     for(unsigned mi = 0; mi < m+2; ++mi) {
        std::cout << boost::format("%-10.5g") % new_m_vec[mi];
     }
     std::cout << "\n";
     for(unsigned mi = 0; mi < m+2; ++mi) {
        std::cout << boost::format("%-10.5g") % new_rate_arr[mi];
     }
     std::cout << "\n";
  }
  // Create the 1D spline directly from the rate array
  RooSpline1D rate_spline("interp_rate_"+key, "", mass_var, 
                        force_template_limit ? m+2 : m, 
                        force_template_limit ? new_m_vec.data() : m_vec.data(),
                        force_template_limit ? new_rate_arr.data() : rate_arr.data(),
                        interp);
 
  if (file) {
    TGraphErrors tmp(m, m_vec.data(), rate_arr.data(), nullptr, rate_unc_arr.data());
    file->WriteTObject(&tmp, "interp_rate_"+key);
  }
  // Collect all terms that will go into the total normalisation:
  //   nominal * systeff_1 * systeff_2 * ... * systeff_N
  RooArgList rate_prod(rate_spline);
  // For each shape systematic build a 1D spline for kappa_hi and kappa_lo
  for (unsigned ssi = 0; ssi < ss; ++ssi) {
    ss_spl_hi_arr[ssi] = std::make_shared<RooSpline1D>("spline_hi_" +
        key + "_" + ss_vec[ssi], "", mass_var, m, m_vec.data(),
        ss_k_hi_arr[ssi].origin(), interp);
    ss_spl_lo_arr[ssi] = std::make_shared<RooSpline1D>("spline_lo_" +
        key + "_" + ss_vec[ssi], "", mass_var, m, m_vec.data(),
        ss_k_lo_arr[ssi].origin(), interp);
    if (file) {
      TGraph tmp_hi(m, m_vec.data(), ss_k_hi_arr[ssi].origin());
      file->WriteTObject(&tmp_hi, "spline_hi_" + key + "_" + ss_vec[ssi]);
      TGraph tmp_lo(m, m_vec.data(), ss_k_lo_arr[ssi].origin());
      file->WriteTObject(&tmp_lo, "spline_lo_" + key + "_" + ss_vec[ssi]);
    }
    // Then build the AsymPow object for each systematic as a function of the
    // kappas and the nuisance parameter
    ss_asy_arr[ssi] = std::make_shared<AsymPow>("systeff_" +
        key + "_" + ss_vec[ssi], "",
        *(ss_spl_lo_arr[ssi]), *(ss_spl_hi_arr[ssi]),
         *(ss_scale_var_arr[ssi]));
    rate_prod.add(*(ss_asy_arr[ssi]));
  }
  // Same procedure for the lms normalisation systematics: build the splines
  // then the AsymPows and add to the rate_prod list
  for (unsigned lmsi = 0; lmsi < lms; ++lmsi) {
    lms_spl_hi_arr[lmsi] = std::make_shared<RooSpline1D>("spline_hi_" +
        key + "_" + lms_vec[lmsi], "", mass_var, m, m_vec.data(),
        lms_k_hi_arr[lmsi].origin(), interp);
    lms_spl_lo_arr[lmsi] = std::make_shared<RooSpline1D>("spline_lo_" +
        key + "_" + lms_vec[lmsi], "", mass_var, m, m_vec.data(),
        lms_k_lo_arr[lmsi].origin(), interp);
    lms_asy_arr[lmsi] = std::make_shared<AsymPow>("systeff_" +
        key + "_" + lms_vec[lmsi], "", *(lms_spl_lo_arr[lmsi]),
        *(lms_spl_hi_arr[lmsi]), *(lms_var_arr[lmsi]));
    rate_prod.add(*(lms_asy_arr[lmsi]));
  }
  // We'll come back to this rate_prod a bit later.
 
  // Now some fun with binning. We anticipate that the process histograms we
  // have been supplied could have a finer binning than is actually wanted for
  // the analysis (and the fit), in order to avoid known problems with the RMS
  // of a peaking distribution not being morphed smoothly from mass point to
  // mass point if the binning is too wide. The RooMorphingPdf will handle
  // re-binning on the fly, but we have to tell it how to rebin. To do this we
  // assume the observed data histogram has the target binning.
  TH1F data_hist = cb_bp.GetObservedShape();
  TH1F proc_hist = cb_bp.GetShape();
  // The x-axis variable has to be called "CMS_th1x", as this is what
  // text2workspace will use for all the normal processes
  RooRealVar xvar("CMS_th1x", "CMS_th1x", 0,
                 static_cast<float>(data_hist.GetNbinsX()));
  xvar.setBins(data_hist.GetNbinsX());
 
  // Create a second x-axis variable, named specific to the bin, that will be
  // for the finer-binned input
  RooRealVar morph_xvar(("CMS_th1x_"+bin).c_str(), "", 0,
                 static_cast<float>(proc_hist.GetNbinsX()));
  // We're not going to need roofit to evaluate anything as a function of this
  // morphing x-axis variable, so we set it constant
  morph_xvar.setConstant();
  morph_xvar.setBins(proc_hist.GetNbinsX());
  if (verbose) {
    xvar.Print();
    morph_xvar.Print();
  }
 
  // Now we can build the array of vertical-interpolation pdfs (the same that
  // text2workspace builds for every process), one per mass point
  multi_array<std::shared_ptr<FastVerticalInterpHistPdf2>, 1> vpdf_arr(
      extents[m]);
  RooArgList vpdf_list;
 
  TString vert_name = key + "_";
 
  // Follow what ShapeTools.py does and set the smoothing region
  // to the minimum of all of the shape scales
  double qrange = 1.;
  for (unsigned ssi = 0; ssi < ss; ++ssi) {
    if (ss_scale_arr[ssi] < qrange) qrange = ss_scale_arr[ssi];
  }
 
  for (unsigned mi = 0; mi < m; ++mi) {
    // Construct it with the right binning, the right histograms and the right
    // scaling parameters
    vpdf_arr[mi] = std::make_shared<FastVerticalInterpHistPdf2>(
        vert_name + m_str_vec[mi] + "_vmorph", "", morph_xvar, *(list_arr[mi]),
        ss_list, qrange, 0);
    // Add it to a list that we'll supply to the RooMorphingPdf
    vpdf_list.add(*(vpdf_arr[mi]));
  }
  TString morph_name = key + "_morph";
  // At long last, we can build our pdf, giving it:
  //   xvar :       the fixed "CMS_th1x" x-axis variable with uniform integer binning
  //   mass_var:    the floating mass value for the interpolation
  //   vpdf_list:   the list of vertical-interpolation pdfs
  //   m_vec:       the corresponding list of mass points
  //   allow_morph: if false will just evaluate to the closest pdf in mass
  //   data_hist.GetXaxis(): The original (non-uniform) target binning
  //   proc_hist.GetXaxis(): The original (non-uniform) morphing binning
  RooMorphingPdf morph_pdf(morph_name, "", xvar, mass_var, vpdf_list,
                           m_vec, allow_morph, *(data_hist.GetXaxis()),
                           *(proc_hist.GetXaxis()));
  // And we can make the final normalisation product
  // The value of norm_postfix is very important. text2workspace will only look
  // for for a term with the pdf name + "_norm". But it might be the user wants
  // to add even more terms to this total normalisation, so we give them the option
  // of using some other suffix.
  RooProduct morph_rate(morph_name + "_" + TString(norm_postfix), "",
                        rate_prod);
 
  // Dump even more plots
  if (file) MakeMorphDebugPlots(&morph_pdf, &mass_var, m_vec, file, &data_hist);
 
  // Load our pdf and norm objects into the workspace
  ws.import(morph_pdf, RooFit::RecycleConflictNodes());
  ws.import(morph_rate, RooFit::RecycleConflictNodes());
 
  if (!verbose) RooMsgService::instance().setGlobalKillBelow(backup_msg_level);
 
  // Now we can clean up the CH instance a bit
  // We only need one entry for each Process or Systematic now, not one per mass point
  // We'll modify the first mass point to house our new pdf and norm, and drop
  // the rest.
  std::string mass_min = m_str_vec.at(0);
 
  // Dump Process entries for other mass points
  cb.FilterProcs([&](ch::Process const* p) {
    return p->bin() == bin && p->process() == process && p->mass() != mass_min;
  });
  // Dump Systematic entries for other mass points, but only if the type is shape
  // or a lnN that we found varied as a function of the mass. We've already built
  // these uncertainties into our normalisation term. Constant lnN systematics can
  // remain in the datacard and be added by text2workspace
  cb.FilterSysts([&](ch::Systematic const* n) {
    return (n->bin() == bin && n->process() == process) &&
           ((n->mass() != mass_min) || (n->type() == "shape") ||
            (lms_set.count(n->name())));
  });
  // With the remaining Process entry (should only be one if we did this right),
  // Make the mass generic ("*"), drop the TH1 and set the rate to 1.0, as this
  // will now be read from our norm object
  cb.ForEachProc([&](ch::Process * p) {
    if (p->bin() == bin && p->process() == process) {
      p->set_mass("*");
      p->set_shape(nullptr, false);
      p->set_rate(1.0);
    }
  });
  // Just declare the mass to be generic for the remaining systematics
  cb.ForEachSyst([&](ch::Systematic * n) {
    if (n->bin() == bin && n->process() == process) {
      n->set_mass("*");
    }
  });
 
  // And we're done, but note that we haven't populated the Process entry with
  // the PDF or norm objects we created, as we assume the user has more work to
  // do on their RooWorkspace before copying into the CH instance. Once they've
  // imported the workspace:
  //
  //    cb.AddWorkspace(ws);
  //
  // They can populate the Process entries, e.g. if all signals were replaced
  // with Morphing PDFs:
  // 
  //   cb.cp().signals().ExtractPdfs(cb, "htt", "$BIN_$PROCESS_morph");
  return std::string(morph_name);
}
 
void MakeMorphDebugPlots(RooMorphingPdf* pdf, RooAbsReal* mass,
                         std::vector<double> const& masses, TFile* f, TH1 *ref_bins) {
  RooRealVar *rrv = dynamic_cast<RooRealVar*>(mass);
  if (!rrv) return;
  f->cd();
  f->mkdir(pdf->GetName());
  gDirectory->cd(pdf->GetName());
  for (unsigned i = 0; i < masses.size(); ++i) {
    rrv->setVal(masses[i]);
    TH1 * h = pdf->createHistogram("CMS_th1x");
    h->AddDirectory(false);
    TH1 * h2 = nullptr;
    if (ref_bins) { 
      h2 = (TH1*)ref_bins->Clone();
      h2->Reset();
      for (int x = 1; x <= h2->GetNbinsX(); ++x) {
        h2->SetBinContent(x, h->GetBinContent(x));
      }
    } else {
      h2 = h;
      h = nullptr;
    }
    h2->AddDirectory(false);
    h2->SetName(TString::Format("actual_point_%g", masses[i]));
    gDirectory->WriteTObject(h2);
    if (h) delete h;
    if (h2) delete h2;
  }
  double m = masses.front();
  double step = 1;
  if (((masses.back() - masses.front()) / step) > 100.) step = step * 10.;
  if (((masses.back() - masses.front()) / step) > 100.) step = step * 10.;
  if (((masses.back() - masses.front()) / step) < 10.) step = step/10.;
  while (m <= masses.back()) {
    rrv->setVal(m);
    TH1* hm = pdf->createHistogram("CMS_th1x");
    hm->AddDirectory(false);
    TH1 * hm2 = nullptr;
    if (ref_bins) { 
      hm2 = (TH1*)ref_bins->Clone();
      hm2->Reset();
      for (int x = 1; x <= hm2->GetNbinsX(); ++x) {
        hm2->SetBinContent(x, hm->GetBinContent(x));
      }
    } else {
      hm2 = hm;
      hm = nullptr;
    }
    hm2->AddDirectory(false);
    hm2->SetName(TString::Format("morph_point_%g", m));
    //It struggles with m=0, instead adds really small number close to 0
    if(fabs(m) < 1E-5) hm2->SetName(TString::Format("morph_point_0"));
    gDirectory->WriteTObject(hm2);
    if (hm) delete hm;
    if (hm2) delete hm2;
    m += step;
  }
  f->cd();
}
}