Splitter.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
#include <map>
#include <vector>
#include <TF1.h>
 
 
//If compiled within BASF2
#ifdef _PACKAGE_
#include <framework/logging/Logger.h>
#else
#ifndef B2FATAL
#define B2FATAL(arg) { std::cout << arg << std::endl; std::exit(1);}
#define B2ASSERT(str, cond) { if((cond) == false) {std::cout << __FILE__ <<", line "<<__LINE__ << std::endl << "" << #cond << " failed" << std::endl;   std::exit(1);} }
#endif
#define B2INFO(arg) { std::cout << arg << std::endl;}
#define B2WARNING(arg) { std::cout << "WARNING : "<< arg << std::endl;}
#endif
 
 
 
 
namespace Belle2 {
  // The data are divided into many small atoms, where atom has some specific target time length (few mins)
  // but cannot (by definition) span over two runs.
  // The task of the Splitter class is to cluster atoms into the calibration intervals according to th lossFunction
  // In the next step the atoms in each calibration interval are clustered into the calibration subintervals
 
  struct Atom {
    double t1 = 0; 
    double t2 = 0; 
    int nEv = 0; 
  };
 
 
 
  struct ExpRun {
    int exp; 
    int run; 
    ExpRun(int Exp, int Run) : exp(Exp), run(Run) {}
  };
 
 
  struct ExpRunEvt {
    int exp; 
    int run; 
    int evt; 
    ExpRunEvt(int Exp, int Run, int Evt) : exp(Exp), run(Run), evt(Evt) {}
  };
 
 
 
  inline bool operator!=(ExpRun a, ExpRun b) { return (a.exp != b.exp || a.run != b.run); }
 
  inline bool operator<(ExpRun a, ExpRun b) {return ((a.exp < b.exp) || (a.exp == b.exp && a.run < b.run));}
 
  std::map<ExpRun, std::pair<double, double>> filter(const std::map<ExpRun, std::pair<double, double>>& runs, double cut,
                                                     std::map<ExpRun, std::pair<double, double>>& runsRemoved);
 
  std::pair<int, int> getStartEndIndexes(int nIntervals,  std::vector<int> breaks, int indx);
 
 
  inline std::vector<Atom> slice(std::vector<Atom> vec, int s, int e)
  {
    return std::vector<Atom>(vec.begin() + s, vec.begin() + e + 1);
  }
 
 
  std::vector<std::map<ExpRun, std::pair<double, double>>> breaks2intervalsSep(const std::map<ExpRun,
      std::pair<double, double>>& runsMap, const std::vector<Atom>& currVec, const std::vector<int>& breaks);
 
 
 
  class Splitter {
  public:
 
    Splitter() : lossFun(nullptr)  {}
 
 
 
 
    static std::pair<double, double> getStartEnd(std::vector<std::map<ExpRun, std::pair<double, double>>> res)
    {
      return {res.front().begin()->second.first,
              res.back().rbegin()->second.second};
    }
 
 
    static std::vector<double> getBreaks(std::vector<std::map<ExpRun, std::pair<double, double>>> res)
    {
      std::vector<double> breaks;
      for (int k = 0; k < int(res.size()) - 1; ++k) {
        double e = res.at(k).rbegin()->second.second; //end of the previous interval
        double s = res.at(k + 1).begin()->second.first; //start of the next interval
        breaks.push_back((e + s) / 2.);  //the break time is in between
      }
      return breaks;
    }
 
 
    static std::map<ExpRun, std::pair<double, double>> mergeIntervals(std::map<ExpRun, std::pair<double, double>> I1,
                                                    std::map<ExpRun, std::pair<double, double>> I2);
 
 
    template<typename Evt>
    std::vector<std::vector<std::map<ExpRun, std::pair<double, double>>>>  getIntervals(
      const std::map<ExpRun, std::pair<double, double>>& runs,  std::vector<Evt> evts,
      TString lossFunctionOuter, TString lossFunctionInner, double atomSize = 3. / 60)
    {
      //sort events by time
      std::sort(evts.begin(), evts.end(), [](const Evt & e1, const Evt & e2) { return (e1.t > e2.t); });
 
      // Divide into small intervals
      std::vector<std::pair<double, double>> smallRuns = splitToSmall(runs, atomSize);
 
      std::vector<Atom> atoms = createAtoms(smallRuns, evts);
 
 
      // Calculate breaks for the calib. intervals
      lossFun = toTF1(lossFunctionOuter);
      std::vector<int> breaksSize = dynamicBreaks(atoms);
 
      std::vector<std::vector<std::map<ExpRun, std::pair<double, double>>>> splitsRun; //split intervals with runNumber info
 
 
      // Calculate breaks for the calib. subintervals
      lossFun = toTF1(lossFunctionInner);
      std::vector<int> breaksVtx;
      for (int i = 0; i < int(breaksSize.size()) + 1; ++i) { //loop over all calib. intervals
        int s, e;
        std::tie(s, e) = getStartEndIndexes(atoms.size(), breaksSize, i);
 
        // Store only atoms lying in the current calib. interval
        auto currVec = slice(atoms, s, e);
 
        // Get optimal breaks for particular calib. interval
        std::vector<int> breaksSmall = dynamicBreaks(currVec);
 
        // Incorporate the runNumber information
        auto splitSeg = breaks2intervalsSep(runs, currVec, breaksSmall);
 
        // Put the vector of subintervals (representing division of the interval)
        // into the vector
        splitsRun.push_back(splitSeg);
 
        if (s != 0) breaksVtx.push_back(s - 1);
        for (auto b : breaksSmall)
          breaksVtx.push_back(b + s);
        if (e != int(breaksSize.size())) breaksVtx.push_back(e);
      }
 
      return splitsRun;
    }
 
 
  private:
 
    std::vector<int> dynamicBreaks(const std::vector<Atom>&  runs);
 
 
 
    double getMinLoss(const std::vector<Atom>&  vec,   int e, std::vector<int>& breaks);
 
 
    double lossFunction(const std::vector<Atom>&  vec, int s, int e) const;
 
 
    static std::vector<std::pair<double, double>> splitToSmall(std::map<ExpRun, std::pair<double, double>> runs,
                                                               double intSize = 1. / 60);
 
    template<typename Evt>
    inline std::vector<Atom> createAtoms(const std::vector<std::pair<double, double>>& atomsTimes,  const std::vector<Evt>& evts)
    {
      std::vector<Atom> atoms(atomsTimes.size());
 
      // Store start/end times to atoms
      for (unsigned a = 0; a < atoms.size(); ++a) {
        atoms[a].t1 = atomsTimes[a].first;
        atoms[a].t2 = atomsTimes[a].second;
      }
 
 
      // count the number of events in each atom
      for (unsigned i = 0; i < evts.size(); ++i) {
        for (unsigned a = 0; a < atoms.size(); ++a)
          if (atoms[a].t1 <= evts[i].t  && evts[i].t < atoms[a].t2)
            ++atoms[a].nEv;
 
      }
 
 
 
      return atoms;
    }
 
 
    TF1* toTF1(TString LossString)
    {
      LossString.ReplaceAll("rawTime", "[0]");
      LossString.ReplaceAll("netTime", "[1]");
      LossString.ReplaceAll("maxGap",  "[2]");
      LossString.ReplaceAll("nEv",     "[3]");
 
      return (new TF1("lossFun", LossString));
    }
 
 
    TF1*   lossFun; 
 
 
    std::vector<std::pair<double, std::vector<int>>> cache;
 
 
  };
 
 
  template<typename Evt>
  inline std::map<ExpRun, std::pair<double, double>> getRunInfo(const std::vector<Evt>& evts)
  {
    std::map<ExpRun, std::pair<double, double>> runsInfo;
 
    for (auto& evt : evts) {
      int Exp = evt.exp;
      int Run = evt.run;
      double time = evt.t;
      if (runsInfo.count(ExpRun(Exp, Run))) {
        double tMin, tMax;
        std::tie(tMin, tMax) = runsInfo.at(ExpRun(Exp, Run));
        tMin = std::min(tMin, time);
        tMax = std::max(tMax, time);
        runsInfo.at(ExpRun(Exp, Run)) = {tMin, tMax};
      } else {
        runsInfo[ExpRun(Exp, Run)] = {time, time};
      }
 
    }
    return runsInfo;
  }
 
 
  template<typename Evt>
  inline ExpRunEvt getPosition(const std::vector<Evt>& events, double tEdge)
  {
    ExpRunEvt evt(-1, -1, -1);
    double tBreak = -1e10;
    for (auto& e : events) {
      if (e.t < tEdge) {
        if (e.t > tBreak) {
          tBreak = e.t;
          evt =  ExpRunEvt(e.exp, e.run, e.evtNo);
        }
      }
    }
    return evt;
  }
 
// convert splitPoints [in UTC time] to expRunEvt
  template<typename Evt>
  std::vector<ExpRunEvt> convertSplitPoints(const std::vector<Evt>& events, std::vector<double> splitPoints)
  {
 
    std::vector<ExpRunEvt>  breakPos;
    for (auto p : splitPoints) {
      auto pos = getPosition(events, p);
      breakPos.push_back(pos);
    }
    return breakPos;
  }
 
 
}