RichFill.cpp

Go to the documentation of this file.

#include "Containers/Rich.h"

#include "RichCharge.h"
#include "richidOff.h"
#include "richradidOff.h"
#include "richrecOff.h"

namespace NAIA {

bool RichBaseData::Fill(RichRingR *ringPtr, RichRingBR *ringBPtr) {

  if (!ringPtr && !ringBPtr)
    return false;

  RichBaseDataR data;

  if (ringPtr) {
    data.isNaF = ringPtr->IsNaF();

    data.beta[Rich::BetaType::CIEMAT] = ringPtr->getBeta();
    data.beta[Rich::BetaType::SIMPLE] = ringPtr->Beta;
    data.beta[Rich::BetaType::REFIT] = ringPtr->BetaRefit;
    data.beta[Rich::BetaType::LIP] = 0;

    data.ebeta[Rich::BetaType::CIEMAT] = ringPtr->getBetaError();
    data.ebeta[Rich::BetaType::SIMPLE] = ringPtr->ErrorBeta;
    data.ebeta[Rich::BetaType::REFIT] = ringPtr->getBetaError();
    data.ebeta[Rich::BetaType::LIP] = 0;

    data.isGood = ringPtr->IsGood();
    data.Pos[0] = ringPtr->getTrackEmissionPoint()[0];
    data.Pos[1] = ringPtr->getTrackEmissionPoint()[1];
    data.hasBetaLIP = false;
  }
  if (ringBPtr) {
    data.beta[Rich::BetaType::LIP] = ringBPtr->Beta;
    data.ebeta[Rich::BetaType::LIP] = 0;
    data.hasBetaLIP = true;
  }

  m_beta.push_back(data);

  return true;
}

bool RichPlusData::Fill(AMSEventR *evPtr, RichRingR *ringPtr, RichRingBR *ringPtrB, TrTrackR *trackPtr) {
  RichAdditionalData data;

  data.NHitsTot = evPtr->NRichHit();
  data.NPETotUncorr = RichHitR::getCollectedPhotoElectrons();
  if (ringPtr) {
    data.NHitsUsed = ringPtr->getUsedHits();
    data.NPEExpUncorr = ringPtr->getExpectedPhotoElectrons(false);
    data.NPEExpCorr = ringPtr->getExpectedPhotoElectrons();
    data.NPEMeasUncorr = ringPtr->getPhotoElectrons(false);
    data.NPEMeasCorr = ringPtr->getPhotoElectrons();
    data.Probability = ringPtr->getProb();
    data.RingPMTs = ringPtr->NpColPMT.size();
    data.RingPMTs2 = ringPtr->getPMTs();
    data.RichPMTs = RichHitR::getPMTs() - RichHitR::getPMTs(false);
    data.ChargeConsistency = ringPtr->getPMTChargeConsistency();
    data.RichCharge = ringPtr->getCharge2Estimate();
    data.RichCharge2 = ringPtr->NpExp > 0 ? sqrt(ringPtr->NpCol / ringPtr->NpExp) : -1; // expected   collected

    FillHyp1Data(data, evPtr, ringPtr, trackPtr);
  }
  if (ringPtrB)
    data.LIP_ringProb = ringPtrB->ProbKolm;

  m_data.push_back(data);
  return true;
}

void RichPlusData::FillHyp1Data(RichAdditionalData &data, AMSEventR *evPtr, RichRingR *ringPtr, TrTrackR *trackPtr) {
  // get the PMT calibration, if available
  RichPMTCalib *pmt_calib = (!m_isMC) ? RichPMTCalib::getHead() : nullptr;

  struct pmt_info {
    int pmt = -1;
    int nhit = 0;
    float np = 0;
    std::array<float, 3> coo{0, 0, 0};
    float dist = 0;
  };

  std::map<unsigned int, pmt_info> pmt_info_map;
  std::vector<int> crossed_pmt[2]; // primary/secondary
  std::vector<RichHitR *> rich_hit_on_ring_list;
  std::map<int, RichHitR *> rich_hit_offline_map;

  // fill pmt info map
  for (const RichHitR &hit : evPtr->RichHit()) {
    int pmt = int(hit.Channel / 16);

    // discard bad PMTs
    if ((pmt_calib) && (find(begin(pmt_calib->BadPMTs), end(pmt_calib->BadPMTs), pmt) != end(pmt_calib->BadPMTs)))
      continue;

    float np = hit.Npe;
    pmt_info_map[pmt].pmt = pmt;
    pmt_info_map[pmt].nhit++;
    pmt_info_map[pmt].np += np;
    for (int icoo = 0; icoo < 3; icoo++)
      pmt_info_map[pmt].coo[icoo] += np * hit.Coo[icoo];
  }

  for (auto &it : pmt_info_map) {
    for (int icoo = 0; icoo < 3; icoo++)
      it.second.coo[icoo] /= it.second.np;
    if (trackPtr) {
      double rigidity;
      AMSPlane plm(TVector3(0, 0, it.second.coo[2]), TVector3(1, 0, 0), TVector3(0, 1, 0));
      trackPtr->Interpolate(plm, rigidity);
      TVector3 pglob = plm.getP();
      it.second.dist = sqrt(pow(it.second.coo[0] - pglob.x(), 2.) + pow(it.second.coo[1] - pglob.y(), 2.));
    }
    if (it.second.np < 5)
      continue;
    crossed_pmt[(it.second.dist < 3.5) ? 0 : 1].push_back(it.first);
  }
  int index = 0;
  std::vector<std::pair<int, pmt_info>> pmt_info_vec;
  for (const auto &it : pmt_info_map) {
    pmt_info_vec.emplace_back(it.first, it.second);
  }
  std::sort(begin(pmt_info_vec), end(pmt_info_vec),
            [](const std::pair<int, pmt_info> &elem1, const std::pair<int, pmt_info> &elem2) {
              return elem1.second.np > elem2.second.np;
            });
  std::array<float, 5> pmt_np_uncorr;
  std::array<float, 5> pmt_dist;
  for (const auto &it : pmt_info_vec) {
    if (index > 4)
      continue;
    pmt_np_uncorr[index] = it.second.np;
    pmt_dist[index] = it.second.dist;
    index++;
  }

  // selected ring
  for (int ihit = 0; ihit < ringPtr->getUsedHits(false); ihit++)
    rich_hit_on_ring_list.push_back(evPtr->pRichHit(ringPtr->iRichHit(ihit)));
  // beta hypothesys
  int which = 0;
  for (RichOffline::RichRawEvent *hit = new RichOffline::RichRawEvent(evPtr); hit; hit = hit->next(), which++) {
    rich_hit_offline_map.insert(std::make_pair(which, hit->getpointer()));
  }
  // loop on beta hypothesys
  for (int i = 0; i < ringPtr->RawBetas(); i++) {
    int ihit = ringPtr->HitBeta(i);
    RichHitR *rich_hit = rich_hit_offline_map[ihit];
    if (!rich_hit)
      continue;

    int pmt = int(rich_hit->Channel / 16);
    bool is_bad_pmt =
        (pmt_calib) && (find(begin(pmt_calib->BadPMTs), end(pmt_calib->BadPMTs), pmt) != end(pmt_calib->BadPMTs));

    // we don't have the offline occupancy map used in the dbar analysis. Disable this for now
    bool is_bad_occupancy = false;
    // bool is_bad_occupancy = (pRichOcc) ? (!pRichOcc->IsGood(rich_hit->Channel)) : true;

    bool is_crossed_pri = find(begin(crossed_pmt[0]), end(crossed_pmt[0]), pmt) != end(crossed_pmt[0]);
    bool is_crossed_sec = find(begin(crossed_pmt[1]), end(crossed_pmt[1]), pmt) != end(crossed_pmt[1]);
    bool is_not_in_selected_ring =
        (find(begin(rich_hit_on_ring_list), end(rich_hit_on_ring_list), rich_hit) == end(rich_hit_on_ring_list));
    int status2 = ((is_bad_pmt) ? 0x1 : 0) + ((is_bad_occupancy) ? 0x2 : 0) + ((is_crossed_pri) ? 0x4 : 0) +
                  ((is_crossed_sec) ? 0x8 : 0) + ((is_not_in_selected_ring) ? 0x10 : 0);

    // integrate beta = 1 hypothesys with only good hits
    if ((rich_hit->Status >> 31) & 0x1)
      continue;
    if ((status2 & 0xf) != 0)
      continue;

    int kind_of_tile = ringPtr->IsNaF();
    std::array<float, 2> simple_resolution{1.2e-3, 4e-3};
    for (int ibeta = 0; ibeta < 2; ibeta++) {
      if (std::fabs(ringPtr->RawBeta(i, ibeta) - 1) > 3. * simple_resolution[kind_of_tile])
        continue;

      data.NHitsHyp1Tot[ibeta]++;
      if ((ringPtr->UsedBeta(i) < 0) || (ringPtr->UsedBeta(i) > 1)) {
        data.NHitsHyp1OutOfRing[ibeta]++;
        data.NPEHyp1UncorrTot += rich_hit->Npe;
      }
    }
  }

  // Apparently re-clustering requires a valid TrTrackR. I don't know why
  // this wasn't the case in the DBar ntuple
  if (trackPtr) {
    unsigned int nclus = ringPtr->ClusterizeZ1();
    int dummy1;
    float dummy2, mean;
    if (nclus > 0) {
      for (int iclu = 0; iclu < std::min(10u, nclus); iclu++) {
        ringPtr->GetClusters(iclu, dummy1, mean, dummy2);
        if (std::fabs(mean - ringPtr->getBeta()) > 0.01)
          data.NBadClusters++;
      }
    }
  }

  for (int is = 0; is < 5; is++) {
    if ((pmt_np_uncorr[is] > 5) && (pmt_dist[is] > 3.5))
      data.NSecondaryHits++;
  }
}

} // namespace NAIA