finetune.h

Go to the documentation of this file.

#pragma once
#include "alcor.h"
 
enum type_t
{
  alcor_hit = 1,
  trigger_tag = 9,
  start_spill = 7,
  end_spill = 15
};
 
//  Addition of fine tune class
namespace alcor
{
  //  [0] Slope,  [1] Intercept
  typedef std::array<float, 2> calibration_unit;
 
  //  Class for the fine calibration
  class finecalibration
  {
  public:
    //  --- Variables ---
    //  --- cindex: tdc + eo_channel*4 + chip*4*32
    //  calibration_scheme[device][cindex]  [0] Intercept
    //                                      [1] Slope
    const std::vector<std::string> calibration_csv_fields = {"device", "cindex", "Intercept", "Slope"};
    std::map<std::string, std::map<int, calibration_unit>> calibration_database;
 
    //  --- Generators ---
    finecalibration() {};
    finecalibration(std::string filename) { read_calibration(filename); };
 
    //  --- Methods ---
    //  --- Getters
    float get_intercept(std::string device, int cindex) { return calibration_database[device][cindex][0]; }
    float get_slope(std::string device, int cindex) { return calibration_database[device][cindex][1]; }
    //  --- Setters
    void set_intercept(std::string device, int cindex, float intercept) { calibration_database[device][cindex][0] = intercept; }
    void set_slope(std::string device, int cindex, float slope) { calibration_database[device][cindex][1] = slope; }
    //  --- Calculate phase
    float get_phase(std::string device, int cindex, float fine) { return fine * get_slope(device, cindex) + get_intercept(device, cindex); }
    //  --- Calculate correction parameters from pulsed data
    std::vector<float> get_guesstimate_correction_from_pulser(alcor::t_hit previous_hit, alcor::t_hit current_hit);
    void calculate_correction_from_pulser(std::string filename, std::string device, int cindex);
    //  --- I/O calibration scheme from file
    void read_calibration(std::string filename);
    void write_calibration(std::string filename);
 
  private:
  };
}
 
std::map<int, std::array<float, 2>> get_guesstimate_correction_from_pulser(Tree *data_tree, alcor::t_hit &current_hit, std::vector<int> cindex_target = {0, 1, 2, 3})
{
  //  Run on data tree
  std::map<int, std::pair<int, std::vector<alcor::t_hit>>> calibration_map;
  std::map<int, std::array<float, 2>> result_map;
  for (int iev = 0; iev < nev; ++iev)
  {
    data_tree->GetEntry(iev);
    auto current_cindex = current_hit.tdc + current_hit.channel * 4 + ((int)(current_hit.fifo / 4)) * 4 * 32;
    if (!std::find(cindex_target.begin(), cindex_target.end(), current_cindex) != cindex_target.end())
      continue;
    calibration_map[current_cindex].push_back({iev, current_hit});
  }
  return result_map;
}
 
void alcor::finecalibration::calculate_correction_from_pulser(std::string filename, std::string device, int cindex)
{
  //  Load pulsed signals from data file
  alcor::t_hit current_hit;
  alcor::t_hit previous_hit;
  auto data_tree = load_alcor_data_tree(filename, current_hit);
  auto nev = data_tree->GetEntries();
 
  //  Minimisation function
  auto periodic_sigma_chi2 = [&](const double *parameters)
  {
    //  Define result
    double current_chi2 = 0.;
 
    //  Do the analysis spill by spill
    bool has_first_hit_in_spill = false;
    int current_spill_contributors = 0;
 
    //  Run on data tree
    for (int iev = 0; iev < nev; ++iev)
    {
      data_tree->GetEntry(iev);
 
      //  Start of spill reset
      if (current_hit.type == start_spill)
      {
        has_first_hit_in_spill = false;
        continue;
      }
 
      //  Operate only if alcro hit events
      if (current_hit.type != alcor_hit)
        continue;
 
      //  Select the requested channel
      if (current_hit.column * 4 + current_hit.pixel != (int)(cindex / 4))
        continue;
 
      //  Calculate corrected time
      int coarse = current_hit.coarse + alcor::rollover_to_coarse * current_hit.rollover;
      double current_correction = parameters[current_hit.tdc] + (double)current_hit.fine * parameters[current_hit.tdc + 4];
      current_hit.time = (double)coarse - current_correction;
 
      if (!has_first_hit_in_spill)
      {
        previous_hit = current_hit;
        has_first_hit_in_spill = true;
        continue;
      }
 
      //  Compute the delta of times
      double current_delta = current_hit.time - previous_hit.time - parameters[8];
 
      // this is a safety measure, because sometimes there is some crap
      if (fabs(coarse - (previous_hit.coarse + alcor::rollover_to_coarse * previous_hit.rollover) - parameters[8]) < 5)
      {
        current_chi2 += current_delta * current_delta;
        ++current_spill_contributors;
      }
 
      //  Assign current to previous hit and loop continues
      previous_hit = current_hit;
    }
    
    return current_chi2;
  };
 
  ROOT::Math::Functor minimisation_function(periodic_sigma_chi2, 9);
  ROOT::Fit::Fitter fitter;
  double pStart[9] = {0.5, 0.5, 0.5, 0.5, 0.0156, 0.0156, 0.0156, 0.0156, 320.};
  fitter.SetFCN(minimisation_function, pStart);
  fitter.Config().ParSettings(0).SetName("off_0");
  fitter.Config().ParSettings(1).SetName("off_1");
  fitter.Config().ParSettings(2).SetName("off_2");
  fitter.Config().ParSettings(3).SetName("off_3");
  fitter.Config().ParSettings(4).SetName("iif_0");
  fitter.Config().ParSettings(5).SetName("iif_1");
  fitter.Config().ParSettings(6).SetName("iif_2");
  fitter.Config().ParSettings(7).SetName("iif_3");
  fitter.Config().ParSettings(8).SetName("period");
  fitter.Config().ParSettings(0).Fix();
 
  bool ok = fitter.FitFCN();
  auto result = fitter.Result();
  result.Print(std::cout);
  double par[9], pare[9];
  auto hParam = new TH1F("hParam", "", 9, 0., 9.);
  for (int ipar = 0; ipar < 9; ++ipar)
  {
    par[ipar] = result.Parameter(ipar);
    pare[ipar] = result.ParError(ipar);
  }
};
 
void alcor::finecalibration::read_calibration(std::string filename)
{
  std::ifstream data_stream(filename);
  std::string current_line;
  while (std::getline(data_stream, current_line))
  {
    //  Skip comment characters
    if (current_line[0] == '#' || current_line[0] == ' ')
      continue;
    //  Read database
    std::stringstream string_in_stream(current_line);
    std::map<std::string, std::string> data_by_field;
    std::string current_data;
    for (auto current_field : alcor::finecalibration::calibration_csv_fields)
    {
      string_in_stream >> current_data;
      data_by_field[current_field] = current_data;
    }
    this->set_intercept(data_by_field["device"], std::stoi(data_by_field["cindex"]), std::stof(data_by_field["Intercept"]));
    this->set_slope(data_by_field["device"], std::stoi(data_by_field["cindex"]), std::stof(data_by_field["Slope"]));
  }
};
 
void alcor::finecalibration::write_calibration(std::string filename)
{
  std::ofstream data_stream(filename);
  data_stream << "#\tCalibration file for Fine Tune\n";
  auto current_time = std::chrono::system_clock::now();
  std::time_t current_time_stamp = std::chrono::system_clock::to_time_t(current_time);
  data_stream << "#\tGenerated on " << std::ctime(&current_time_stamp);
  data_stream << "#\t## Device ## Calibration index ## Intercept ## Slope ##\n";
  for (auto [device, calibration_of_device] : calibration_database)
  {
    for (auto [cindex, current_calibration_unit] : calibration_of_device)
    {
      data_stream << device << "\t" << cindex << "\t" << current_calibration_unit[0] << "\t" << current_calibration_unit[1] << "\n";
    }
  }
};