eclCosmicECollectorModule.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/modules/eclCosmicECollector/eclCosmicECollectorModule.h>
 
/* ECL headers. */
#include <ecl/dataobjects/ECLDigit.h>
#include <ecl/dataobjects/ECLElementNumbers.h>
#include <ecl/dbobjects/ECLCrystalCalib.h>
#include <ecl/geometry/ECLGeometryPar.h>
#include <ecl/geometry/ECLNeighbours.h>
 
/* Basf2 headers. */
#include <framework/dataobjects/EventMetaData.h>
#include <trg/ecl/TrgEclMapping.h>
 
/* ROOT headers. */
#include <TH2F.h>
 
using namespace std;
using namespace Belle2;
using namespace ECL;
 

//                 Register the Module
REG_MODULE(eclCosmicECollector);

//                 Implementation

 
eclCosmicECollectorModule::eclCosmicECollectorModule() : CalibrationCollectorModule(), m_ECLExpCosmicESame("ECLExpCosmicESame"),
  m_ECLExpCosmicEDifferent("ECLExpCosmicEDifferent"), m_ElectronicsCalib("ECLCrystalElectronics"),
  m_CosmicECalib("ECLCrystalEnergyCosmic")
{
  setDescription("Calibration Collector Module for ECL single crystal energy calibration using cosmic rays");
  setPropertyFlags(c_ParallelProcessingCertified);
  addParam("minCrysE", m_minCrysE, "Minimum energy in GeV for a crystal to be considered hit", 0.010);
  addParam("mockupL1", m_mockupL1, "Calculate energy per trigger cell in lieu of trigger simulation", false);
  addParam("trigThreshold", m_trigThreshold, "Minimum energy in GeV per trigger cell to mock up L1 trigger", 0.1);
}
 

void eclCosmicECollectorModule::prepare()
{
  B2INFO("eclCosmicECollector: Experiment = " << m_evtMetaData->getExperiment() << "  run = " << m_evtMetaData->getRun());

  m_eclDigitArray.isRequired();

  auto EnvsCrysSameRing = new TH2F("EnvsCrysSameRing", "Normalized energy vs crystal ID, same theta ring;crystal ID;E/Expected",
                                   ECLElementNumbers::c_NCrystals,
                                   0, ECLElementNumbers::c_NCrystals, 99, .025, 2.5);
  registerObject<TH2F>("EnvsCrysSameRing", EnvsCrysSameRing);
 
  auto EnvsCrysDifferentRing = new TH2F("EnvsCrysDifferentRing",
                                        "Normalized energy vs crystal ID, different theta rings;crystal ID;E/Expected", ECLElementNumbers::c_NCrystals, 0,
                                        ECLElementNumbers::c_NCrystals, 99, .025, 2.5);
  registerObject<TH2F>("EnvsCrysDifferentRing", EnvsCrysDifferentRing);
 
  auto ExpEvsCrysSameRing = new TH1F("ExpEvsCrysSameRing",
                                     "Sum expected energy vs crystal ID, same theta ring;crystal ID;Energy (GeV)", ECLElementNumbers::c_NCrystals, 0,
                                     ECLElementNumbers::c_NCrystals);
  registerObject<TH1F>("ExpEvsCrysSameRing", ExpEvsCrysSameRing);
 
  auto ExpEvsCrysDifferentRing = new TH1F("ExpEvsCrysDifferentRing",
                                          "Sum expected energy vs crystal ID, different theta rings;crystal ID;Energy (GeV)", ECLElementNumbers::c_NCrystals, 0,
                                          ECLElementNumbers::c_NCrystals);
  registerObject<TH1F>("ExpEvsCrysDifferentRing", ExpEvsCrysDifferentRing);
 
  auto ElecCalibvsCrysSameRing = new TH1F("ElecCalibvsCrysSameRing",
                                          "Sum electronics calib const vs crystal ID, same theta ring;crystal ID", ECLElementNumbers::c_NCrystals, 0,
                                          ECLElementNumbers::c_NCrystals);
  registerObject<TH1F>("ElecCalibvsCrysSameRing", ElecCalibvsCrysSameRing);
 
  auto ElecCalibvsCrysDifferentRing = new TH1F("ElecCalibvsCrysDifferentRing",
                                               "Sum electronics calib const vs crystal ID, different theta rings;crystal ID", ECLElementNumbers::c_NCrystals, 0,
                                               ECLElementNumbers::c_NCrystals);
  registerObject<TH1F>("ElecCalibvsCrysDifferentRing", ElecCalibvsCrysDifferentRing);
 
  auto InitialCalibvsCrysSameRing = new TH1F("InitialCalibvsCrysSameRing",
                                             "Sum initial cosmic calib const vs crystal ID, same theta ring;crystal ID", ECLElementNumbers::c_NCrystals, 0,
                                             ECLElementNumbers::c_NCrystals);
  registerObject<TH1F>("InitialCalibvsCrysSameRing", InitialCalibvsCrysSameRing);
 
  auto InitialCalibvsCrysDifferentRing = new TH1F("InitialCalibvsCrysDifferentRing",
                                                  "Sum initial cosmic calib const vs crystal ID, different theta rings;crystal ID", ECLElementNumbers::c_NCrystals, 0,
                                                  ECLElementNumbers::c_NCrystals);
  registerObject<TH1F>("InitialCalibvsCrysDifferentRing", InitialCalibvsCrysDifferentRing);
 
  auto CalibEntriesvsCrysSameRing = new TH1F("CalibEntriesvsCrysSameRing",
                                             "Entries in calib vs crys histograms, same theta ring;crystal ID;Entries per crystal", ECLElementNumbers::c_NCrystals, 0,
                                             ECLElementNumbers::c_NCrystals);
  registerObject<TH1F>("CalibEntriesvsCrysSameRing", CalibEntriesvsCrysSameRing);
 
  auto CalibEntriesvsCrysDifferentRing = new TH1F("CalibEntriesvsCrysDifferentRing",
                                                  "Entries in calib vs crys histograms, different theta rings;crystal ID;Entries per crystal", ECLElementNumbers::c_NCrystals, 0,
                                                  ECLElementNumbers::c_NCrystals);
  registerObject<TH1F>("CalibEntriesvsCrysDifferentRing", CalibEntriesvsCrysDifferentRing);
 
  auto RawDigitAmpvsCrys = new TH2F("RawDigitAmpvsCrys", "Digit Amplitude vs crystal ID;crystal ID;Amplitude",
                                    ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals, 100, 0,
                                    2000);
  registerObject<TH2F>("RawDigitAmpvsCrys", RawDigitAmpvsCrys);

  B2INFO("Input parameters to eclCosmicECollector:");
  B2INFO("trigThreshold: " << m_trigThreshold);
  B2INFO("minCrysE: " << m_minCrysE);

  FirstSet.resize(8737);
  EperCrys.resize(ECLElementNumbers::c_NCrystals);
  HitCrys.resize(ECLElementNumbers::c_NCrystals);
  EnergyPerTC.resize(576);
  TCperCrys.resize(ECLElementNumbers::c_NCrystals);
 

  if (m_ECLExpCosmicESame.hasChanged()) {ExpCosmicESame = m_ECLExpCosmicESame->getCalibVector();}
  if (m_ECLExpCosmicEDifferent.hasChanged()) {ExpCosmicEDifferent = m_ECLExpCosmicEDifferent->getCalibVector();}
  if (m_ElectronicsCalib.hasChanged()) {ElectronicsCalib = m_ElectronicsCalib->getCalibVector();}
  if (m_CosmicECalib.hasChanged()) {CosmicECalib = m_CosmicECalib->getCalibVector();}

  for (int ic = 1; ic < 9000; ic += 1000) {
    B2INFO("DB constants for cellID=" << ic << ": ExpCosmicESame = " << ExpCosmicESame[ic - 1] << " ExpCosmicEDifferent = " <<
           ExpCosmicEDifferent[ic - 1] << " ElectronicsCalib = " << ElectronicsCalib[ic - 1] << " CosmicECalib = " << CosmicECalib[ic - 1]);
  }

  for (int ic = 0; ic < ECLElementNumbers::c_NCrystals; ic++) {
    if (ElectronicsCalib[ic] <= 0) {B2FATAL("eclCosmicECollector: ElectronicsCalib = " << ElectronicsCalib[ic] << " for crysID = " << ic);}
    if (ExpCosmicESame[ic] == 0) {B2FATAL("eclCosmicECollector: ExpCosmicESame = 0 for crysID = " << ic);}
    if (ExpCosmicEDifferent[ic] == 0) {B2FATAL("eclCosmicECollector: ExpCosmicEDifferent = 0 for crysID = " << ic);}
    if (CosmicECalib[ic] == 0) {B2FATAL("eclCosmicECollector: CosmicECalib = 0 for crysID = " << ic);}
  }


  const short m_crystalsPerRing[69] = {
    48, 48, 64, 64, 64, 96, 96, 96, 96, 96, 96, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 96, 96, 96, 96, 96, 64, 64, 64
  };

  short firstCrystal[69] = {};
  for (int it = 1; it < 69; it++) {firstCrystal[it] = firstCrystal[it - 1] + m_crystalsPerRing[it - 1];}

  std::vector<short int> ThetaIDCrys(ECLElementNumbers::c_NCrystals);
  std::vector<short int> PhiIDCrys(ECLElementNumbers::c_NCrystals);
  int cID = 0;
  for (int it = 0; it < 69; it++) {
    for (int ic = 0; ic < m_crystalsPerRing[it]; ic++) {
      ThetaIDCrys.at(cID) = it;
      PhiIDCrys.at(cID) = ic;
      cID++;
    }
  }

  TrgEclMapping* trgecl_obj = new TrgEclMapping();
  for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
    TCperCrys[crysID] = trgecl_obj->getTCIdFromXtalId(crysID + 1) - 1;
  }

 
  int ThetaID = 0;
  double n2Overn0 = (double)m_crystalsPerRing[ThetaID + 2] / (double)m_crystalsPerRing[ThetaID];
  for (int phiID = 0; phiID < m_crystalsPerRing[ThetaID]; phiID++) {
    int crysID = firstCrystal[ThetaID] + phiID;

    int phiIDA = phiID + 1;
    if (phiIDA == m_crystalsPerRing[ThetaID]) {phiIDA = 0;}
    int phiIDB = phiID - 1;
    if (phiIDB == -1) {phiIDB = m_crystalsPerRing[ThetaID] - 1;}
    FirstSet[crysID] = CenterCrys.size();
    CenterCrys.push_back(crysID);
    NeighbourA.push_back(firstCrystal[ThetaID] + phiIDA);
    NeighbourB.push_back(firstCrystal[ThetaID] + phiIDB);

    double dphiB = n2Overn0 * phiID + 0.0001;
    phiIDB = dphiB;
    CenterCrys.push_back(crysID);
    NeighbourA.push_back(firstCrystal[ThetaID + 1] + phiID);
    NeighbourB.push_back(firstCrystal[ThetaID + 2] + phiIDB);
 
    phiIDB++;
    CenterCrys.push_back(crysID);
    NeighbourA.push_back(firstCrystal[ThetaID + 1] + phiID);
    NeighbourB.push_back(firstCrystal[ThetaID + 2] + phiIDB);
  }


  ECLNeighbours* myNeighbours4 = new ECLNeighbours("F", 0.95);
 
  for (int crysID = firstCrystal[1]; crysID < firstCrystal[68]; crysID++) {
    std::vector<short int> neighbours = myNeighbours4->getNeighbours(crysID + 1);

    int nA = -1;
    int nB = -1;
    std::vector<short int> nextThetaNeighbours;
    std::vector<short int> previousThetaNeighbours;
    for (auto& ID1 : neighbours) {
      int temp0 = ID1 - 1;
      if (temp0 != crysID && ThetaIDCrys[temp0] == ThetaIDCrys[crysID] && nA == -1) {
        nA = temp0;
      } else if (temp0 != crysID && ThetaIDCrys[temp0] == ThetaIDCrys[crysID] && nA != -1) {
        nB = temp0;
      }

      if (ThetaIDCrys[temp0] == ThetaIDCrys[crysID] + 1) {nextThetaNeighbours.push_back(temp0);}
      if (ThetaIDCrys[temp0] == ThetaIDCrys[crysID] - 1) {previousThetaNeighbours.push_back(temp0);}
    }

    if (nA >= 0 && nB >= 0) {
      FirstSet[crysID] = CenterCrys.size();
      CenterCrys.push_back(crysID);
      NeighbourA.push_back(nA);
      NeighbourB.push_back(nB);
    } else {
      B2FATAL("No neighbour pair with the same thetaID for crysID = " << crysID);
    }

    for (auto& IDn : nextThetaNeighbours) {
      for (auto& IDp : previousThetaNeighbours) {
        CenterCrys.push_back(crysID);
        NeighbourA.push_back(IDn);
        NeighbourB.push_back(IDp);
      }
    }
  }

 
  ThetaID = 68;
  for (int phiID = 0; phiID < m_crystalsPerRing[ThetaID]; phiID++) {
    int crysID = firstCrystal[ThetaID] + phiID;

    int phiIDA = phiID + 1;
    if (phiIDA == m_crystalsPerRing[ThetaID]) {phiIDA = 0;}
    int phiIDB = phiID - 1;
    if (phiIDB == -1) {phiIDB = m_crystalsPerRing[ThetaID] - 1;}
    FirstSet[crysID] = CenterCrys.size();
    CenterCrys.push_back(crysID);
    NeighbourA.push_back(firstCrystal[ThetaID] + phiIDA);
    NeighbourB.push_back(firstCrystal[ThetaID] + phiIDB);

    CenterCrys.push_back(crysID);
    NeighbourA.push_back(firstCrystal[ThetaID - 1] + phiID);
    NeighbourB.push_back(firstCrystal[ThetaID - 2] + phiID);
  }

  FirstSet[ECLElementNumbers::c_NCrystals] = CenterCrys.size();
}
 

void eclCosmicECollectorModule::collect()
{

  if (iEvent == 0) {
    for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
      getObjectPtr<TH1F>("ExpEvsCrysSameRing")->Fill(crysID + 0.001, ExpCosmicESame[crysID]);
      getObjectPtr<TH1F>("ElecCalibvsCrysSameRing")->Fill(crysID + 0.001, ElectronicsCalib[crysID]);
      getObjectPtr<TH1F>("InitialCalibvsCrysSameRing")->Fill(crysID + 0.001, CosmicECalib[crysID]);
      getObjectPtr<TH1F>("CalibEntriesvsCrysSameRing")->Fill(crysID + 0.001);
 
      getObjectPtr<TH1F>("ExpEvsCrysDifferentRing")->Fill(crysID + 0.001, ExpCosmicEDifferent[crysID]);
      getObjectPtr<TH1F>("ElecCalibvsCrysDifferentRing")->Fill(crysID + 0.001, ElectronicsCalib[crysID]);
      getObjectPtr<TH1F>("InitialCalibvsCrysDifferentRing")->Fill(crysID + 0.001, CosmicECalib[crysID]);
      getObjectPtr<TH1F>("CalibEntriesvsCrysDifferentRing")->Fill(crysID + 0.001);
 
    }
  }

  /* Metadata; check if required DB objects have changed */
  if (iEvent % 10000 == 0) {B2INFO("eclCosmicECollector: iEvent = " << iEvent);}
  iEvent++;
 
  if (m_ECLExpCosmicESame.hasChanged()) {B2FATAL("eclCosmicECollector: ECLExpCosmicESame has changed");}
  if (m_ECLExpCosmicEDifferent.hasChanged()) {B2FATAL("eclCosmicECollector: ECLExpCosmicEDifferent has changed");}
  if (m_ElectronicsCalib.hasChanged()) {B2FATAL("eclCosmicECollector: ElectronicsCalib has changed");}
  if (m_CosmicECalib.hasChanged()) {
    B2DEBUG(9, "eclCosmicECollector: new values for ECLCosmicECalib");
    CosmicECalib = m_CosmicECalib->getCalibVector();
  }

  memset(&EperCrys[0], 0, EperCrys.size()*sizeof EperCrys[0]);
  memset(&EnergyPerTC[0], 0, EnergyPerTC.size()*sizeof EnergyPerTC[0]);
 
 
  for (auto& eclDigit : m_eclDigitArray) {
    int crysID = eclDigit.getCellId() - 1;

    EperCrys[crysID] = eclDigit.getAmp() * abs(CosmicECalib[crysID]) * ElectronicsCalib[crysID];
    EnergyPerTC[TCperCrys[crysID]] += EperCrys[crysID];
    getObjectPtr<TH2F>("RawDigitAmpvsCrys")->Fill(crysID + 0.001, eclDigit.getAmp());
  }
  for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {HitCrys[crysID] = EperCrys[crysID] > m_minCrysE;}

  if (m_mockupL1) {
    float maxTrig = 0.;
    for (int tc = 0; tc < 576; tc++) {
      if (EnergyPerTC[tc] > maxTrig) {maxTrig = EnergyPerTC[tc];}
    }
    bool ECLtrigger = maxTrig > m_trigThreshold;
    if (!ECLtrigger) {return;}
  }

  for (int i = 0; i < FirstSet[ECLElementNumbers::c_NCrystals]; i++) {
    if (HitCrys[NeighbourA[i]] && HitCrys[NeighbourB[i]]) {
      int crysID = CenterCrys[i];
      if (i == FirstSet[crysID]) {

        getObjectPtr<TH2F>("EnvsCrysSameRing")->Fill(crysID + 0.001, EperCrys[crysID] / abs(ExpCosmicESame[crysID]));
        getObjectPtr<TH1F>("ExpEvsCrysSameRing")->Fill(crysID + 0.001, ExpCosmicESame[crysID]);
        getObjectPtr<TH1F>("ElecCalibvsCrysSameRing")->Fill(crysID + 0.001, ElectronicsCalib[crysID]);
        getObjectPtr<TH1F>("InitialCalibvsCrysSameRing")->Fill(crysID + 0.001, CosmicECalib[crysID]);
        getObjectPtr<TH1F>("CalibEntriesvsCrysSameRing")->Fill(crysID + 0.001);
      } else {
        getObjectPtr<TH2F>("EnvsCrysDifferentRing")->Fill(crysID + 0.001, EperCrys[crysID] / abs(ExpCosmicEDifferent[crysID]));
        getObjectPtr<TH1F>("ExpEvsCrysDifferentRing")->Fill(crysID + 0.001, ExpCosmicEDifferent[crysID]);
        getObjectPtr<TH1F>("ElecCalibvsCrysDifferentRing")->Fill(crysID + 0.001, ElectronicsCalib[crysID]);
        getObjectPtr<TH1F>("InitialCalibvsCrysDifferentRing")->Fill(crysID + 0.001, CosmicECalib[crysID]);
        getObjectPtr<TH1F>("CalibEntriesvsCrysDifferentRing")->Fill(crysID + 0.001);
      }
    }
  }
}