NtpSelector.cpp

Go to the documentation of this file.

#include "Utility/Logging/Logging.h"

#include "NtpMaker/NtpSelector.h"
#include "NtpMaker/NtpTools.h"

#include "Tofrec02_ihep.h"

namespace NAIA {
void NtpSelector::UBegin() {
  m_logger->info("UBegin =============================================");

  auto stpw = m_benchmark.MakeStopwatch("NtpSelector::UBegin");

  m_logger->info("Creating output file {}", m_outputFilename);
  m_outputFile = std::make_unique<TFile>(m_outputFilename.c_str(), "recreate");
  m_outputFile->SetCompressionLevel(7);

  m_logger->info("Creating NAIA event tree");
  m_chain = std::make_unique<NAIAChain>(NAIAChain::AccessMode::Write);
  m_chain->m_event.SetMC(m_isMC);
  m_chain->SetupBranches(m_isMC);
  m_logger->info("NAIA branches created");
}

bool NtpSelector::Notify() {
  static bool firstCall = true;
  if (!firstCall) {
    m_logger->debug("Switching file in chain: {}", m_tree->GetCurrentFile()->GetName());
    // reset file info containers
    m_chain->m_fileInfo.Clear();
    if (m_isMC)
      m_chain->m_fileInfoMC.Clear();
  } else {
    firstCall = false;
  }
  return AMSEventR::Notify();
}

void NtpSelector::ProcessFileInfo() {
  m_logger->debug("Called NtpSelector::ProcessFileInfo");

  m_chain->m_fileInfo.Clear();
  if (m_isMC)
    m_chain->m_fileInfoMC.Clear();

  m_chain->m_fileInfo.FileName = m_tree->GetCurrentFile()->GetName();

  m_chain->m_fileInfo.Run = Run();
  m_chain->m_fileInfo.EventNo.first = Event();
  m_chain->m_fileInfo.EventNo.second = Event();

  m_chain->m_fileInfo.UTCTime.first = fHeader.Time[0];
  m_chain->m_fileInfo.UTCTime.second = fHeader.Time[0];

  if (m_isMC) {
    m_chain->m_fileInfoMC.EventNo.first = Event();
    m_chain->m_fileInfoMC.EventNo.second = Event();
    MCEventgR *primary = (nMCEventg() != 0) ? pMCEventg(0) : nullptr;
    m_chain->m_fileInfoMC.Charge = (primary) ? static_cast<int>(primary->Charge) : 0;
    m_chain->m_fileInfoMC.Mass = (primary) ? primary->Mass : 0;
    m_chain->m_fileInfoMC.MomentumRange.first = std::numeric_limits<float>::max();
    m_chain->m_fileInfoMC.MomentumRange.second = 0;
    auto *datacard = m_tree->GetCurrentFile()->Get<TObjString>("DataCards");
    if (!datacard) {
      m_logger->error("No Datacard Object available!");
    }
    m_chain->m_fileInfoMC.DatacardPID = NtpTools::GetDatacardValue<int>(datacard, "PART");
    m_chain->m_fileInfoMC.DatacardMomentumRange.first = NtpTools::GetDatacardValue<float>(datacard, "PMIN");
    m_chain->m_fileInfoMC.DatacardMomentumRange.second = NtpTools::GetDatacardValue<float>(datacard, "PMAX");
    m_chain->m_fileInfoMC.DatacardNGen = NtpTools::GetDatacardValue<unsigned int>(datacard, "TRIG");

    if (m_chain->m_fileInfo.FileName.find(".l1.") != std::string::npos) {
      m_chain->m_fileInfoMC.Focus = MCFileInfo::SimFocus::L1;
    } else if (m_chain->m_fileInfo.FileName.find(".l19.") != std::string::npos) {
      m_chain->m_fileInfoMC.Focus = MCFileInfo::SimFocus::L19;
    } else if (m_chain->m_fileInfo.FileName.find(".l1a9.") != std::string::npos) {
      m_chain->m_fileInfoMC.Focus = MCFileInfo::SimFocus::L19;
    } else if (m_chain->m_fileInfo.FileName.find(".tb.") != std::string::npos) {
      m_chain->m_fileInfoMC.Focus = MCFileInfo::SimFocus::TB;
    }
  } else {
    // not MC
    // [VF] here we reproduce the logic of AMSEventR::IsBadRun. Can't use that since the "reason" parameter is thrown
    // away and we actually need it
    if (getsetup())
      m_chain->m_fileInfo.BadRunTag = getsetup()->IsBadRun(m_chain->m_fileInfo.BadRunReason, UTime(), Run());
    else
      m_chain->m_fileInfo.BadRunTag = 2;
  }
}

bool NtpSelector::UProcessCut() {
  if (m_chain->m_fileInfo.Run == 0) {
    // this happens when a file is first loaded
    ProcessFileInfo();
  }

  // fill FileInfo stuff...
  m_chain->m_fileInfo.NEvents++;
  m_chain->m_fileInfo.EventNo.second = Event();
  m_chain->m_fileInfo.UTCTime.second = fHeader.Time[0];

  if (m_isMC) {
    m_chain->m_fileInfoMC.EventNo.first = std::min(Event(), m_chain->m_fileInfoMC.EventNo.first);
    m_chain->m_fileInfoMC.EventNo.second = std::max(Event(), m_chain->m_fileInfoMC.EventNo.second);

    MCEventgR *primary = (nMCEventg() != 0) ? pMCEventg(0) : nullptr;
    float p = (primary) ? primary->Momentum : 0;
    m_chain->m_fileInfoMC.MomentumRange.first = std::min(p, m_chain->m_fileInfoMC.MomentumRange.first);
    m_chain->m_fileInfoMC.MomentumRange.second = std::max(p, m_chain->m_fileInfoMC.MomentumRange.second);
  }

  // quickly discard entries we don't care about...
  if (Entry() < firstentry) {
    m_logger->trace("Skip entry {} < {}", Entry(), firstentry);
    return false;
  } else if (Entry() > lastentry) {
    // TODO: check why Abort doesn't work...
    m_logger->trace("Skip entry {} > {}", Entry(), lastentry);
    return false;
  }

  static bool got_it = false;
  if (!got_it) {
    m_logger->debug("UProcessCut - First event processed - Entry: {},  Run: {}, Event: {}", Entry(), Run(), Event());
  } else if (Entry() % 10000 == 0) {
    m_logger->debug("UProcessCut - Event processed - Entry: {},  Run: {}, Event: {}", Entry(), Run(), Event());
  }
  got_it = true;

  m_logger->trace("Entry: {}", Entry());

  m_chain->Clear();

  // cut really bad stuff
  if ((nDaqEvent() < 1) || (nLevel1() < 1)) {
    return false;
  }

  auto stpw = m_benchmark.MakeStopwatch("Event processing");

  // temporary event processing for prototyping
  static unsigned long lastSec = 0;
  //  int ibetah_select = -1;
  ParticleR *particle = nullptr;

  if (m_isMC) {
    auto stpw = m_benchmark.MakeStopwatch("MCTruth::Fill");
    m_chain->m_event.mcTruthBase.Fill(this);
    m_chain->m_event.mcTruthPlus.Fill(this);
  }

  {
    auto stpw = m_benchmark.MakeStopwatch("(Header+EventSummary)::Fill");
    m_chain->m_event.header.Fill(this);
    m_chain->m_event.evSummary.Fill(this);
  }

  // RTI
  if (!m_isMC && lastSec != m_chain->m_event.header.UTCTime) {
    auto stpw = m_benchmark.MakeStopwatch("RTIInfo::Fill");

    // don't fill the first time we enter here
    if (lastSec > 0) {
      m_chain->m_rtiInfo.FillEventVariables();
      m_chain->FillRTI();
    }
    m_chain->m_rtiInfo.Clear();

    Logging::MuteOutput();
    GetRTI(m_rti);
    m_chain->m_rtiInfo.Fill(&m_rti, m_chain->m_event.header.RunTag);
    Logging::UnmuteOutput();
    lastSec = m_chain->m_event.header.UTCTime;
  }

  if (!m_isMC) {
    auto stpw = m_benchmark.MakeStopwatch("RTIInfo::AccumulateEventVariables");
    Logging::MuteOutput();
    m_chain->m_rtiInfo.AccumulateEventVariables(this);
    Logging::UnmuteOutput();
  }

  if (pDaqEvent(0)) {
    auto stpw = m_benchmark.MakeStopwatch("DAQ::Fill");
    m_chain->m_event.daq.Fill(*pDaqEvent(0));
  }

  // temporary! To be reviewed
  particle = pParticle(0);

  if (particle) {

    if (particle->pBetaH()) {
      m_chain->m_event.header.SetMaskCategory(NAIA::Category::HasTof);
    }
    // TofBase
    {
      auto stpw = m_benchmark.MakeStopwatch("TofBase::Fill");
      m_chain->m_event.tofBase.Fill(particle->pBeta(), particle->pBetaH());
      float UTOFCharge = m_chain->m_event.tofBase.Charge[NAIA::Tof::ChargeType::Upper];
      if (std::fabs(UTOFCharge - 1) < 0.5) {
        m_chain->m_event.SetMaskCategory(NAIA::Category::Charge1_Tof);
      } else if (std::fabs(UTOFCharge - 2) < 0.5) {
        m_chain->m_event.SetMaskCategory(NAIA::Category::Charge2_Tof);
      } else if (UTOFCharge > 2.5) {
        m_chain->m_event.SetMaskCategory(NAIA::Category::ChargeGT2_Tof);
      }
    }
    // TofPlus
    {
      auto stpw = m_benchmark.MakeStopwatch("TofPlus::Fill");
      m_chain->m_event.tofPlus.Fill(particle->pBeta(), particle->pBetaH());
    }

    if (particle->pTrTrack()) {
      if (m_isMC)
        m_chain->m_event.mcTruthPlus.FillTrMCHit(particle->pTrTrack());

      m_chain->m_event.header.SetMaskCategory(NAIA::Category::HasTrack);

      // TrTrackBase
      {
        auto stpw = m_benchmark.MakeStopwatch("TrTrackBase::Fill");

        m_logger->trace("Track no. = {}", particle->iTrTrack());

        if (MatchAllBits(m_chain->m_event.header.Mask(), NAIA::Category::HasTof) &&
            ContainsKeys(m_chain->m_event.tofBase.Beta, NAIA::Tof::BetaType::BetaH)) {
          m_chain->m_event.trTrackBase.SetBeta(m_chain->m_event.tofBase.Beta[NAIA::Tof::BetaType::BetaH]);
        }

        m_chain->m_event.trTrackBase.Fill(particle->pTrTrack());

        float innerCharge = m_chain->m_event.trTrackBase.InnerCharge[NAIA::TrTrack::ChargeRecoType::STD];
        if (std::fabs(innerCharge - 1) < 0.5) {
          m_chain->m_event.SetMaskCategory(NAIA::Category::Charge1_Trk);
        } else if (std::fabs(innerCharge - 2) < 0.5) {
          m_chain->m_event.SetMaskCategory(NAIA::Category::Charge2_Trk);
        } else if (innerCharge > 2.5) {
          m_chain->m_event.SetMaskCategory(NAIA::Category::ChargeGT2_Trk);
        }
      }
      // TrTrackPlus
      {
        auto stpw = m_benchmark.MakeStopwatch("TrTrackPlus::Fill");
        m_chain->m_event.trTrackPlus.Fill(particle->pTrTrack());
      }
      // TrdKBase
      {
        auto stpw = m_benchmark.MakeStopwatch("TrdKBase::Fill");

        Logging::MuteOutput();
        auto trdK = std::make_unique<TrdKCluster>(this, particle->pTrTrack(), 0);
        Logging::UnmuteOutput();

        if (trdK) {
          m_chain->m_event.header.SetMaskCategory(NAIA::Category::HasTrd);
          m_chain->m_event.trdKBase.Fill(trdK.get(), particle->pTrTrack());
        }
      }

      // we still need to fill electron-related track fits
    }

    if (particle->pEcalShower()) {
      m_chain->m_event.header.SetMaskCategory(NAIA::Category::HasEcal);
      // EcalBase
      {
        auto stpw = m_benchmark.MakeStopwatch("EcalBase::Fill");
        m_chain->m_event.ecalBase.Fill(particle->pEcalShower());
      }
      // EcalPlus
      {
        auto stpw = m_benchmark.MakeStopwatch("EcalPlus::Fill");
        m_chain->m_event.ecalPlus.Fill(particle->pEcalShower());
      }
    }

    // fill standalone containers
    {
      auto stpw = m_benchmark.MakeStopwatch("DefineStandalone");
      DefineStandAlone();
    }
    if (m_betahStPtr) {
      m_chain->m_event.header.SetMaskCategory(NAIA::Category::HasTofStandalone);
      // TofBaseStandalone
      {
        auto stpw = m_benchmark.MakeStopwatch("TofBaseStandalone::Fill");
        m_chain->m_event.tofBaseSt.Fill(nullptr, m_betahStPtr, true);
      }
      // TofPlusStandalone
      {
        auto stpw = m_benchmark.MakeStopwatch("TofPlusStandalone::Fill");
        m_chain->m_event.tofPlusSt.Fill(nullptr, m_betahStPtr, true);
      }

      // TrdKBase
      {
        auto stpw = m_benchmark.MakeStopwatch("TrdKBaseStandalone::Fill");

        Logging::MuteOutput();
        auto trdK = std::make_unique<TrdKCluster>();
        if (m_betahStPtr->pTrdTrack()) {
          trdK->Build(m_betahStPtr->pTrdTrack());
        } else {
          trdK->Build(m_betahStPtr);
        }
        Logging::UnmuteOutput();

        if (trdK) {
          m_chain->m_event.header.SetMaskCategory(NAIA::Category::HasTrdStandalone);
          m_chain->m_event.trdKBaseSt.Fill(trdK.get(), nullptr);
        }
      }
    }

    if (particle->pRichRing()) {
      m_chain->m_event.header.SetMaskCategory(NAIA::Category::HasRich);
      // RichBase
      {
        auto stpw = m_benchmark.MakeStopwatch("RichBase::Fill");
        m_chain->m_event.richBase.Fill(particle->pRichRing(), particle->pRichRingB());
      }
      // RichPlus
      {
        auto stpw = m_benchmark.MakeStopwatch("RichPlus::Fill");
        m_chain->m_event.richPlus.Fill(this, particle->pRichRing(), particle->pTrTrack());
      }
    }

    // unbiased ExtHit
    if (MatchAllBits(m_chain->m_event.header.Mask(), NAIA::Category::HasTofStandalone)) {
      auto stpw = m_benchmark.MakeStopwatch("UnbExtHitBase::Fill");
      Logging::MuteOutput();
      m_chain->m_event.extHitBase.Fill(this);
      Logging::UnmuteOutput();
    }

    // only at the very end we screw up gbatch's tracker fits internals and fill electron-related variables
    if (particle->pTrTrack() && MatchAllBits(m_chain->m_event.header.Mask(), Category::HasTrd)) {
      // condition for Kalman electron fit: ECAL EnergyD > 0.4 GeV
      bool refitKalman = MatchAllBits(m_chain->m_event.header.Mask(), Category::HasEcal) &&
                         m_chain->m_event.ecalBase.Energy[Ecal::EnergyRecoType::EnergyD] > 0.4;
      {
        auto stpw = m_benchmark.MakeStopwatch("TrTrackBase::FillElectronVars");
        m_chain->m_event.trTrackBase.FillElectronVars(particle->pTrTrack(), refitKalman);
      }
      {
        auto stpw = m_benchmark.MakeStopwatch("TrTrackPlus::FillElectronVars");
        float zEcalCOG = MatchAllBits(m_chain->m_event.header.Mask(), Category::HasEcal)
                             ? static_cast<float>(m_chain->m_event.ecalBase.GetCOG().z())
                             : std::numeric_limits<float>::max();
        m_chain->m_event.trTrackPlus.FillElectronVars(particle->pTrTrack(), zEcalCOG, refitKalman);
      }
    }

    m_logger->trace("Event mask: {}", to_string_binary<32>(m_chain->m_event.header.Mask()));
  } else {
    // do selection without particle here?
  }

  m_chain->Fill();

  return true;
}

void NtpSelector::UProcessFill() {}

void NtpSelector::UTerminate() {
  // last call for accumulated info
  m_chain->FillFileInfo();

  m_chain->m_rtiInfo.FillEventVariables();
  m_chain->FillRTI();

  m_logger->info("UTerminate =============================================");
  // m_chain->Print();

  m_outputFile->cd();
  NAIA::Version versionHeader;
  m_outputFile->WriteTObject(&versionHeader, "VersionHeader");
  m_chain->Write();

  m_benchmark.Print(m_logger);
}

void NtpSelector::DefineStandAlone() {
  m_betahStPtr = nullptr;

  // dummy variables for gbatch terrible interface
  int nlayers;
  float qrms;

  unsigned int ib = 0;
  double prev_tofcharge = 0;
  // Rebuild Tof with standalone reconstruction
  auto oldOpt = TofRecH::BuildOpt;
  TofRecH::BuildOpt = 3110;
  m_logger->trace("Rebuilding Tof with opt {}", TofRecH::BuildOpt);
  TofRecH::ReBuild();
  m_logger->trace("{} BetaH objects found", nBetaH());

  // Search Max Q standalone BetaH
  for (BetaHR &betah : BetaH()) {
    auto tofGoodPatt = 1000 * betah.IsGoodQPathL(0) + 100 * betah.IsGoodQPathL(1) + 10 * betah.IsGoodQPathL(2) +
                       1 * betah.IsGoodQPathL(3);
    float temp_tofcharge = betah.GetQ(nlayers, qrms, 2, TofClusterHR::DefaultQOptIon, tofGoodPatt);
    if (temp_tofcharge <= 0)
      continue;

    m_logger->trace("StAl BetaH n. {} - Charge {} (Chi2C: {} - Chi2T: {})", ib, temp_tofcharge, betah.GetChi2C(),
                    betah.GetChi2T());
    if (temp_tofcharge > prev_tofcharge) {
      m_betahStPtr = &betah;
      prev_tofcharge = temp_tofcharge;
    }
    ib++;
  }

  // search for same charge, better spatial chisq.
  for (BetaHR &betah : BetaH()) {
    auto tofGoodPatt_h = 1000 * m_betahStPtr->IsGoodQPathL(0) + 100 * m_betahStPtr->IsGoodQPathL(1) +
                         10 * m_betahStPtr->IsGoodQPathL(2) + 1 * m_betahStPtr->IsGoodQPathL(3);
    auto tofGoodPatt = 1000 * betah.IsGoodQPathL(0) + 100 * betah.IsGoodQPathL(1) + 10 * betah.IsGoodQPathL(2) +
                       1 * betah.IsGoodQPathL(3);

    int high_charge = floor(m_betahStPtr->GetQ(nlayers, qrms, 2, TofClusterHR::DefaultQOptIon, tofGoodPatt_h) + 0.5);
    int temp_charge = floor(betah.GetQ(nlayers, qrms, 2, TofClusterHR::DefaultQOptIon, tofGoodPatt) + 0.5);

    if ((high_charge == temp_charge) && (betah.GetChi2C() < m_betahStPtr->GetChi2C())) {
      m_betahStPtr = &betah;
    }
  }

  if (m_betahStPtr) {
    m_logger->trace("Selected StAl BetaH: - Charge {} (Chi2C: {})",
                    m_betahStPtr->GetQ(nlayers, qrms, 2, TofClusterHR::DefaultQOptIon), m_betahStPtr->GetChi2C());
  }

  TofRecH::BuildOpt = oldOpt;
  m_logger->trace("Rebuilding Tof with opt {}", TofRecH::BuildOpt);
  TofRecH::ReBuild();
}

} // namespace NAIA