RTIInfoFill.cpp

Go to the documentation of this file.

#include "Containers/RTIInfo.h"

#include "TMath.h"

// gbatch headers
#include "GM_SubLibrary.h"
#include "root.h"

#include <algorithm>

namespace NAIA {

static constexpr double Re = 6372.795477598e5; // cm

template <typename ValueType, long unsigned int N>
constexpr int FindBin(ValueType value, const std::array<ValueType, N> &bins) {
  auto binIt = std::lower_bound(begin(bins), end(bins), value, std::less_equal<ValueType>());
  return binIt == end(bins) ? -1 : std::distance(begin(bins), binIt) - 1;
};

template <typename ValueType>
int FindBin(ValueType value, const std::map<unsigned int, std::pair<ValueType, ValueType>> &bins) {
  auto binIt = std::find_if(begin(bins), end(bins),
                            [&value](auto it) { return (value >= it.second.first && value < it.second.second); });
  return binIt == end(bins) ? -1 : binIt->first;
};

auto GoodParticle = [](ParticleR *part) {
  BetaHR *beta = part->pBetaH();

  TrTrackR *trk = part->pTrTrack();
  if (!trk)
    return false;

  constexpr int refit = 31;
  constexpr int fit_algo = 1;
  int id_max = trk->iTrTrackPar(fit_algo, 0, refit);
  if (id_max < 0)
    return false;

  bool hasinnery = false;
  bool L2y = false, L34y = false, L56y = false, L78y = false;
  TrRecHitR *trkHit = nullptr;
  trkHit = trk->GetHitLJ(2);
  if (trkHit)
    if (trkHit->pTrCluster('y')) {
      L2y = true;
    }
  trkHit = trk->GetHitLJ(3);
  if (trkHit)
    if (trkHit->pTrCluster('y')) {
      L34y = true;
    }
  trkHit = trk->GetHitLJ(4);
  if (trkHit)
    if (trkHit->pTrCluster('y')) {
      L34y = true;
    }
  trkHit = trk->GetHitLJ(5);
  if (trkHit)
    if (trkHit->pTrCluster('y')) {
      L56y = true;
    }
  trkHit = trk->GetHitLJ(6);
  if (trkHit)
    if (trkHit->pTrCluster('y')) {
      L56y = true;
    }
  trkHit = trk->GetHitLJ(7);
  if (trkHit)
    if (trkHit->pTrCluster('y')) {
      L78y = true;
    }
  trkHit = trk->GetHitLJ(8);
  if (trkHit)
    if (trkHit->pTrCluster('y')) {
      L78y = true;
    }
  if (L2y && L34y && L56y && L78y)
    hasinnery = true;

  float chi2 = trk->GetNormChisqY(id_max);
  if (!beta || beta->GetSumHit() < 4 || !trk || id_max < 0 || !hasinnery || chi2 < 0 || chi2 > 10)
    return false;

  return true;
};

void RTIInfo::Fill(AMSSetupR::RTI *rtiPtr, unsigned int runTag) {
  Run = rtiPtr->run;
  RunTag = runTag;
  FirstEvNo = rtiPtr->evno;
  LastEvNo = rtiPtr->evnol;
  LivetimeFraction = rtiPtr->lf;
  MostProbableHeRig = rtiPtr->mphe;
  Theta = rtiPtr->theta;
  Phi = rtiPtr->phi;
  Altitude = rtiPtr->r - Re;
  Zenith = rtiPtr->zenith;
  GalacticLat = rtiPtr->glat;
  GalacticLong = rtiPtr->glong;
  nEvent = rtiPtr->nev;
  nError = rtiPtr->nerr;
  nTrigger = rtiPtr->ntrig;
  nHWError = rtiPtr->nhwerr;
  nEventParticle = rtiPtr->npart;
  nTRDHitAvg = rtiPtr->mtrdh;
  nTrClusterAvg = rtiPtr->mtrcl;
  good = rtiPtr->good;
  UTCTime = rtiPtr->utime;

  for (unsigned int i = 0; i < 4; i++) {
    for (unsigned int j = 0; j < 2; j++) {
      MaxStoermerCutoff[i][j] = rtiPtr->cf[i][j];
      MaxIGRFCutoff[i][j] = rtiPtr->cfi[i][j];
      if (i < 2) {
        MinIGRFCutoff[i][j] = rtiPtr->cfl[i][j];
      }
    }
  }

  AMSPoint pn1, pn9, pd1, pd9;
  AMSEventR::GetRTIdL1L9(0, pn1, pd1, UTCTime, 60);
  AMSEventR::GetRTIdL1L9(1, pn9, pd9, UTCTime, 60);
  MeanAlignDiffExtLayer[0] = {static_cast<float>(pd1.x()), static_cast<float>(pd1.y()), static_cast<float>(pd1.z())};
  MeanAlignDiffExtLayer[1] = {static_cast<float>(pd9.x()), static_cast<float>(pd9.y()), static_cast<float>(pd9.z())};

  // again, this should be illegal
  double GMC[3];
  GM_GeoMagSphericalCoo(UTCTime, Altitude, Theta, Phi, GMC);

  MagThetaSphCoo = static_cast<float>(GMC[2]);
  MagPhiSphCoo = static_cast<float>(GMC[1]);
  MagAltSphCoo = static_cast<float>(GMC[0]);

  double lambda2 = cos(MagThetaSphCoo * TMath::DegToRad()) * cos(MagThetaSphCoo * TMath::DegToRad());
  LShell = (MagAltSphCoo / Re) / lambda2;

  AMSSetupR amssetup;
  AMSSetupR::DSPError dsperr;
  DSPError = amssetup.getDSPError(dsperr, UTCTime);
  DSPErrorJINJ = false;
  for (int inode = 0x80; inode < 0x88; inode++)
    if (dsperr.SearchNA(inode))
      DSPErrorJINJ = true;
  DSPErrorJLV1 = false;
  for (int inode = 0x88; inode < 0x8C; inode++)
    if (dsperr.SearchNA(inode))
      DSPErrorJLV1 = true;
  DSPErrorJINF = false;
  for (int inode = 0x34; inode < 0x3E; inode++)
    if (dsperr.SearchNA(inode))
      DSPErrorJINF = true;
  for (int inode = 0x96; inode < 0xCE; inode++)
    if (dsperr.SearchNA(inode))
      DSPErrorJINF = true;
}

void RTIInfo::FillEventVariables() {
  if (m_processedEvents == 0) {
    return;
  }

  auto normalize = [this](auto &value) { value /= m_processedEvents; };

  normalize(nTRDTrack);
  normalize(nTRDHits);

  std::for_each(begin(nTRDLayerHits), begin(nTRDLayerHits), normalize);

  normalize(nEcalShower);
  normalize(nEcalHits);

  normalize(nTrClusterX);
  normalize(nTrClusterY);
  normalize(nTrClusterL1X);
  normalize(nTrClusterL1Y);

  normalize(nTofClusterH);

  normalize(nParticle);
  normalize(nGoodParticle);
  normalize(nPartTrkInn);
  normalize(nPartTrk);

  for (auto &inner : nPartRate)
    std::for_each(begin(inner), end(inner), normalize);
  for (auto &inner : nPartBetaRate)
    std::for_each(begin(inner), end(inner), normalize);
  for (auto &inner : nPartInnRate)
    std::for_each(begin(inner), end(inner), normalize);
  std::for_each(begin(nPartAboveCutoff), end(nPartAboveCutoff), normalize);
  std::for_each(begin(nPartBelowCutoff), end(nPartBelowCutoff), normalize);

  normalize(nTrk);
  normalize(nTrkInn);

  for (auto &inner : nTrkRate)
    std::for_each(begin(inner), end(inner), normalize);
  for (auto &inner : nTrkInnRate)
    std::for_each(begin(inner), end(inner), normalize);

  normalize(nACC);
  std::for_each(begin(TrigPhysBPatt), end(TrigPhysBPatt), normalize);
  std::for_each(begin(TrigJMembPatt), end(TrigJMembPatt), normalize);
  std::for_each(begin(TrigRates), end(TrigRates), normalize);

  std::for_each(begin(TofFlags), end(TofFlags), normalize);

  normalize(TofRaw_nstdc);
  normalize(TofRaw_nftdc);
  normalize(TofRaw_nsumh);
  normalize(TofRaw_nsumsh);

  normalize(AntiRaw_ntdct);
  normalize(AntiRaw_nftdc);

  normalize(nRoomError);

  m_processedEvents = 0;
}

void RTIInfo::AccumulateEventVariables(AMSEventR *evPtr) {
  m_processedEvents++;

  // TRD
  nTRDTrack += evPtr->NTrdTrack();
  nTRDHits += evPtr->nTrdRawHit();

  std::for_each(begin(evPtr->TrdCluster()), end(evPtr->TrdCluster()),
                [this](const TrdClusterR &cl) { nTRDLayerHits[cl.Layer]++; });

  // ECAL
  nEcalShower += evPtr->NEcalShower();
  nEcalHits += evPtr->nEcalHit();

  // TrCluster
  std::for_each(begin(evPtr->TrRawCluster()), end(evPtr->TrRawCluster()), [this](TrRawClusterR &cl) {
    if (cl.GetSide() == 0) {
      nTrClusterX++;
      if (cl.GetLayerJ() == 1)
        nTrClusterL1X++;
    } else {
      nTrClusterY++;
      if (cl.GetLayerJ() == 1)
        nTrClusterL1Y++;
    }
  });

  // TOF
  nTofClusterH += evPtr->nTofClusterH();

  // Particle rates
  std::for_each(begin(evPtr->Particle()), end(evPtr->Particle()), [this](ParticleR &part) {
    nParticle++;

    if (GoodParticle(&part)) {
      nGoodParticle++;

      constexpr int refit = 31;
      constexpr int fit_algo = 1;
      int id_inner = part.pTrTrack()->iTrTrackPar(fit_algo, 3, refit);
      if (id_inner >= 0) {
        nPartTrkInn++;
        float rig = part.pTrTrack()->GetRigidity(id_inner);
        float q = part.pTrTrack()->GetQ(part.pBetaH()->GetBeta());
        int rbin = FindBin(rig, rigbins);
        int qbin = rig < 0 ? 0 : FindBin(q, qbins);
        if (qbin >= 0 && rbin >= 0)
          nPartInnRate[qbin][rbin]++;
      }
      int id_max = part.pTrTrack()->iTrTrackPar(fit_algo, 0, refit);
      if (id_max >= 0) {
        nPartTrk++;
        float rig = part.pTrTrack()->GetRigidity(id_inner);
        float q = part.pTrTrack()->GetQ(part.pBetaH()->GetBeta());
        int rbin = FindBin(rig, rigbins);
        int qbin = rig < 0 ? 0 : FindBin(q, qbins);
        if (qbin >= 0 && rbin >= 0)
          nPartRate[qbin][rbin]++;

        int bbin = FindBin(part.pBetaH()->GetBeta(), betabins);
        nPartBetaRate[qbin][bbin]++;

        float maxrigcutoff = max(fabs(MaxIGRFCutoff[3][0]), fabs(MaxIGRFCutoff[3][1])); // 40 deg
        float minrigcutoff = min(fabs(MinIGRFCutoff[1][0]), fabs(MinIGRFCutoff[1][1])); // 40 deg
        if (qbin >= 0 && fabs(rig) < minrigcutoff)
          nPartBelowCutoff[qbin]++;
        if (qbin >= 0 && fabs(rig) > maxrigcutoff)
          nPartAboveCutoff[qbin]++;
      }
    }
  });

  std::for_each(begin(evPtr->TrTrack()), end(evPtr->TrTrack()), [this](TrTrackR &trtrack) {
    constexpr int refit = 31;
    constexpr int fit_algo = 1; // Kalman fitter

    int id_inner = trtrack.iTrTrackPar(fit_algo, 3, refit);
    if (id_inner >= 0) {
      nTrkInn++;
      float rig = trtrack.GetRigidity(id_inner);
      float q = trtrack.GetQ(1.0);
      int rbin = FindBin(rig, rigbins);
      int qbin = rig < 0 ? 0 : FindBin(q, qbins);
      if (qbin >= 0 && rbin >= 0)
        nTrkInnRate[qbin][rbin]++;
    }
    int id_max = trtrack.iTrTrackPar(fit_algo, 0, refit);
    if (id_max >= 0) {
      nTrk++;
      float rig = trtrack.GetRigidity(id_max);
      float q = trtrack.GetQ(1.0);
      int rbin = FindBin(rig, rigbins);
      int qbin = rig < 0 ? 0 : FindBin(q, qbins);
      if (qbin >= 0 && rbin >= 0)
        nTrkRate[qbin][rbin]++;
    }
  });

  auto *trig = evPtr->pLevel1(0);
  if (trig) {
    trig->RestorePhysBPat();
    auto bpatt = std::bitset<8>(trig->PhysBPatt);
    for (size_t ipatt = 0; ipatt < 8; ipatt++)
      if (bpatt[ipatt])
        TrigPhysBPatt[ipatt]++;

    auto jpatt = std::bitset<16>(trig->JMembPatt);
    for (size_t ipatt = 0; ipatt < 16; ipatt++)
      if (jpatt[ipatt])
        TrigJMembPatt[ipatt]++;

    for (size_t i = 0; i < 19; i++)
      TrigRates[i] += (float)trig->TrigRates[i];

    auto antipatt = std::bitset<8>(trig->AntiPatt);
    for (size_t iacc = 0; iacc < 8; iacc++)
      if (antipatt[iacc])
        nACC++;

    if (trig->TofFlag1 < 0)
      TofFlags[15]++; // No FTC
    else if (trig->TofFlag1 < 15)
      TofFlags[trig->TofFlag1]++;

    std::for_each(begin(evPtr->TofRawSide()), end(evPtr->TofRawSide()), [this](TofRawSideR &tofraw) {
      TofRaw_nstdc += tofraw.nstdc;   // Low Threshold History
      TofRaw_nftdc += tofraw.nftdc;   // Fast Trigger History
      TofRaw_nsumh += tofraw.nsumh;   // High Threshold History
      TofRaw_nsumsh += tofraw.nsumsh; // Super-High Threshold History
    });

    std::for_each(begin(evPtr->AntiRawSide()), end(evPtr->AntiRawSide()), [this](AntiRawSideR &antiraw) {
      AntiRaw_ntdct += antiraw.ntdct; // Low Threshold History
      AntiRaw_nftdc += antiraw.nftdc; // Fast Trigger History
    });

    // Check for room errors
    auto *daq = evPtr->pDaqEvent(0);
    if (daq) {
      bool jinj_room_error = false;
      for (int ijinj = 0; ijinj < 4; ijinj++) {
        if ((daq->JINJStatus[ijinj] >> 8) & 0x2) {
          for (int ijinf = 0; ijinf < 24; ijinf++) {
            if ((((daq->JError[ijinf]) >> 3) == 0x1A) || (((daq->JError[ijinf]) >> 3) == 0x12))
              jinj_room_error = true;
          }
        }
      }
      if (jinj_room_error)
        nRoomError++;
    }
  }
}

} // namespace NAIA