MdADCCalibration.cc

Go to the documentation of this file.

#include <utl/ErrorLogger.h>
#include <utl/UTMPoint.h>
#include <fwk/CentralConfig.h>

#include <evt/Event.h>
#include <mevt/MEvent.h>
#include <mevt/Counter.h>
#include <mevt/Module.h>
#include <utl/Branch.h>
#include <utl/ConfigUtils.h>

#include "MdADCCalibration.h"

using namespace fwk;
using namespace utl;
using namespace MdADCCalibrationAG;

namespace MdADCCalibrationAG {

  VModule::ResultFlag
  MdADCCalibration::Init()
{
  // Configuration reading
  utl::Branch config = fwk::CentralConfig::GetInstance()->GetTopBranch("MdADCCalibration");

  LoadConfig(config, "useOfflineCalibration", fUseOfflineCalibration, 0);
  LoadConfig(config, "useOnlineCalibration", fUseOnlineCalibration, 1);

  // Fill the fUseOfflineCalibration map from the configuration file
  if (fUseOfflineCalibration) {
    utl::Branch calibrationLG = config.GetChild("moduleOfflineCalibrationLG");
    for (utl::Branch b = calibrationLG.GetFirstChild(); b; b = b.GetNextSibling()) {
      utl::AttributeMap attibutes = b.GetAttributes();
      const int counterId = det::VManager::FindComponent<unsigned int>("counterId", attibutes);
      const int moduleId = det::VManager::FindComponent<unsigned int>("moduleId", attibutes);
      const int globalModuleId = counterId*1000 + moduleId;
      const double calibration = b.Get<double>();
      fModuleOfflineCalibrationLG[globalModuleId] = calibration;
    }

    utl::Branch calibrationHG = config.GetChild("moduleOfflineCalibrationHG");
    for (utl::Branch b = calibrationHG.GetFirstChild(); b; b = b.GetNextSibling()) {
      utl::AttributeMap attibutes = b.GetAttributes();
      const int counterId = det::VManager::FindComponent<unsigned int>("counterId", attibutes);
      const int moduleId = det::VManager::FindComponent<unsigned int>("moduleId", attibutes);
      const int globalModuleId = counterId*1000 + moduleId;
      const double calibration = b.Get<double>();
      fModuleOfflineCalibrationHG[globalModuleId] = calibration;
    }

    utl::Branch resolutionLG = config.GetChild("moduleOfflineResolutionLG");
    for (utl::Branch b = resolutionLG.GetFirstChild(); b; b = b.GetNextSibling()) {
      utl::AttributeMap attibutes = b.GetAttributes();
      const int counterId = det::VManager::FindComponent<unsigned int>("counterId", attibutes);
      const int moduleId = det::VManager::FindComponent<unsigned int>("moduleId", attibutes);
      const int globalModuleId = counterId*1000 + moduleId;
      const double resolution = b.Get<double>();
      fModuleOfflineResolutionLG[globalModuleId] = resolution;
    }

    utl::Branch resolutionHG = config.GetChild("moduleOfflineResolutionHG");
    for (utl::Branch b = resolutionHG.GetFirstChild(); b; b = b.GetNextSibling()) {
      utl::AttributeMap attibutes = b.GetAttributes();
      const int counterId = det::VManager::FindComponent<unsigned int>("counterId", attibutes);
      const int moduleId = det::VManager::FindComponent<unsigned int>("moduleId", attibutes);
      const int globalModuleId = counterId*1000 + moduleId;
      const double resolution = b.Get<double>();
      fModuleOfflineResolutionHG[globalModuleId] = resolution;
    }

  }
  return eSuccess;
}

  VModule::ResultFlag
  MdADCCalibration::Run(evt::Event& event)
  {

  if (!event.HasMEvent()) {
    WARNING("\n+++++++++\nEvent without MEvent: nothing to be done. Skipping to next event.\n++++++++++\n");
    return eContinueLoop;
  }

  INFO("Starting");
  std::ostringstream os;

  mevt::MEvent& me = event.GetMEvent();

  for (mevt::MEvent::CounterIterator ic = me.CountersBegin(), ec = me.CountersEnd(); ic != ec; ++ic) {

    mevt::Counter& counter = *ic;
    const int counterId = counter.GetId();

    if (counter.IsRejected()) {
      os.str("");
      os << "Skipping the rejected counter " << counterId << " .";
      INFO(os);
      continue;
    }

    os.str("");
    os << "Processing counter " << counterId;
    INFO(os);


    for (mevt::Counter::ModuleIterator im = counter.ModulesBegin(), em = counter.ModulesEnd(); im != em; ++im) {

      DEBUGLOG("Starting module loop");
      mevt::Module& module = *im;
      const int moduleId = module.GetId();

      os.str("");
      os << "Processing module " << moduleId;
      DEBUGLOG(os);


      if (module.IsRejected()) {
        os.str("");
        os << "Module " << module.GetId() << " of counter " << counterId;
        os << " discarded because flagged as: " << module.GetRejectionReason();
        WARNING(os);
        continue;
      }

      // Load calibration into the event
      LoadModuleCalibration(counter, module);

    } // end module loop
  } // end counter loop
  return eSuccess;
}


VModule::ResultFlag
MdADCCalibration::Finish()
{
  return eSuccess;
}


void MdADCCalibration::LoadModuleCalibration(mevt::Counter& counter, mevt::Module& module)
{


 double calibrationA = 0.0;
 double calibrationB = 0.0;
 double resolutionA = 0.0;
 double resolutionB = 0.0;


 if(fUseOnlineCalibration)
 {
     GetOnlineCalibration(calibrationA, resolutionA, calibrationB, resolutionB, module);
 }

 if (calibrationA == 0.0 && fUseOfflineCalibration)
 {
     GetOfflineCalibration(calibrationA, resolutionA, calibrationB, resolutionB, counter, module);

 }

  // Set results in event
  mevt::ModuleRecData& mRecData = module.GetRecData();

  mRecData.SetMeanChargeMuonHG(calibrationB);
  mRecData.SetStdDevChargeMuonHG(resolutionB);
  mRecData.SetMeanChargeMuonLG(calibrationA);
  mRecData.SetStdDevChargeMuonLG(resolutionA);



  std::ostringstream os;
  os.str("");
  os << "Mean charge muon ADC LG " << calibrationA << ", StdDev  = " << resolutionA
     << "Mean charge muon ADC HG " << calibrationB << ", StdDev  = " << resolutionB;
  INFO(os);
}

void MdADCCalibration::GetOfflineCalibration(double& calibrationA, double& resolutionA, double& calibrationB, double& resolutionB, mevt::Counter& counter, mevt::Module& module)
{
    std::ostringstream os;
    const int globalModuleId = counter.GetId()*1000 + module.GetId();
    if(fModuleOfflineCalibrationLG.find(globalModuleId) != fModuleOfflineCalibrationLG.end())
    {
        calibrationA = fModuleOfflineCalibrationLG[globalModuleId];
    }else
    {
        calibrationA = fModuleOfflineCalibrationLG[0];
        os.str("");
        os << "Module " << module.GetId() << " of counter " << counter.GetId();
        os << " uses global calibration";
        WARNING(os);
    }

    if(fModuleOfflineResolutionLG.find(globalModuleId) != fModuleOfflineResolutionLG.end())
    {
        resolutionA = fModuleOfflineResolutionLG[globalModuleId];
    }else
    {
        resolutionA = fModuleOfflineResolutionLG[0];
    }

    if(fModuleOfflineCalibrationHG.find(globalModuleId) != fModuleOfflineCalibrationHG.end())
    {
        calibrationB = fModuleOfflineCalibrationHG[globalModuleId];
    }else
    {
        calibrationB = fModuleOfflineCalibrationHG[0];
    }

    if(fModuleOfflineResolutionHG.find(globalModuleId) != fModuleOfflineResolutionHG.end())
    {
        resolutionB = fModuleOfflineResolutionHG[globalModuleId];
    }else
    {
        resolutionB = fModuleOfflineResolutionHG[0];
    }
}

void MdADCCalibration::GetOnlineCalibration(double& calibrationA, double& resolutionA, double& calibrationB, double& resolutionB, mevt::Module& module)
{

    if(module.GetId()>0){
        calibrationA = 0.0;
        resolutionA = 0.0;
        calibrationB = 0.0;
        resolutionB = 0.0;
    }


    /*
     * This method has to be implemented. The integration of the online calibration into CDAS is not yet ready, so I do not know what kind of data
     * I will be receiving in this instance.
     *
     * The noise histogram has to be subtracted from the charge histogram to obtain the calibration histogram from which the mean value
     * and StdDev give the calibration and resolution. There are 6 histograms in total, 3 for the low gain channel and 3 for the high gain channel
     *
     *
     * You can check GAP_2019_046 to see how this is performed,
     *
     * */

    /*const int moduleId = module.GetId();
    const vector<CalibrationData> dataCalib = module.GetCalibrationData()
    if(dataCalib.size() > 0)
    {
        int nroBins = 1000;
        double xminA = -1000;
        double xmaxA = 3000;
        double binSizeA = (xmaxA-xminA)/(double) nroBins;

        double xminB = -6000;
        double xmaxB = 10000;
        double binSizeB = (xmaxB-xminB)/(double) nroBins;


        int histogramChargeA[nroBins];
        int histogramChargeB[nroBins];
        int histogramNoiseA[nroBins];
        int histogramNoiseB[nroBins];
        int histogramCalibrationA[nroBins];
        int histogramCalibrationB[nroBins];

        for(calib_iterator calibIt = dataCalib.begin() ; calibIt != dataCalib.end() ; ++calibIt)
        {


        }
    }

    */


}




} // end namespace