eclAutocovarianceCalibrationC1Algorithm.cc

/**************************************************************************
 * 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.                  *
 **************************************************************************/
 
/* Own header. */
#include <ecl/calibration/eclAutocovarianceCalibrationC1Algorithm.h>
 
/* ECL headers. */
#include <ecl/dataobjects/ECLElementNumbers.h>
#include <ecl/dbobjects/ECLCrystalCalib.h>
 
/* ROOT headers. */
#include <TFile.h>
#include <TGraph.h>
#include <TH2F.h>
 
using namespace Belle2;
using namespace ECL;
 
eclAutocovarianceCalibrationC1Algorithm::eclAutocovarianceCalibrationC1Algorithm():
  CalibrationAlgorithm("eclAutocovarianceCalibrationC1Collector")
{
  setDescription(
    "Determine noise threshold for waveforms to be used in computing the covariance matrix"
  );
}
 
CalibrationAlgorithm::EResult eclAutocovarianceCalibrationC1Algorithm::calibrate()
{
 
  gROOT->SetBatch();
 
  auto PPVsCrysID = getObjectPtr<TH2F>("PPVsCrysID");
 
  std::vector<float> cellIDs;
  std::vector<float> PPamps;
  std::vector<float> PPampsErrorVector;
 
  for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
 
    TH1F* hPP = (TH1F*)PPVsCrysID->ProjectionY("hPP", crysID + 1, crysID + 1);
 
    int Total = hPP->GetEntries();
    float fraction = 0.0;
    int counter = 0;
 
    if (Total < m_TotalCountsThreshold) {
      B2INFO("eclAutocovarianceCalibrationC1Algorithm: warning total entries for cell ID " << crysID + 1 << " is only: " << Total <<
             " Requirement is m_TotalCountsThreshold: " << m_TotalCountsThreshold);
      return c_NotEnoughData;
    } else {
      int subTotal = 0;
      while (fraction < m_lowestEnergyFraction) {
        subTotal += hPP->GetBinContent(counter);
        fraction = ((float)subTotal) / ((float)Total);
        counter++;
      }
    }
 
    cellIDs.push_back(crysID + 1);
    PPamps.push_back(counter);
    PPampsErrorVector.push_back(0);
 
    B2INFO("eclAutocovarianceCalibrationC1Algorithm: crysID counter fraction Total " << crysID << " " << counter << " " << fraction <<
           " " << Total);
 
  }
 
  auto gPPVsCellID = new TGraph(cellIDs.size(), cellIDs.data(), PPamps.data());
  gPPVsCellID->SetName("gPPVsCellID");
 
  TString fName = m_outputName;
  TDirectory::TContext context;
  TFile* histfile = new TFile(fName, "recreate");
  histfile->cd();
  PPVsCrysID->Write();
  gPPVsCellID->Write();
  histfile->Close();
  delete histfile;
 
  ECLCrystalCalib* PPThreshold = new ECLCrystalCalib();
  PPThreshold->setCalibVector(PPamps, PPampsErrorVector);
  saveCalibration(PPThreshold, "ECLAutocovarianceCalibrationC1Threshold");
  B2INFO("eclAutocovarianceCalibrationC1Algorithm: successfully stored ECLAutocovarianceCalibrationC1Threshold constants");
 
  return c_OK;
}