KVDataQualityAudit.cpp

#include "KVDataQualityAudit.h"
#include <algorithm>
 
 
 
void KVDataQualityAudit::Add(const KVReconstructedNucleus& N)
{
   // Add this reconstructed nucleus to the audit.
 
   auto tel = telescopes.get_object<idtelescope>(N.GetIdentifyingTelescope()->GetName());
   if (!tel) {
      tel = new idtelescope(N.GetIdentifyingTelescope()->GetName());
      telescopes.Add(tel);
   }
   tel->add(N);
}
 
 
 
 
void KVDataQualityAudit::Print(Option_t* opt) const
{
   for (auto p : telescopes) {
      p->Print(opt);
   }
}
 
 
 
 
Long64_t KVDataQualityAudit::Merge(TCollection* l)
{
   // Used to merge audit objects in different ROOT files (using hadd).
   //
   // This method is called for the object in the first file, the TCollection
   // contains the objects with the same name in all the other files.
 
   TIter next(l);
   TObject* o;
   Long64_t nmerge = 0;
   while ((o = next())) {
      nmerge += merge(dynamic_cast<KVDataQualityAudit*>(o));
   }
   return nmerge;
}
 
 
 
 
KVDataQualityAudit::IDStatus KVDataQualityAudit::CanIdentify(const TString& tel_name, int Z, int& A, double E) const
{
   // Call this method to test whether a given particle could be identified in a given identification
   // telescope according to the data quality audit.
   //
   // \param[in] tel_name name of identification telescope
   // \param[in] Z atomic number \f$Z\f$ of particle
   // \param[in,out] A mass number \f$A\f$ of particle. May be changed by method.
   // \param[in] E total incident energy of particle
   //
   // \returns IDStatus object which indicates if & how well the identification would proceed
 
   if (!HasTelescope(tel_name)) // telescope does not appear in audit, not functioning
      return IDStatus(IDStatus::EStatus::telescope_not_functioning);
 
   auto idtel = GetTelescope(tel_name);
   if (!idtel->HasElement(Z)) // telescope never identifies this Z value
      return IDStatus(IDStatus::EStatus::identification_not_possible_for_this_Z);
 
   auto& elem = idtel->GetElement(Z);
   if (E < elem.GetIdThreshold()) // no identification for this energy
      return IDStatus(IDStatus::EStatus::below_threshold_for_Z_identification);
 
   if (!elem.HasIsotopes()) { // telescope does not provide isotopic identification
      // attribute default mass to particle
      A = elem.get_default_mass();
      return IDStatus(IDStatus::EStatus::Z_identified_but_not_A);
   }
 
   if (elem.HasIsotope(A)) {
      if (E > elem.GetIsotope(A).GetIdThreshold())
         return IDStatus(IDStatus::EStatus::full_isotopic_identification);
 
      // telescope can identify isotope, but we are below threshold for full identification
      int a_min_thresh;
      if (E < elem.get_minimum_isotopic_threshold(a_min_thresh)) { // no isotopic identification possible
         // attribute default mass to particle
         A = elem.get_default_mass();
         return IDStatus(IDStatus::EStatus::Z_identified_but_not_A);
      }
      // look for closest isostope with A'<A with good threshold?
      // [this is modelled on CsI R/L behaviour, where heavier isotopes have higher thresholds]
      if (a_min_thresh < A) {
         int a = A - 1;
         while (a >= a_min_thresh) {
            if (elem.HasIsotope(a) && E > elem.GetIsotope(a).GetIdThreshold()) {
               // attribute isotope mass to particle
               A = a;
               return IDStatus(IDStatus::EStatus::partial_isotopic_identification);
            }
            --a;
         }
         // we shouldn't arrive here?
         return IDStatus(IDStatus::EStatus::weird_result_1);
      }
      // no luck? now look for first isotope with threshold less than particle's energy
      auto seuil_map = elem.get_isotopes_sorted_by_threshold();
      int last_a = 0;
      for (auto& s : seuil_map) {
         if (s.first > E) {
            if (last_a > 0) {
               // use isotopic mass of previous isotope in map, its threshold was <E
               A = last_a;
               return IDStatus(IDStatus::EStatus::partial_isotopic_identification);
            }
            else
               return IDStatus(IDStatus::EStatus::weird_result_3);
         }
         last_a = s.second;
      }
      // we shouldn't arrive here?
      return IDStatus(IDStatus::EStatus::weird_result_2);
   }
 
   // particle's A was never identified by the telescope:
   //   if A > Amax or A < Amin,no identification (in between lines of different Z)
   //   if Amax > A > Amin, look for closest mass in list with good threshold?
 
   // start by filling list of candidate masses with OK thresholds
   auto seuil_map = elem.get_isotopes_sorted_by_threshold();
   std::map<int, int> ok_masses;
   for (auto& s : seuil_map) {
      if (E > s.first) ok_masses[s.second] = 1;
      else break;
   }
   // if ok_masses is empty, the particle's energy is below all isotopic thresholds
   if(ok_masses.empty())
   {
      return IDStatus(IDStatus::EStatus::Z_identified_but_not_A);
   }
   // get largest isotopically resolved mass above threshold
   auto last = ok_masses.rbegin();
   int Amax;
   if(last != ok_masses.rend())
      Amax = (*last).first;
   else
   {
      Error("CanIdentify", "problem with ok_masses map...: tel=%s Z=%d A=%d E=%f", tel_name.Data(), Z, A, E);
      std::cout << "contents of seuil_map:\n";
      for(auto& p : seuil_map)
      {
         std::cout << "seuil_map[" << p.first << "] = " << p.second << std::endl;
      }
      return IDStatus(IDStatus::EStatus::weird_result_4);
   }
   // get smallest isotopically resolved mass above threshold
   auto first = ok_masses.begin();
   auto Amin = (*first).first;
   if (A > Amax) {
      return IDStatus(IDStatus::EStatus::A_greater_than_measured_Amax);
   }
   else if (A < Amin) {
      return IDStatus(IDStatus::EStatus::A_less_than_measured_Amin);
   }
   else {
      int adiff(999), near_a(-1);
      for (auto& p : ok_masses) {
         auto diff = std::abs(p.first - A);
         if (diff < adiff) {
            adiff = diff;
            near_a = p.first;
         }
      }
      A = near_a;
      return IDStatus(IDStatus::EStatus::A_never_measured_but_between_Amin_and_Amax);
   }
   // should never get here?
   return IDStatus(IDStatus::EStatus::weird_result_4);
}
 
 
 
 
void KVDataQualityAudit::add(const KVDataQualityAudit::idtelescope* idt)
{
   // Add copy of idtelescope data to this audit
   telescopes.Add(new idtelescope(*idt));
}
 
 
 
 
Long64_t KVDataQualityAudit::merge(KVDataQualityAudit* other)
{
   // Merge other audit object with this one. Return 1 if merge is successful,
   // 0 in case of problems.
   if (!other) return 0; // possibly called for wrong object type
 
   TIter next_tel(other->GetTelescopeList());
   idtelescope* idt;
   Long64_t nmerge = 0;
   while ((idt = (idtelescope*)next_tel())) {
      if (HasTelescope(idt))
         GetTelescope(idt)->merge(idt);
      else
         add(idt);
      ++nmerge;
   }
   return nmerge;
}
 
 
 
 
void KVDataQualityAudit::isotope::add(const KVReconstructedNucleus& N)
{
   // add an isotopically identified nucleus
   //
   // for calibrated particles, we store the minimum energy (=> threshold)
   A = N.GetA();
   ++counts;
   if (N.IsCalibrated() && ((emin < 0) || (N.GetEnergy() < emin))) emin = N.GetEnergy();
}
 
 
 
 
void KVDataQualityAudit::isotope::print() const
{
   std::cout << "\t\tA=" << A << "\t\t[E>" << emin << "]" << std::endl;
}
 
 
 
 
void KVDataQualityAudit::element::add(const KVReconstructedNucleus& N)
{
   // add an element (Z-identified nucleus)
   //
   // if isotopically identified, add to list of isotopes
   //
   // for calibrated particles, we store the minimum energy (=> threshold)
   if (N.GetZ() > 0) Z = N.GetZ();
   ++counts;
   if (N.IsCalibrated() && ((emin < 0) || (N.GetEnergy() < emin))) emin = N.GetEnergy();
   if (N.IsAMeasured()) {
      isotopes[N.GetA()].add(N);
   }
   else if (N.GetA() > 0) A = N.GetA(); // default mass for Z-only identification
}
 
 
 
 
void KVDataQualityAudit::element::print() const
{
   std::cout << "\tZ=" << (int)Z << "\t\t[E>" << emin << "]";
   if (A > 0) std::cout << " default A=" << A;
   if (!HasIsotopes()) {
      std::cout << std::endl;
      return;
   }
   std::cout << "\t\tIsotopes: " << GetIsotopeList().AsString() << std::endl;
   for (auto i : isotopes) i.second.print();
}
 
 
 
 
void KVDataQualityAudit::element::merge(const KVDataQualityAudit::element& elem)
{
   // Merge the informations in element elem with those in this one:
   //
   //  - we keep the lower of the two threshold values
   //  - we sum the counts for the element
   //  - we merge the lists of isotopes of the two elements
   //  - we keep the default mass, if it has been set
 
   if (emin > 0) {
      if (elem.emin > 0) emin = std::min(emin, elem.emin);
   }
   else if (elem.emin > 0) emin = elem.emin;
   A = std::max(A, elem.A);
   counts += elem.counts;
   // look at isotopes in other's list. if any are not in our list they are added.
   // if any are in both, they are merged
   std::for_each(std::begin(elem.isotopes), std::end(elem.isotopes),
   [&](const std::pair<int, isotope>& isotop) {
      if (HasIsotope(isotop.first)) {
         isotopes[isotop.first].merge(isotop.second);
      }
      else {
         isotopes[isotop.first] = isotop.second;
      }
   }
                );
}
 
 
 
 
double KVDataQualityAudit::element::get_mean_isotopic_mass() const
{
   // calculate and return mean mass of isotopes measured for this element
 
   if (isotopes.empty()) return 0;
   double sum{0}, weight{0};
   std::for_each(std::begin(isotopes), std::end(isotopes),
   [&](const std::pair<int, isotope>& count) {
      sum += count.second.counts;
      weight += count.first * count.second.counts;
   }
                );
   if (sum > 0) return weight / sum;
   return 0;
}
 
 
 
 
int KVDataQualityAudit::element::get_max_isotopic_mass() const
{
   // return max A of isotopes measured for this element
 
   if (isotopes.empty()) return 0;
   auto last = isotopes.rbegin();
   return (*last).first;
}
 
 
 
 
int KVDataQualityAudit::element::get_min_isotopic_mass() const
{
   // return min A of isotopes measured for this element
 
   if (isotopes.empty()) return 0;
   auto first = isotopes.begin();
   return (*first).first;
}
 
 
 
 
double KVDataQualityAudit::element::get_minimum_isotopic_threshold_mev_per_nuc(int& a) const
{
   // \returns the smallest recorded threshold (in MeV/u) for isotopic identification for the element
   // \param[out] a mass number corresponding to minimum
   a = -1;
   if (isotopes.empty()) return -1;
   double minE = 9999;
   for (auto& p : isotopes) {
      double myE = p.second.emin / p.first;
      if ((myE > 0) && (myE < minE)) {
         minE = myE;
         a = p.first;
      }
   }
   if (minE < 9999) return minE;
   // uncalibrated particles: all thresholds are -1 MeV
   return -1;
}
 
 
 
 
double KVDataQualityAudit::element::get_minimum_isotopic_threshold(int& a) const
{
   // \returns the smallest recorded threshold (in MeV) for isotopic identification for the element
   // \param[out] a mass number corresponding to minimum
   a = -1;
   if (isotopes.empty()) return -1;
   double minE = 9999;
   for (auto& p : isotopes) {
      double myE = p.second.emin;
      if ((myE > 0) && (myE < minE)) {
         minE = myE;
         a = p.first;
      }
   }
   if (minE < 9999) return minE;
   // uncalibrated particles: all thresholds are -1 MeV
   return -1;
}
 
 
 
 
std::map<double, int> KVDataQualityAudit::element::get_isotopes_sorted_by_threshold() const
{
   // \returns a map sorted according to increasing identification threshold of isotope mass numbers
 
   std::map<double, int> seuils;
   for (auto& i : isotopes) {
      seuils[i.second.GetIdThreshold()] = i.first;
   }
   return seuils;
}
 
 
 
 
std::map<int, double> KVDataQualityAudit::element::get_isotopic_distribution() const
{
   // \returns probability distribution for each isotope as a map between isotope A and P(A)
 
   double sum{0};
   std::for_each(std::begin(isotopes), std::end(isotopes),
   [&](const std::pair<int, isotope>& count) {
      sum += count.second.counts;
   }
                );
   std::map<int, double> proba;
   for (auto& p : isotopes) {
      proba[p.first] = p.second.counts / sum;
   }
   return proba;
}
 
 
 
 
void KVDataQualityAudit::idtelescope::add(const KVReconstructedNucleus& N)
{
   elements[N.GetZ()].add(N);
}
 
 
 
 
void KVDataQualityAudit::idtelescope::Print(Option_t*) const
{
   std::cout << GetName() << " :" << std::endl;
   std::cout << "Elements: " << GetElementList().AsString() << std::endl;
   for (auto e : elements) e.second.print();
}
 
 
 
 
void KVDataQualityAudit::idtelescope::merge(const KVDataQualityAudit::idtelescope* idt)
{
   // Merge data in idtelescope idt with this one
 
   // look at elements in other's list. if any are not in our list they are added.
   // if any are in both, they are merged
   std::for_each(std::begin(idt->elements), std::end(idt->elements),
   [&](const std::pair<int, element>& elem) {
      if (HasElement(elem.first)) {
         elements[elem.first].merge(elem.second);
      }
      else {
         elements[elem.first] = elem.second;
      }
   }
                );
}
 
 
 
 
std::map<int, double> KVDataQualityAudit::idtelescope::get_element_distribution() const
{
   // \returns probability distribution for all elements identified by telescope as a map between Z and P(Z)
 
   double sum{0};
   std::for_each(std::begin(elements), std::end(elements),
   [&](const std::pair<int, element>& count) {
      sum += count.second.counts;
   }
                );
   std::map<int, double> proba;
   for (auto& p : elements) {
      proba[p.first] = p.second.counts / sum;
   }
   return proba;
}
 
 
 
 
double KVDataQualityAudit::idtelescope::get_mean_Z() const
{
   // \returns mean Z measured in telescope
   double sum{0}, weight{0};
   std::for_each(std::begin(elements), std::end(elements),
   [&](const std::pair<int, element>& count) {
      sum += count.second.counts;
      weight += count.second.counts * count.first;
   }
                );
   return weight / sum;
}
 
 
 
 
int KVDataQualityAudit::idtelescope::get_max_Z() const
{
   // \returns maximum Z measured in telescope
   auto last = elements.rbegin();
   return (*last).first;
}
 
 
 
int KVDataQualityAudit::idtelescope::get_max_Z_with_isotopes() const
{
   // \returns maximum Z for which isotopic masses were measured
   auto it = elements.rbegin();
   while (it != elements.rend()) {
      if ((*it).second.HasIsotopes()) return (*it).first;
      ++it;
   }
   return 0;
}
 
 
 
int KVDataQualityAudit::idtelescope::get_min_Z() const
{
   // \returns minimum Z measured in telescope
   auto first = elements.begin();
   return (*first).first;
}
 
 
 
 
double KVDataQualityAudit::idtelescope::get_mean_A() const
{
   // \returns mean A measured in telescope
   double sum{0}, weight{0};
   std::for_each(std::begin(elements), std::end(elements),
   [&](const std::pair<int, element>& count) {
      sum += count.second.counts;
      weight += count.second.counts * count.second.get_mean_isotopic_mass();
   }
                );
   return weight / sum;
}
 
 
 
 
int KVDataQualityAudit::idtelescope::get_max_A() const
{
   // \returns maximum A measured in telescope
   int max = 0;
   for (auto& e : elements) {
      max = std::max(max, e.second.get_max_isotopic_mass());
   }
   return max;
}
 
 
 
 
int KVDataQualityAudit::idtelescope::get_min_A() const
{
   // \returns smallest A measured in telescope
   int min = 999;
   for (auto& e : elements) {
      if (auto other_A = e.second.get_min_isotopic_mass())
         min = std::min(min, other_A);
   }
   return min;
}
 
 
 
 
ClassImp(KVDataQualityAudit)