KVIDZAFromZGrid.cpp

//Created by KVClassFactory on Tue Mar  8 10:00:16 2016
//Author: Diego Gruyer
 
#include "KVIDZAFromZGrid.h"
#include "TMultiGraph.h"
#include "KVIDZALine.h"
#include "TCanvas.h"
 
ClassImp(KVIDZAFromZGrid)
ClassImp(interval)
ClassImp(interval_set)
 
 
 
 
KVIDZAFromZGrid::KVIDZAFromZGrid()
{
   // Default constructor
   // Grid is declared as a 'ZOnlyGrid' by default (this is internal mechanics)
 
   init();
   fFits.SetOwner();
   SetOnlyZId();
   fIgnoreMassID = false;
}
 
 
 
 
 
void KVIDZAFromZGrid::Copy(TObject& obj) const
{
   // This method copies the current state of 'this' object into 'obj'
   // You should add here any member variables, for example:
   //    (supposing a member variable KVIDZAFromZGrid::fToto)
   //    CastedObj.fToto = fToto;
   // or
   //    CastedObj.SetToto( GetToto() );
 
   KVIDZAGrid::Copy(obj);
   KVIDZAFromZGrid& i = (KVIDZAFromZGrid&)obj;
   i.fPIDRange = fPIDRange;
   i.fPIDRangeZList = fPIDRangeZList;
   fTables.Copy(i.fTables);
   fFits.Copy(i.fFits);
   i.fIgnoreMassID = fIgnoreMassID;
}
 
 
 
 
 
void KVIDZAFromZGrid::ReadFromAsciiFile(std::ifstream& gridfile)
{
   fPIDRange = kFALSE;
   KVIDZAGrid::ReadFromAsciiFile(gridfile);
 
   if (GetParameters()->HasParameter("PIDRANGE")) {
      fPIDRange = kTRUE;
      LoadPIDRanges();
   }
 
   // if <PARAMETER> IgnoreMassID=1 appears in file, we are only using the PID intervals to clean
   // up a messy de-e plot, not to give mass identification. particles will only be identified in Z.
   if (GetParameters()->HasParameter("IgnoreMassID") && GetParameters()->GetIntValue("IgnoreMassID") == 1)
      fIgnoreMassID = true;
   else
      fIgnoreMassID = false;
 
   check_pidranges_and_massfits();
}
 
 
 
 
void KVIDZAFromZGrid::WriteToAsciiFile(std::ofstream& gridfile)
{
   ExportToGrid();
   KVIDZAGrid::WriteToAsciiFile(gridfile);
}
 
 
 
 
void KVIDZAFromZGrid::LoadPIDRanges()
{
   fPIDRangeZList.Clear();
   KVIDentifier* id = 0;
   TIter it(GetIdentifiers());
   while ((id = (KVIDentifier*)it())) {
      int zz = id->GetZ();
      if (!GetParameters()->HasParameter(Form("PIDRANGE%d", zz))) continue;
      KVString mes = GetParameters()->GetStringValue(Form("PIDRANGE%d", zz));
      if (mes.IsWhitespace()) continue;
      int type = (mes.Contains(",") ? 2 : 1);
      interval_set* itv = new interval_set(zz, type);
      itv->SetName(GetName());
      mes.Begin("|");
      while (!mes.End()) {
         KVString tmp = mes.Next();
         tmp.Begin(":");
         int aa = tmp.Next().Atoi();
         KVString val = tmp.Next();
         double pidmin, pidmax, pid;
         if (type == 1) itv->add(aa, val.Atof());
         else if (type == 2) {
            val.Begin(",");
            pidmin = val.Next().Atof();
            pid = val.Next().Atof();
            pidmax = val.Next().Atof();
            itv->add(aa, pid, pidmin, pidmax);
//            itv->add(aa, pid, pid-0.02, pid+0.02);
         }
      }
      fPIDRangeZList.Add(zz);
      fTables.Add(itv);
   }
   fPIDRange = kTRUE;
   GetParameters()->SetValue("PIDRANGE", fPIDRangeZList.AsString());
}
 
 
 
 
void KVIDZAFromZGrid::ResetPIDRanges()
{
   fTables.Clear();
   fPIDRange = kFALSE;
   fPIDRangeZList.Clear();
}
 
 
 
 
void KVIDZAFromZGrid::ReloadPIDRanges()
{
   fTables.Clear();
   LoadPIDRanges();
}
 
 
 
 
interval_set* KVIDZAFromZGrid::GetIntervalSet(int zint) const
{
   interval_set* itv = 0;
   TIter it(&fTables);
   while ((itv = (interval_set*)it())) if (itv->GetZ() == zint) return itv;
   return 0;
}
 
 
 
 
void KVIDZAFromZGrid::AddIntervalSet(interval_set* ivs)
{
   // add an interval set to the grid, updating the corresponding parameters
   fPIDRangeZList.Add(ivs->GetZ());
   fTables.Add(ivs);
}
 
 
 
 
void KVIDZAFromZGrid::RemoveIntervalSet(int zint)
{
   // Remove interval set for given Z from grid
   RemoveIntervalSet(GetIntervalSet(zint));
}
 
 
 
 
void KVIDZAFromZGrid::RemoveIntervalSet(interval_set* itvs)
{
   // Remove interval set from grid, update parameters related to interval sets
   auto z_int = itvs->GetZ();
   fPIDRangeZList.Remove(z_int); // remove Z from PIDRANGE list
   fTables.Remove(itvs);
   delete itvs;
}
 
 
 
 
void KVIDZAFromZGrid::PrintPIDLimits()
{
//    ((interval_set*)fTables.At(12))->fIntervals.ls();
 
//   for (int zz = fZminInt; zz <= fZmaxInt; zz++) {
//      Info("PrintPIDLimits", "Z=%2d    [%.4lf  %.4lf]", zz, ((interval_set*)fTables.At(zz - fZminInt))->fPIDmins.at(0),
//           ((interval_set*)fTables.At(zz - fZminInt))->fPIDmaxs.at(((interval_set*)fTables.At(zz - fZminInt))->fNPIDs - 1));
//      }
}
 
 
 
 
void KVIDZAFromZGrid::ClearPIDIntervals()
{
   fPIDRangeZList.Clear();
   if (GetParameters()->HasParameter("PIDRANGE")) GetParameters()->RemoveParameter("PIDRANGE");
   for (int ii = 1; ii < 50; ii++) {
      if (GetParameters()->HasParameter(Form("PIDRANGE%d", ii))) GetParameters()->RemoveParameter(Form("PIDRANGE%d", ii));
   }
}
 
 
 
 
int KVIDZAFromZGrid::is_inside(double pid) const
{
   // Look for a set of mass-interval definitions in which the given PID
   // falls (PID from linearisation of Z identification).
   //
   // In principle this should be the set corresponding to Z=nint(PID),
   // but if not Z+/-1 are also tried.
   //
   // Returns the value of Z for the set found (or 0 if no set found)
 
   int zint = TMath::Nint(pid);
   interval_set* it = GetIntervalSet(zint);
   if (it) {
      if (it->is_inside(pid)) return zint;
      else if (it->is_above(pid)) {
 
         it = GetIntervalSet(zint + 1);
         if (it && it->is_inside(pid)) return zint + 1;
         else return 0;
      }
      else {
         it = GetIntervalSet(zint - 1);
         if (it && it->is_inside(pid)) return zint - 1;
         else return 0;
      }
   }
   else return 0;
}
 
 
 
 
void KVIDZAFromZGrid::Initialize()
{
   // General initialisation method for identification grid.
   //
   // This method MUST be called once before using the grid for identifications.
   //
   //  + The ID lines are sorted.
   //  + The natural line widths of all ID lines are calculated.
   //  + The line with the largest Z (Zmax line) is found.
 
   SetOnlyZId();
   KVIDZAGrid::Initialize();
 
   // set up mass fits (if any)
   fFits.Clear();
   if (GetParameters()->HasStringParameter("MASSFITS")) {
      KVNumberList zlist(GetParameters()->GetStringValue("MASSFITS"));
      for (auto z : zlist) {
         auto massfit = GetParameters()->GetTStringValue(Form("MASSFIT_%d", z));
         massfit.ReplaceAll(":", "=");
         KVNameValueList fitparams;
         fitparams.Set(massfit);
         fFits.Add(new KVMultiGaussIsotopeFit(z, fitparams));
      }
   }
}
 
 
 
 
void KVIDZAFromZGrid::check_pidranges_and_massfits()
{
   // Check for inconsistencies between 'PIDRANGE' and 'MASSFIT' parameters
   //
   // We also check & signal the inclusion of isotopes which are only resonances:
   //  5Li, 6Be, 8Be, 9B, etc.
 
   if(fPar.HasParameter("MASSFITS"))
   {
      KVNumberList massfits(GetParameters()->GetStringValue("MASSFITS"));
      // check same as PIDRANGE
      KVNumberList pidrange(GetParameters()->GetStringValue("PIDRANGE"));
      // exclude from massfits any missing pidrange
      auto massfits_copy = massfits;
      massfits.Inter(pidrange);
      if(massfits != massfits_copy)
      {
         KVError::Error(this, "KVIDZAFromZGrid::check_pidranges_and_massfits",
                 "Grid %s has inconsistent MASSFITS (%s) and PIDRANGE (%s) - correcting...",
                 GetName(), massfits_copy.AsString(), pidrange.AsString());
         GetParameters()->SetValue("MASSFITS", massfits.AsString());
         // remove any fits which do not have a corresponding pidrange
         for(auto z : massfits)
         {
            if(!GetParameters()->HasParameter(Form("PIDRANGE%d",z)))
            {
               KVError::Warning(this, "KVIDZAFromZGrid::check_pidranges_and_massfits",
                       "Grid %s : removing fit for Z=%d...",
                       GetName(), z);
               GetParameters()->RemoveParameter(Form("MASSFITS_%d",z));
            }
         }
      }
      // for each MASSFIT_Z, check the 'Alist' compared to the masses given in PIDRANGEZ
      // we signal any case where the list of fitted masses is not a subset of the pidrange masses
      massfits = GetParameters()->GetStringValue("MASSFITS");
      for (auto z : massfits) {
         if(pidrange.Contains(z))
         {
            auto massfit = GetParameters()->GetTStringValue(Form("MASSFIT_%d", z));
            massfit.ReplaceAll(":", "=");
            KVNameValueList fitparams;
            fitparams.Set(massfit);
            KVNumberList alist(fitparams.GetStringValue("Alist"));
            for(auto A : alist)
            {
               if(KVNucleus(z,A).IsResonance())
                  KVError::Error(this, "check_pidranges_and_massfits",
                                  "Grid %s includes resonance %s", GetName(), KVNucleus(z,A).GetSymbol() );
            }
            KVNumberList pidalist;
            KVString pidr = GetParameters()->GetStringValue(Form("PIDRANGE%d", z));
            pidr.Begin("|");
            while(!pidr.End())
            {
               auto pipi = pidr.Next();
               pipi.Begin(":");
               pidalist.Add(pipi.Next().Atoi());
            }
            auto alist_orig = alist;
            alist.Inter(pidalist);
            if(alist != alist_orig)
            {
               KVError::Warning(this, "check_pidranges_and_massfits",
                       "Grid %s has inconsistent masses for Z=%d for MASSFIT (%s) and PIDRANGE (%s)",
                       GetName(), z, alist_orig.AsString(), pidalist.AsString());
            }
         }
      }
   }
}
 
 
 
 
bool KVIDZAFromZGrid::MassIdentificationFromMultiGaussFit(KVMultiGaussIsotopeFit* fitfunc, KVIdentificationResult* idr) const
{
   double P;
   auto A = fitfunc->GetA(idr->PID, P);
   if (A) {
      idr->A = *A;
      idr->PID = fitfunc->GetInterpolatedA(idr->PID);
      if (P > 0.5) idr->IDquality = KVIDZAGrid::kICODE0; // probability of A is >50%
      else idr->IDquality = KVIDZAGrid::kICODE3;// OK, slight ambiguity of A
   }
   else {
      // no A returned => background noise
      idr->IDquality = KVIDZAGrid::kICODE5;
      return false;
   }
   return true;
}
 
 
 
 
void KVIDZAFromZGrid::Identify(Double_t x, Double_t y, KVIdentificationResult* idr) const
{
   // Fill the KVIdentificationResult object with the results of identification for point (x,y)
   // corresponding to some physically measured quantities related to a reconstructed nucleus.
   //
   // If identification is successful, idr->IDOK = true.
   // In this case, idr->Aident and idr->Zident indicate whether isotopic or only Z identification
   // was acheived.
   //
   // In case of unsuccessful identification, idr->IDOK = false,
   // BUT idr->Zident and/or idr->Aident may be true: this is to indicate which kind of
   // identification was attempted but failed (this changes the meaning of the quality code)
 
   idr->Aident = idr->Zident = kFALSE;
 
   KVIDZAGrid::Identify(x, y, idr);
   idr->Zident = kTRUE; // meaning Z identification was attempted, even if it failed
   if (!idr->IDOK) return;
 
   bool have_pid_range_for_Z = fPIDRange && fPIDRangeZList.Contains(idr->Z);
   auto mass_fit_for_Z = GetMultiGaussFit(idr->Z);
   bool have_mass_fit_for_Z = (mass_fit_for_Z != nullptr);
   bool mass_id_success = false;
 
   if ((have_mass_fit_for_Z || have_pid_range_for_Z)
         && (!fHasMassIDRegion || idr->HasFlag(GetName(), "MassID"))) { // if a mass ID region is defined, we must be inside it
      // try mass identification
      if (have_mass_fit_for_Z)
         mass_id_success = MassIdentificationFromMultiGaussFit(mass_fit_for_Z, idr);
      else
         mass_id_success = (DeduceAfromPID(idr) > 0);
      if (mass_id_success) {
         // mass identification was at least attempted
         // make sure grid's quality code is consistent with KVIdentificationResult
         const_cast<KVIDZAFromZGrid*>(this)->fICode = idr->IDquality;
         idr->Aident = kTRUE; // meaning A identification was attempted, even if it failed
      }
      else {
         // the pid falls outside of any mass ranges for a Z which has assigned isotopes
         // therefore although the Z identification was good, we cannot consider this
         // particle to be identified
         const_cast<KVIDZAFromZGrid*>(this)->fICode = kICODE4;
         idr->IDquality = fICode; // otherwise identfication result quality code is not coherent with comment (see below)
      }
      idr->IDOK = (fICode < kICODE4);
   }
 
   // ignore isotopic successful isotopic identification if fIgnoreMassID=true
   if (fIgnoreMassID && idr->IDOK && idr->Aident) idr->Aident = false;
 
   // set comments in identification result
   switch (fICode) {
      case kICODE0:
         idr->SetComment("ok");
         break;
      case kICODE1:
         if (mass_id_success) idr->SetComment("slight ambiguity of A, which could be larger");
         else idr->SetComment("slight ambiguity of Z, which could be larger");
         break;
      case kICODE2:
         if (mass_id_success) idr->SetComment("slight ambiguity of A, which could be smaller");
         else idr->SetComment("slight ambiguity of Z, which could be smaller");
         break;
      case kICODE3:
         if (mass_id_success) idr->SetComment("slight ambiguity of A, which could be larger or smaller");
         else idr->SetComment("slight ambiguity of Z, which could be larger or smaller");
         break;
      case kICODE4:
         if (mass_id_success) idr->SetComment("point is outside of mass identification range");
         else idr->SetComment("Z identification correct but no mass identification");
         break;
      case kICODE5:
         if (mass_id_success) idr->SetComment("point is in between two isotopes A & A+2 (e.g. 5He, 8Be, 9B)");
         else idr->SetComment("point is in between two lines of different Z, too far from either to be considered well-identified");
         break;
      case kICODE6:
         idr->SetComment("(x,y) is below first line in grid");
         break;
      case kICODE7:
         idr->SetComment("(x,y) is above last line in grid");
         break;
      default:
         idr->SetComment("no identification: (x,y) out of range covered by grid");
   }
}
 
 
 
 
 
double KVIDZAFromZGrid::DeduceAfromPID(KVIdentificationResult* idr) const
{
   // First look for a set of mass intervals in which the PID of the identification result falls,
   // if there is one (see KVIDZAFromZGrid::is_inside).
   // If an interval set is found for a Z different to the original identification, idr->Z is changed.
   // Then call interval_set::eval for the mass interval for this Z.
 
   int zint = is_inside(idr->PID);
   if (!zint) return -1;
   if (zint != idr->Z) idr->Z = zint;
 
   double res = 0.;
   interval_set* it = GetIntervalSet(zint);
   if (it) res = it->eval(idr);
   return res;
}
 
 
 
 
 
void KVIDZAFromZGrid::ExportToGrid()
{
   ClearPIDIntervals();
 
   interval_set* itvs = 0;
   TIter npid(GetIntervalSets());
   while ((itvs = (interval_set*)npid())) {
      if (!itvs->GetNPID()) continue;
      fPIDRangeZList.Add(itvs->GetZ());
   }
   GetParameters()->SetValue("PIDRANGE", fPIDRangeZList.AsString());
 
   itvs = 0;
   TIter next(GetIntervalSets());
   while ((itvs = (interval_set*)next())) {
      if (!itvs->GetNPID()) continue;
      KVString par = Form("PIDRANGE%d", itvs->GetZ());
      KVString val = "";
      interval* itv = 0;
      TIter ni(itvs->GetIntervals());
      while ((itv = (interval*)ni())) {
         val += Form("%d:%lf,%lf,%lf|", itv->GetA(), itv->GetPIDmin(), itv->GetPID(), itv->GetPIDmax());
         if(KVNucleus(itvs->GetZ(),itv->GetA()).IsResonance())
            KVError::Error(this, "ExportToGrid", "Grid %s has interval set for Z=%d including resonance %s",
                           GetName(), itvs->GetZ(), KVNucleus(itvs->GetZ(),itv->GetA()).GetSymbol());
      }
      val.Remove(val.Length() - 1);
      GetParameters()->SetValue(par.Data(), val.Data());
   }
   // now remove any massfits which do not have a corresponding pidrange
   if(GetParameters()->HasParameter("MASSFITS"))
   {
      KVNumberList massfits(GetParameters()->GetStringValue("MASSFITS"));
      auto massfits_copy = massfits;
      massfits.Inter(fPIDRangeZList);
      if(massfits != massfits_copy)
      {
         ::Info("KVIDZAFromZGrid::ExportToGrid()",
                "Grid %s : will correct MASSFITS(%s) to be same as PIDRANGE(%s)",
                GetName(), massfits.AsString(), fPIDRangeZList.AsString());
         GetParameters()->SetValue("MASSFITS", fPIDRangeZList.AsString());
         massfits = fPIDRangeZList;
      }
      // remove any fit not in MASSFITS/PIDRANGE list of Z
      for(auto z : massfits)
      {
         if(!GetParameters()->HasParameter(Form("PIDRANGE%d",z)))
         {
            ::Info("KVIDZAFromZGrid::ExportToGrid()",
                   "Grid %s : removing fit for Z=%d for which no PIDRANGE exists (%s)",
                   GetName(), z, fPIDRangeZList.AsString());
            GetParameters()->RemoveParameter(Form("MASSFITS_%d",z));
         }
      }
   }
}
 
 
 
 
 
 
double interval_set::eval(KVIdentificationResult* idr)
{
   double pid = idr->PID;
   if (pid < 0.5) return 0.;
   // calculate interpolated mass from PID
   double res = fPIDs.Eval(pid);
   int ares = 0;
 
   if (fType == KVIDZAFromZGrid::kIntType) {
 
      // look for mass interval PID is in
      // in case it falls between two intervals remember also the interval
      // immediately to the left & right of the PID
      interval* left_int(nullptr), *right_int(nullptr);
      interval* inter;
      TIter it(&fIntervals);
      while ((inter = (interval*)it())) {
         if (inter->is_inside(pid)) {
            ares = inter->GetA();
            break;
         }
         else if (inter->is_left_of(pid)) {
            left_int = inter;
         }
         else if (!right_int && inter->is_right_of(pid)) {
            right_int = inter;
         }
      }
      if (ares != 0) {
         // the PID is inside a defined mass interval
         idr->A = ares;
         idr->PID = res;
         idr->IDquality = KVIDZAGrid::kICODE0;
      }
      else {
         // the PID is not inside a defined mass interval
         //
         // * if it is in between two consecutive masses i.e. A and A+1 then it is
         //   Z- and A-identified with a slight ambiguity of A
         // * if it is in between two non-consecutive masses i.e. A and A+2 then it
         //   is not identified (e.g. 5He, 8Be, 9B)
         if (!right_int || !left_int) {
            // case where no left or right interval were found
            // to prevent from crashes but should not appen
            idr->A = ares;
            idr->PID = res;
            idr->IDquality = KVIDZAGrid::kICODE5;
         }
         else {
            int dA = right_int->GetA() - left_int->GetA();
            if (dA == 1) {
               // OK, slight ambiguity of A
               ares = TMath::Nint(res);
               idr->A = ares;
               idr->PID = res;
               idr->IDquality = KVIDZAGrid::kICODE3;
            }
            else {
               // in a hole where no isotopes should be (e.g. 5He, 8Be, 9B)
               idr->A = ares;
               idr->PID = res;
               idr->IDquality = KVIDZAGrid::kICODE5;
            }
         }
      }
   }
   else {
      ares = TMath::Nint(res);
      idr->A = ares;
      idr->PID = res;
      if (ares > fPIDs.GetX()[0] && ares < fPIDs.GetX()[fNPIDs - 1]) {
         idr->IDquality = KVIDZAGrid::kICODE0;
      }
      else {
         idr->IDquality = KVIDZAGrid::kICODE4;
      }
   }
   return res;
}
 
 
 
 
bool interval_set::is_inside(double pid)
{
   if (fType != KVIDZAFromZGrid::kIntType) return kTRUE;
 
//   Info("is_inside","min: %d max:%d npids:%d", ((interval*)fIntervals.At(0))->GetA(), ((interval*)fIntervals.At(fNPIDs-1))->GetA(), fNPIDs);
 
   if (pid > ((interval*)fIntervals.At(0))->GetPIDmin() && pid < ((interval*)fIntervals.At(fNPIDs - 1))->GetPIDmax()) return kTRUE;
   else return kFALSE;
}
 
 
 
 
bool interval_set::is_above(double pid)
{
   if (fType != KVIDZAFromZGrid::kIntType) return kTRUE;
 
   if (pid > ((interval*)fIntervals.At(fNPIDs - 1))->GetPIDmax()) return kTRUE;
   else return kFALSE;
}
 
 
 
 
 
TString interval_set::GetListOfMasses()
{
   if (!GetNPID()) return "-";
   KVNumberList alist;
   for (int ii = 0; ii < GetNPID(); ii++) alist.Add(((interval*)fIntervals.At(ii))->GetA());
   return alist.AsString();
}
 
 
 
 
interval_set::interval_set(int zz, int type)
{
   fType = type;
   fZ = zz;
   fNPIDs = 0;
}
 
 
 
 
void interval_set::add(int aa, double pid, double pidmin, double pidmax)
{
   if (fType == KVIDZAFromZGrid::kIntType && !(pid > pidmin && pid < pidmax))
   {
      // sort into correct order: pidmin, pid, pidmax
      std::map<double,int> toto;
      ++toto[pidmin];
      ++toto[pid];
      ++toto[pidmax];
      auto it = std::begin(toto);
      pidmin = it->first;
      ++it;
      pid = it->first;
      ++it;
      pidmax = it->first;
   }
 
   fPIDs.SetPoint(fNPIDs, pid, aa);
   if (fType == KVIDZAFromZGrid::kIntType) {
      if (pid) fIntervals.AddLast(new interval(fZ, aa, pid, pidmin, pidmax));
   }
   fNPIDs++;
}