/src/methodComp.cpp

/**
 * methodComp.cpp - Compare PHPI method with the peak fitting method.
 *
 * author: Christoph Weinsheimer
 * email: christoph.weinsheimer@desy.de
 */
#include <iostream>

#include <fast/utils.h>
#include <fast/processor.h>
#include <fast/pulsefct.h>

#include "TFile.h"
#include "TMath.h"
#include "TTree.h"

#include "version.h"


using namespace fast;
using ConfigParser = fast::utils::ConfigParser;


int main(int argc, char* argv[]) {
  common::printInfo(argc, argv);
  utils::printLibVersion();
  utils::Logger::instance()->setLevel(utils::LogLevel::INFO);

  if (argc != 4) {
    std::cout << "Syntax: methodComp <input> <config> <output>" << std::endl;
    return 1;
  }

  // Parse configuration
  auto config = ConfigParser::parseJson(argv[2]);
  const Channel ch = Channel::ChA;

  // Read input file
  TFile infile(argv[1], "READ");
  std::shared_ptr<TTree> inTree((TTree*) infile.Get("alazar_data"));

  // Create output file
  TFile outfile(argv[3], "RECREATE");
  auto outTree = std::make_shared<TTree>("processed_data", "processed_data");
  auto plotDir = outfile.mkdir("plots");
  plotDir->cd();

  // Create peakfinding unit
  std::shared_ptr<PeakFindingUnit> pfu = nullptr;
  double threshold = config["threshold"]["value"];
  auto noise_ref = ConfigParser::parseJson(
    utils::expandEnv(config["threshold"]["reference"]));
  std::stringstream param_key;
  param_key << std::fixed << std::setprecision(2);

  switch (config["method"].get<unsigned int>()) {
    case PeakFindingUnit::DERIV:
      param_key << config["params"]["cutOff"].get<float>();
      threshold *= noise_ref["deriv"][param_key.str()]["sigma"].get<double>();
      pfu = std::make_shared<DerivProcessor>(
        "L2", ch, config["params"]["cutOff"], threshold
      );
      break;
    case PeakFindingUnit::XCORR:
      param_key << config["params"]["tMinFactor"].get<float>();
      threshold *= noise_ref["xcorr"][param_key.str()]["sigma"].get<double>();
      pfu = std::make_shared<XCorrProcessor>(
        "L2", ch, PulseShape::PHOTON_55mV(), config["params"]["tMinFactor"],
        threshold
      );
      break;
  }
  pfu->setTree(outTree);

  // Create Processors
  PulseFitProcessor pu_fit("fit", ch, 0.005, pfu);
  pu_fit.setTree(outTree);
  PHPIProcessor     pu_phpi("phpi", ch, true);
  pu_phpi.setTree(outTree);

  // Create TreeRecord for truth record part
  struct TruthRecord : public TreeRecord {
    double truth_amplitude;
    double truth_tRise;
    double truth_tFall;
    double truth_tOffset;
    BranchMap branchmap{
      {"truth_amplitude_0", &truth_amplitude},
      {"truth_tRise_0",     &truth_tRise},
      {"truth_tFall_0",     &truth_tFall},
      {"truth_tOffset_0",   &truth_tOffset},
    };
  };

  // Create TreeReader
  TreeReader<TruthRecord>  truthReader(inTree);
  TreeReader<AlazarRecord> treeReader(inTree);

  // Create additional branches
  double fit_integral;
  double fit_dAmplitude, fit_dIntegral, phpi_dAmplitude, phpi_dIntegral;
  outTree->Branch("fit_integral",    &fit_integral);
  outTree->Branch("fit_dAmplitude",  &fit_dAmplitude);
  outTree->Branch("fit_dIntegral",   &fit_dIntegral);
  outTree->Branch("phpi_dAmplitude", &phpi_dAmplitude);
  outTree->Branch("phpi_dIntegral",  &phpi_dIntegral);

  // Event loop
  for (size_t i = 0; i < treeReader.size(); ++i) {
    auto event = treeReader.getEntry(i);
    auto truth = truthReader.getEntry(i);

    auto fitRecord  = pu_fit(event);
    auto phpiRecord = pu_phpi(event);

    fit_dAmplitude  = (fitRecord.amplitude[0] - truth.truth_amplitude)
                      / truth.truth_amplitude;
    phpi_dAmplitude = (phpiRecord.ph - truth.truth_amplitude)
                      / truth.truth_amplitude;

    auto truth_integral = pulsefct::pulse_integral(truth.truth_amplitude,
                                                   truth.truth_tRise,
                                                   truth.truth_tFall);
    fit_integral  = pulsefct::pulse_integral(fitRecord.amplitude[0],
                                             fitRecord.tRise[0],
                                             fitRecord.tFall[0]);

    fit_dIntegral  = (fit_integral - truth_integral) / truth_integral;
    phpi_dIntegral = (phpiRecord.pi - truth_integral) / truth_integral;

    if (!fitRecord.fitResult->IsValid()) {
      utils::LogWarning("methodComp")
        << "Fit failed! Event " + std::to_string(i);
      fitRecord.graph->Write(utils::add_suffix("failed", i).c_str());
    }

    outTree->Fill();
  }

  outfile.cd();
  outTree->Write();

  return 0;
}