KVDataQualityAudit.h

#ifndef __KVDATAQUALITYAUDIT_H
#define __KVDATAQUALITYAUDIT_H
 
#include "KVBase.h"
 
#include <KVNumberList.h>
#include <KVReconstructedNucleus.h>
#include <KVUniqueNameList.h>
 
class KVDataQualityAudit : public KVBase {
private:
   KVUniqueNameList telescopes;
 
   Long64_t merge(KVDataQualityAudit*);
 
public:
   KVDataQualityAudit() : KVDataQualityAudit("") {}
   KVDataQualityAudit(const Char_t* name, const Char_t* title = "")
      : KVBase(name, title)
   {
      telescopes.SetOwner();
   }
   virtual ~KVDataQualityAudit() {}
   void Add(const KVReconstructedNucleus& N);
   void Print(Option_t* opt = "") const override;
   const KVSeqCollection* GetTelescopeList() const
   {
      return &telescopes;
   }
   Bool_t HasTelescope(const TString& tel_name) const
   {
      return telescopes.FindObject(tel_name) != nullptr;
   }
   Long64_t Merge(TCollection*);
 
   struct isotope {
      Double_t counts{0};
      Float_t emin{-1.0};
      UShort_t A{0};
      void add(const KVReconstructedNucleus& N);
      void print() const;
      void merge(const isotope& isoto)
      {
         if (emin > 0) {
            if (isoto.emin > 0) emin = std::min(emin, isoto.emin);
         }
         else if (isoto.emin > 0) emin = isoto.emin;
         counts += isoto.counts;
      }
      Float_t GetIdThreshold() const
      {
         return emin;
      }
 
      ClassDef(isotope, 1)
   };
   struct element {
      Double_t counts{0};
      Float_t emin{-1.0};
      UShort_t A{0};// default (calculated) mass given when no isotopic identification available
      UChar_t Z{0};
      std::map<int, isotope> isotopes;
      void add(const KVReconstructedNucleus& N);
      void print() const;
      void merge(const element&);
      Bool_t HasIsotopes() const
      {
         return !isotopes.empty();
      }
      Bool_t HasIsotope(int A) const
      {
         return isotopes.find(A) != std::end(isotopes);
      }
      const isotope& GetIsotope(int A) const
      {
         try {
            return isotopes.at(A);
         }
         catch (...) {
            std::cerr << "Error in <KVDataQualityAudit::GetIsotope>: No isotope with A=" << A << std::endl;
            throw;
         }
      }
      KVNumberList GetIsotopeList() const
      {
         KVNumberList l;
         for (auto& p : isotopes) l.Add(p.first);
         return l;
      }
      double get_mean_isotopic_mass() const;
      int get_max_isotopic_mass() const;
      std::map<int, double> get_isotopic_distribution() const;
      int get_min_isotopic_mass() const;
      double get_minimum_isotopic_threshold_mev_per_nuc(int&) const;
      double get_minimum_isotopic_threshold(int&) const;
      int get_default_mass() const
      {
         return A;
      }
      Float_t GetIdThreshold() const
      {
         return emin;
      }
      std::map<double, int> get_isotopes_sorted_by_threshold() const;
 
      ClassDef(element, 1)
   };
   class idtelescope : public TNamed {
      std::map<int, element> elements;
   public:
      idtelescope(const TString& name = "") : TNamed(name, name) {}
      ROOT_COPY_CTOR(idtelescope, TNamed)
      ROOT_COPY_ASSIGN_OP(idtelescope)
      void Copy(TObject& o) const override
      {
         TNamed::Copy(o);
         idtelescope& other = (idtelescope&)o;
         other.elements = elements;
      }
      void add(const KVReconstructedNucleus& N);
      void Print(Option_t* opt = "") const override;
      void merge(const idtelescope*);
      Bool_t HasElement(int Z) const
      {
         return elements.find(Z) != std::end(elements);
      }
      const element& GetElement(int Z) const
      {
         try {
            return elements.at(Z);
         }
         catch (...) {
            std::cerr << "Error in <KVDataQualityAudit::GetElement>: No element with Z=" << Z << std::endl;
            throw;
         }
      }
      KVNumberList GetElementList() const
      {
 
         KVNumberList l;
         for (auto& p : elements) l.Add(p.first);
         return l;
      }
      std::map<int, double> get_element_distribution() const;
      double get_mean_Z() const;
      int get_max_Z() const;
      double get_mean_A() const;
      int get_max_A() const;
      int get_min_A() const;
      int get_min_Z() const;
      int get_max_Z_with_isotopes() const;
 
      ClassDefOverride(idtelescope, 1)
   };
 
   Bool_t HasTelescope(const idtelescope* idt) const
   {
      return HasTelescope(idt->GetName());
   }
   idtelescope* GetTelescope(const TString& tel_name) const
   {
      return (idtelescope*)telescopes.FindObject(tel_name);
   }
   idtelescope* GetTelescope(const idtelescope* idt) const
   {
      return GetTelescope(idt->GetName());
   }
 
   class IDStatus {
   public:
      enum class EStatus {
         telescope_not_functioning,
         below_threshold_for_Z_identification,
         identification_not_possible_for_this_Z,
         Z_identified_but_not_A,
         partial_isotopic_identification,
         A_greater_than_measured_Amax,
         A_less_than_measured_Amin,
         A_never_measured_but_between_Amin_and_Amax,
         full_isotopic_identification,
         weird_result_1,
         weird_result_2,
         weird_result_3,
         weird_result_4
      };
      IDStatus(EStatus s) : status{s} {}
      operator bool() const
      {
         return (status == EStatus::Z_identified_but_not_A)
                || (status == EStatus::full_isotopic_identification)
                || (status == EStatus::partial_isotopic_identification)
                || (status == EStatus::A_never_measured_but_between_Amin_and_Amax);
      }
      bool operator==(EStatus s) const
      {
         return status == s;
      }
   private:
      EStatus status;
   public:
      friend std::ostream& operator<<(std::ostream& os, const IDStatus& st_ob)
      {
         os << "ID Status: ";
         switch (st_ob.status) {
            case EStatus::telescope_not_functioning:
               os << "telescope not functioning";
               break;
            case EStatus::below_threshold_for_Z_identification:
               os << "below threshold for Z identification";
               break;
            case EStatus::identification_not_possible_for_this_Z:
               os << "identification not possible for this Z";
               break;
            case EStatus::Z_identified_but_not_A:
               os << "Z identified but not A";
               break;
            case EStatus::partial_isotopic_identification:
               os << "partial isotopic identification";
               break;
            case EStatus::A_greater_than_measured_Amax:
               os << "A greater than measured Amax";
               break;
            case EStatus::A_less_than_measured_Amin:
               os << "A less than measured Amin";
               break;
            case EStatus::A_never_measured_but_between_Amin_and_Amax:
               os << "A_never_measured_but_between_Amin_and_Amax";
               break;
            case EStatus::full_isotopic_identification:
               os << "full isotopic identification";
               break;
            case EStatus::weird_result_1:
               os << "weird result 1";
               break;
            case EStatus::weird_result_2:
               os << "weird result 2";
               break;
            case EStatus::weird_result_3:
               os << "weird result 3";
               break;
            case EStatus::weird_result_4:
               os << "weird result 4";
               break;
         }
         return os;
      }
   };
 
   IDStatus CanIdentify(const TString& tel_name, int Z, int& A, double E) const;
 
private:
   void add(const idtelescope*);
 
   ClassDefOverride(KVDataQualityAudit, 1) //Audit of experimental data identification and calibrations
};
 
#endif