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"
#include "tkdcards.h"
extern int StopSelector;

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

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

  m_logger->info("Creating output file {}", m_outputFilename);
  std::string output_file_option = m_reprocess_rti ? "update" : "recreate";
  m_outputFile = std::make_unique<TFile>(m_outputFilename.c_str(), output_file_option.c_str());
  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);
  { // propagate MC flag to containers
    m_chain->m_event.evSummary.SetMC(m_isMC);
    m_chain->m_event.richPlus.SetMC(m_isMC);
    m_chain->m_event.tofBase.SetMC(m_isMC);
    m_chain->m_event.extHitBase.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 (StopSelector)
    Abort("Aborting TSelector Loop\n");
  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) {
    SetDefaultMCTuningParameters();
    TRCLFFKEY.ClusterCofGOpt = 1;
    TRFITFFKEY.Zshift = -1;

    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);
  } else {
    TrLinearEtaDB::SetLinearCluster();
    TRFITFFKEY.Zshift = 2;
  }

  TRFITFFKEY.ErcHeY = 0;

  // quickly discard entries we don't care about...
  if (Entry() < firstentry) {
    m_logger->trace("Skip entry {} < {}", Entry(), firstentry);
    return false;
  } else if (Entry() > lastentry) {
    // FIXME: check is Paolo's patch to gbatch is accepted upstream. Then set StopSelector
    // StopSelector = true;
    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");

  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();
  }

  // NOTE: We stop here if we only have to process RTIInfo :)
  if (m_reprocess_rti)
    return true;

  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) {

    // TofBase
    {
      auto stpw = m_benchmark.MakeStopwatch("TofBase::Fill");
      bool has_tof = m_chain->m_event.tofBase.Fill(particle->pBeta(), particle->pBetaH());
      if (has_tof) {
        m_chain->m_event.header.SetMaskCategory(NAIA::Category::HasTof);
      }
      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());

      // 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]);
        }

        bool has_track = m_chain->m_event.trTrackBase.Fill(particle->pTrTrack());
        if (has_track) {
          m_chain->m_event.header.SetMaskCategory(NAIA::Category::HasTrack);
        }

        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());
      }
      // SecondTrTrackBase
      if (NTrTrack() > 1) {
        auto stpw = m_benchmark.MakeStopwatch("SecondTrTrackBase::Fill");
        TrTrackR *second_track = nullptr;
        for (TrTrackR &track : TrTrack()) {
          if (&track == particle->pTrTrack() || track.GetNhits() < 5) {
            continue;
          }

          if (!second_track || second_track->GetRigidity() < track.GetRigidity()) {
            second_track = &track;
          }
        }

        if (second_track) {
          m_chain->m_event.secondTrTrackBase.Fill(second_track, TrTrackFillType::SecondTrack);
        }
      }
      // Standalone TrTrack
      {
        auto stpw = m_benchmark.MakeStopwatch("TrTrackBaseStandalone::Fill");
        m_chain->m_event.trTrackBaseSt.Fill(particle->pTrTrack(), TrTrackFillType::Standalone);

        float innerCharge_St = m_chain->m_event.trTrackBaseSt.InnerCharge[NAIA::TrTrack::ChargeRecoType::STD];
        if (std::fabs(innerCharge_St - 1) < 0.5) {
          m_chain->m_event.SetMaskCategory(NAIA::Category::Charge1_Trk_St);
        } else if (std::fabs(innerCharge_St - 2) < 0.5) {
          m_chain->m_event.SetMaskCategory(NAIA::Category::Charge2_Trk_St);
        } else if (innerCharge_St > 2.5) {
          m_chain->m_event.SetMaskCategory(NAIA::Category::ChargeGT2_Trk_St);
        }
      }

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

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

        if (trdK) {
          // set beta for TRD filling
          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.trdKBase.SetBeta(m_chain->m_event.tofBase.Beta[NAIA::Tof::BetaType::BetaH]);
            m_chain->m_event.trdKBase.Type = TrdK::BuildType::TrkTrack;
          }

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

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

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

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

    // fill standalone containers
    {
      auto stpw = m_benchmark.MakeStopwatch("DefineStandalone");
      DefineStandAlone();
    }

    if (m_betahStPtr) {
      // TofBaseStandalone
      {
        auto stpw = m_benchmark.MakeStopwatch("TofBaseStandalone::Fill");
        bool has_tof_standalone = m_chain->m_event.tofBaseSt.Fill(nullptr, m_betahStPtr, true);
        if (has_tof_standalone) {
          m_chain->m_event.header.SetMaskCategory(NAIA::Category::HasTofStandalone);
        }

        float UTOFCharge_St = m_chain->m_event.tofBaseSt.ChargeNoPL[NAIA::Tof::ChargeType::Upper];
        if (std::fabs(UTOFCharge_St - 1) < 0.5) {
          m_chain->m_event.SetMaskCategory(NAIA::Category::Charge1_Tof_St);
        } else if (std::fabs(UTOFCharge_St - 2) < 0.5) {
          m_chain->m_event.SetMaskCategory(NAIA::Category::Charge2_Tof_St);
        } else if (UTOFCharge_St > 2.5) {
          m_chain->m_event.SetMaskCategory(NAIA::Category::ChargeGT2_Tof_St);
        }
      }
      // TofPlusStandalone
      {
        auto stpw = m_benchmark.MakeStopwatch("TofPlusStandalone::Fill");
        m_chain->m_event.tofPlusSt.Fill(nullptr, m_betahStPtr, true);
      }

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

        Logging::MuteOutput();
        auto trdK = std::make_unique<TrdKCluster>();
        if (m_trdTrackStPtr) {
          m_logger->trace("Building standalone TRD from TrdTrack");
          trdK->Build(m_trdTrackStPtr);
          m_chain->m_event.trdKBaseSt.Type = TrdK::BuildType::TrdTrack;
        } else {
          m_logger->trace("Building standalone TRD from BetaH");
          trdK->Build(m_betahStPtr);
          m_chain->m_event.trdKBaseSt.Type = TrdK::BuildType::BetaH;
        }
        Logging::UnmuteOutput();

        if (trdK) {
          // set beta for TRD filling
          if (MatchAllBits(m_chain->m_event.header.Mask(), NAIA::Category::HasTofStandalone) &&
              ContainsKeys(m_chain->m_event.tofBaseSt.Beta, NAIA::Tof::BetaType::BetaH)) {
            m_chain->m_event.trdKBaseSt.SetBeta(m_chain->m_event.tofBaseSt.Beta[NAIA::Tof::BetaType::BetaH]);
          }

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

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

      bool use_trd = MatchAllBits(m_chain->m_event.header.Mask(), NAIA::Category::HasTrdStandalone);

      Logging::MuteOutput();
      m_chain->m_event.extHitBase.Fill(this, use_trd);
      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");
        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.trTrackBase.FillElectronVars(particle->pTrTrack(), zEcalCOG, refitKalman);
      }
      {
        auto stpw = m_benchmark.MakeStopwatch("TrTrackPlus::FillElectronVars");
        m_chain->m_event.trTrackPlus.FillElectronVars(particle->pTrTrack(), 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();

  if (m_reprocess_rti) {
    m_logger->info("Overwriting RTIInfo tree");
    m_chain->m_rti.tree->BuildIndex("RTIInfo.UTCTime");
    m_outputFile->WriteTObject(m_chain->m_rti.tree, "RTI", "Overwrite");
  } else {
    m_logger->info("Writing chain to file");
    NAIA::Version versionHeader;
    m_outputFile->WriteTObject(&versionHeader, "VersionHeader");
    m_chain->Write();
  }
  m_benchmark.Print(m_logger);
}

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

  // 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());

  auto getMaxQBetaH = [this](bool matchTrd = true) {
    int ib = -1;
    BetaHR *maxQBetaH = nullptr;

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

    double prev_tofcharge = 0;

    for (BetaHR &betah : BetaH()) {
      ++ib;
      float temp_tofcharge = betah.GetQ(nlayers, qrms);
      if (temp_tofcharge <= 0)
        continue;

      if (matchTrd) {
        if (!betah.pTrdTrack())
          continue;
        // count TRD track hits
        // NOTE(vformato): Is this the right thing?
        unsigned int nTrdHits{0};
        for (unsigned int i = 0; i < betah.pTrdTrack()->NTrdSegment(); ++i) {
          nTrdHits += betah.pTrdTrack()->pTrdSegment(i)->NTrdCluster();
        }

        if (nTrdHits < 15)
          continue;

        m_logger->trace("StAl BetaH n. {} - Charge {} (Chi2C: {} - Chi2T: {}) - TRD hits: {}", ib, temp_tofcharge,
                        betah.GetChi2C(), betah.GetChi2T(), nTrdHits);
      } else {
        m_logger->trace("StAl BetaH n. {} - Charge {} (Chi2C: {} - Chi2T: {})", ib, temp_tofcharge, betah.GetChi2C(),
                        betah.GetChi2T());
      }

      if (temp_tofcharge > prev_tofcharge) {
        maxQBetaH = &betah;
        prev_tofcharge = temp_tofcharge;
      }
    }

    return maxQBetaH;
  };

  m_betahStPtr = getMaxQBetaH(false);
  m_trdTrackStPtr = m_betahStPtr ? m_betahStPtr->pTrdTrack() : nullptr;

  // just for safety, restore old opt
  TofRecH::BuildOpt = oldOpt;
}

} // namespace NAIA