KVIDTelescope.cpp

/***************************************************************************
$Id: KVIDTelescope.cpp,v 1.52 2009/05/05 15:54:04 franklan Exp $
Author : $Author: franklan $
                          KVIDTelescope.cpp  -  description
                             -------------------
    begin                : Wed Jun 18 2003
    copyright            : (C) 2003 by J.D Frankland
    email                : frankland@ganil.fr
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
#include "TROOT.h"
#include "KVIDTelescope.h"
#include "KVGroup.h"
#include "KVNucleus.h"
#include "KVReconstructedNucleus.h"
#include "KVIDGraph.h"
#include "KVIDGrid.h"
#include "Riostream.h"
#include "TPluginManager.h"
#include "KVMultiDetArray.h"
#include "KVDataSet.h"
#include "KVIDGridManager.h"
#include "KVIDZALine.h"
#include "KVIDCutLine.h"
#include "KVIdentificationResult.h"
#include "TMath.h"
#include "TClass.h"
#include "TH2.h"
#include "KVParticleCondition.h"
 
#include <KVDetectorSignalExpression.h>
#include <KVIDZAGrid.h>
 
using namespace std;
 
ClassImp(KVIDTelescope)
TEnv* KVIDTelescope::fgIdentificationBilan = nullptr;
 
 
 
KVIDTelescope::KVIDTelescope()
   : fGroup(nullptr)
{
   init();
}
 
 
 
 
void KVIDTelescope::init()
{
   //default init
   fDetectors.SetCleanup(kTRUE);
   fIDGrids.SetCleanup(kTRUE);
}
 
 
 
 
void KVIDTelescope::Initialize(void)
{
   // Default initialisation for ID telescopes.
   // If telescope has at least 1 grid then it is ready to identify
   // particles after initialising the grid(s) (kReadyForID=true);
   // otherwise kReadyForID is set to kFALSE, unless the current dataset (if defined)
   // has been declared to have no associated identification/calibration parameters,
   // in which case kReadyForID is by default set to kTRUE (for filtering simulations).
   //
   // In order to enable mass identification for certain telescopes without a dedicated
   // implementation (e.g. for simulating array response), put the following lines
   // in your .kvrootrc:
   //
   //   [dataset].[telescope label].MassID:   yes
   //
   // If you want to limit mass identification to certain values of Z and/or A,
   // add the following line:
   //
   //   [dataset].[telescope label].MassID.Validity:    [expression]
   //
   // where [expression] is some valid C++ boolean expression involving Z and/or A,
   // for example
   //
   //   [dataset].[telescope label].MassID.Validity:  (Z>3)&&(A<20)
   //
   //For identifications using more than one grid, the default behaviour is to try identification
   //with each grid in turn until a successful identification is obtained. The order in which
   //the grids should be tried should be specified by a variable with the following format:
   //
   //~~~~~~~~~~~~~~~~
   //[Dataset].[telescope label].GridOrder:  [Grid1],[Grid2],...
   //~~~~~~~~~~~~~~~~
 
   ResetBit(kReadyForID);
 
   // for datasets with no calib/ident infos, all id telescopes work
   if (gDataSet && !gDataSet->HasCalibIdentInfos()) {
      SetBit(kReadyForID);
   }
   else { // for datasets with calib/ident infos, we need a grid & all detectors working
      // looping over detectors to check they are working
      // if one of them is not -> set kReadyForID to false
      TIter it(GetDetectors());
      KVDetector* det = 0;
      while ((det = (KVDetector*)it())) if (!det->IsOK()) {
            ResetBit(kReadyForID);
            return;
         }
 
      if (GetIDGrid()) {
         KVIDGraph* gr;
         TIter it(GetListOfIDGrids());
         bool ok = kTRUE;
         KVUniqueNameList tmp_list;// for re-ordering grids
         bool mass_id = false;
         while ((gr = (KVIDGraph*)it())) {
            tmp_list.Add(gr);
            if (gr->HasMassIDCapability()) mass_id = true;
            gr->Initialize();
            // make sure both x & y axes' signals are well set up
            if (!fGraphCoords[gr].fVarX || !fGraphCoords[gr].fVarY) {
               ok = kFALSE;
               Warning("Initialize",
                       "ID tel. %s: grid %s has undefined VarX(%s:%p) or VarY(%s:%p) - WILL NOT USE",
                       GetName(), gr->GetName(), gr->GetVarX(), fGraphCoords[gr].fVarX, gr->GetVarY(), fGraphCoords[gr].fVarY);
            }
         }
         // set to true if at least one grid can provide mass identification
         SetHasMassID(mass_id);
         // if more than one grid, need to re-order them according to [Dataset].[telescope label].GridOrder
         if (GetListOfIDGrids()->GetEntries() > 1 && gDataSet) {
            KVString grid_list = gDataSet->GetDataSetEnv(Form("%s.GridOrder", GetLabel()));
            ok = kFALSE;
            if (grid_list == "")
               Warning("Initialize", "ID telescope %s has %d grids but no %s variable defined",
                       GetName(), GetListOfIDGrids()->GetEntries(), Form("%s.GridOrder", GetLabel()));
            else if (grid_list.GetNValues(",") != GetListOfIDGrids()->GetEntries())
               Warning("Initialize", "ID telescope %s has %d grids but %d grids appear in variable %s",
                       GetName(), GetListOfIDGrids()->GetEntries(), grid_list.GetNValues(","), Form("%s.GridOrder", GetLabel()));
            else {
               fIDGrids.Clear();
               grid_list.Begin(",");
               while (!grid_list.End()) {
                  auto gr_name = grid_list.Next();
                  auto gr_ob = tmp_list.FindObject(gr_name);
                  if (!gr_ob) {
                     Info("Initialize", "IDtel=%s grid %s missing", GetName(), gr_name.Data());
                  }
                  else {
                     fIDGrids.Add(gr_ob);
                  }
               }
               ok = kTRUE;
            }
         }
         if (ok) SetBit(kReadyForID);
      }
   }
 
   if (gDataSet) {
      SetHasMassID(gDataSet->GetDataSetEnv(Form("%s.MassID", GetLabel()), kFALSE));
      KVString valid;
      if ((valid = gDataSet->GetDataSetEnv(Form("%s.MassID.Validity", GetLabel()), "")) != "") {
         valid.ReplaceAll("Z", "_NUC_->GetZ()");
         valid.ReplaceAll("A", "_NUC_->GetA()");
         fMassIDValidity.reset(new KVParticleCondition(valid));
      }
   }
}
 
 
 
 
void KVIDTelescope::AddDetector(KVDetector* d)
{
   // Add a detector to the telescope.
   //
   // Detectors must be added in the order they will be hit by impinging particles,
   // with the last detector being the one particles stopped in the telescope will stop in.
   // i.e. dE1, dE2, ..., Eres
   //
   // Update name of telescope to "ID_[name of 1st detector]_[name of 2nd detector]_ ... _[name of last detector]"
 
   if (d) {
      fDetectors.Add(d);
      if (GetSize() > 1) {
         TString old_name = GetName();
         old_name += Form("_%s", GetDetectors()->Last()->GetName());
         SetName(old_name);
      }
      else SetName(Form("ID_%s", GetDetector(1)->GetName()));
      //d->AddIDTelescope(this);  <= caused multiple copies to exist in detector's list
   }
   else {
      Warning("AddDetector", "Called with null pointer");
   }
}
 
 
 
 
void KVIDTelescope::Print(Option_t* opt) const
{
   // print out telescope structure
   //if opt="fired" only fired detectors are printed
 
   TIter next(GetDetectors());
   KVDetector* obj;
 
   if (!strcmp(opt, "fired")) {
      while ((obj = (KVDetector*) next())) {
 
         if (obj->Fired() || obj->GetEnergy())
            obj->Print("data");
      }
   }
   else {
      cout << "\n" << opt << "Structure of KVIDTelescope object: " <<
           GetName() << " " << GetType() << endl;
      cout << opt <<
           "--------------------------------------------------------" <<
           endl;
      while ((obj = (KVDetector*) next())) {
         cout << opt << "Detector: " << obj->GetName() << endl;
      }
   }
}
 
 
 
 
 
KVDetector* KVIDTelescope::GetDetector(const Char_t* name) const
{
   // Return a pointer to the detector in the telescope with the name "name".
 
   KVDetector* tmp = (KVDetector*) GetDetectors()->FindObject(name);
   if (!tmp)
      Warning("GetDetector(const Char_t *name)",
              "Detector %s not found in telescope %s", name, GetName());
   return tmp;
}
 
 
 
 
 
KVGroup* KVIDTelescope::GetGroup() const
{
   return fGroup;
}
 
 
 
 
 
void KVIDTelescope::SetGroup(KVGroup* kvg)
{
   fGroup = kvg;
}
 
 
 
 
UInt_t KVIDTelescope::GetGroupNumber()
{
   return (GetGroup() ? GetGroup()->GetNumber() : 0);
}
 
 
 
 
Bool_t KVIDTelescope::Identify(KVIdentificationResult* idr, Double_t x, Double_t y)
{
   //Default identification method.
   //
   //Works for ID telescopes for which one or more identification grids are defined, each
   //with VARX/VARY parameters corresponding to a KVDetectorSignal or KVDetectorSignalExpression
   //associated with one or other of the detectors constituting the telescope.
   //
   //For identifications using more than one grid, the default behaviour is to try identification
   //with each grid in turn until a successful identification is obtained. The order in which
   //the grids should be tried should be specified by a variable with the following format:
   //
   //~~~~~~~~~~~~~~~~
   //[Dataset].[Array Name].[ID type].GridOrder:  [Grid1],[Grid2],...
   //~~~~~~~~~~~~~~~~
   //
   //where the name of each grid is given as "VARY_VARX". If no variable defining the order is found,
   //the grids will be tried in the order they were found in the file containing the grids for this
   //telescope.
   //
   // The KVIdentificationResult is first Clear()ed; then it is filled with IDtype = GetType()
   // of this identification telescope, IDattempted = true, and the results of the identification
   // procedure.
 
   idr->Clear();
   idr->IDattempted = true;
   idr->SetIDType(GetType());
 
   KVIDGraph* grid;
   TIter it(GetListOfIDGrids());
   while ((grid = (KVIDGraph*)it())) { //loop over grids in order given by [Dataset].[Array Name].[ID type].GridOrder:
      Double_t de, e;
      auto fired = GetIDGridCoords(e, de, grid, x, y);// true if all signals required by e & de fired/present in event
      idr->SetGridName(grid->GetName());
      if (fired && grid->IsIdentifiable(e, de, &idr->Rejecting_Cut)) {
         grid->Identify(e, de, idr);
         if (idr->IDOK) break; // stop on first successful identification
      }
      else {
         // either all signals required to calculate coordinates in grid are not present/fired,
         // or particle rejected by cut in grid, in which case idr->Rejecting_Cut contains its name.
         idr->IDOK = kFALSE;
         idr->IDquality = KVIDZAGrid::kICODE8;
      }
   }
   idr->IDcode = GetIDCode();
 
   return kTRUE;
}
 
 
 
 
 
void KVIDTelescope::SetIDGrid(KVIDGraph* grid)
{
   // Add an identification grid to the list of grids used by this telescope.
   //
   // If the grid's VARX and VARY parameters are set and contain the names of valid
   // detector signals (see formatting rules below) they will be used by
   // GetIDGridXCoord() and GetIDGridYCoord() to return the coordinates
   // needed to perform particle identification using the grid.
   //
   // The name of the grid is set to "VARY_VARX" (just the signal names, not the detector
   // label part - see below). This value will be stored in the
   // KVIdentificationResult corresponding to an attempted identification of a
   // KVReconstructedNucleus by this grid.
   //
   // VARX/VARY Formatting
   //
   // To be valid, grid VARX/Y parameters should be set as follows:
   //
   //~~~~~~~~~~~~~~~~~~
   //      [signal name]
   //  or  [det_label]::[signal name]
   //~~~~~~~~~~~~~~~~~~
   //
   // where
   //
   //~~~~~~~~~~~~~~~~~~
   //    [det_label] (optional)  = detector label i.e. string returned by KVDetector::GetLabel()
   //                    method for detector. By default, VARX is assumed to be the Eres detector
   //                    or last detector and VARY the DE detector or first detector
   //    [signal_name] = name of a signal defined for the detector, possibly depending
   //                    on availability of calibration
   //
   //    To see all available signals for a detector, use
   //
   //         KVDetector::GetListOfDetectorSignals()
   //~~~~~~~~~~~~~~~~~~
 
   if (grid) {
      fIDGrids.Add(grid);
      KVString det_labels_x, det_labels_y;
      KVDetectorSignal* xx = GetSignalFromGridVar(grid->GetVarX(), "X", det_labels_x);
      KVDetectorSignal* yy = GetSignalFromGridVar(grid->GetVarY(), "Y", det_labels_y);
      GraphCoords gc;
      gc.fVarX = xx;
      gc.fDetLabelsX = det_labels_x;
      gc.fVarY = yy;
      gc.fDetLabelsY = det_labels_y;
      fGraphCoords[grid] = gc;
      TString grid_name;
      if (xx && yy) {
         grid_name.Form("%s_%s", yy->GetName(), xx->GetName());
         grid->SetName(grid_name);
      }
   }
}
 
 
 
 
KVDetectorSignal* KVIDTelescope::GetSignalFromGridVar(const KVString& var, const KVString& axe, KVString& det_labels)
{
   // Deduce & return pointer to detector signal from grid VARX/VARY parameter
   //
   // To be valid, grid VARX/Y parameters should be set using mathematical expressions
   // which use the following references to detector signals for the telescope:
   //
   //~~~~~~~~~~~~~~~~~~
   //      [signal name]
   //  OR  [det_label]::[signal name]
   //~~~~~~~~~~~~~~~~~~
   //
   // where
   //
   //~~~~~~~~~~~~~~~~~~
   //    [signal_name] = name of a signal defined for 1 of the detectors of the telescope
   //    [det_label]   = optional detector label i.e. string returned by
   //                    KVDetector::GetLabel() method for detector
   //~~~~~~~~~~~~~~~~~~
   //
   // If `[det_label]` is not given, we assume for `VARX` the last (E) detector,
   // while for `VARY` we assume the first (dE) detector. If this telescope has only
   // one detector, we use it for both variables.
   //
   //    To see all available signals for a detector, use
   //
   //~~~~~~~~~~~~~~~~~~
   //         KVDetector::GetListOfDetectorSignals()
   //~~~~~~~~~~~~~~~~~~
   //
   // #### Example ####
   // Imagine a telescope which combines 2 detectors, with labels SI and CSI (in that order, i.e. SI is the dE detector,
   // CSI is the residual energy detector). The following cases are valid:
   //
   //~~~~~
   // VARX = Energy                              : use 'Energy' signal of CSI detector (default for VARX)
   // VARY = ADC                                 : use 'ADC' signal of SI detector (default for VARY)
   // VARY = CSI::Energy                         : use 'Energy' signal of CSI detector (overrides default for VARY, SI)
   // VARX = (Q3-Q3Fast)/(0.8*Q3)                : uses 'Q3' and 'Q3Fast' signals of CSI detector (default for VARX)
   // VARY = 1.5*SI::ADC - CSI::Q3Fast/CSI::Q3   : combination of signals from both detectors
   //~~~~~
   //
   // The 'det_labels' string will be filled with a comma-separated list of the labels of each detector used
   // in the expressions.
   //
   // \note if any signals are not defined, they will be evaluated as zero
 
   if (var == "") {
      Warning("GetSignalFromGridVar",
              "No VAR%s defined for grid for telescope %s. KVIDTelescope-derived class handling identification must override GetIDMapX/GetIDMapY",
              axe.Data(), GetName());
      return nullptr;
   }
 
   KVString dum = var;
   KVDetector* det = nullptr;
   KVDetectorSignal* ds(nullptr);
   KVString sig_type;
 
   // check if VARX/Y is a mathematical expression
   Bool_t is_expression = var.GetNValues("+-*/()") > 1;
   // examine each term in expression (this will split the elements in a mathematical expression,
   // or just take the whole expression if no math operators are present)
   var.Begin("+-*/()");
   bool explicit_det_ref = false;
   bool multidetexpr = false;
   KVString det_label;
   while (!var.End()) {
      KVString t = var.Next();
      // check if we have an explicit reference to a detector
      if (t.Contains("::")) {
         t.Begin("::");
         det_label = t.Next();
         auto _det = (KVDetector*)GetDetectors()->FindObjectByLabel(det_label);
         if (_det) {
            explicit_det_ref = true;
            sig_type = t.Next();
            // check signal is defined for detector
            if (_det->GetDetectorSignal(sig_type)) {
               if (det_labels.Length()) {
                  if (!det_labels.Contains(det_label)) det_labels += Form(",%s", det_label.Data());
               }
               else
                  det_labels = det_label;
            }
            else {
               Warning("GetSignalFromGridVar", "signal '%s' not found in det '%s' -> replaced by '0'", sig_type.Data(), _det->GetName());
               dum.ReplaceAll(Form("%s::%s", _det->GetLabel(), sig_type.Data()), "0");
            }
         }
         if (det && (_det != det)) multidetexpr = true; // more than 1 detector is referenced
         det = _det;
      }
   }
   // treat all cases with explicit detector references
   if (explicit_det_ref) {
      if (!multidetexpr) {
         // only 1 detector is directly referenced
         if (!is_expression) {
            // it is not a mathematical expression: this is just the name of a detector signal
            return det->GetDetectorSignal(sig_type);
         }
         else {
            // math expression with explicit reference to signal(s) of 1 detector
            sig_type = var;
            // remove all explicit references to detector: 'det' pointer has been set
            sig_type.ReplaceAll(Form("%s::", det_label.Data()), "");
            // expression will be handled below
         }
      }
      else {
         // use specific constructor for explicit detector references
         ds = new KVDetectorSignalExpression(var, dum, GetDetectors());
         if (!ds->IsValid()) {
            delete ds;
            ds = nullptr;
            Error("GetSignalFromGridVar",
                  "Problem initialising ID-grid %s coordinate for telescope %s."
                  " Check definition of VAR%s for grid (=%s)",
                  axe.Data(), GetName(), axe.Data(), var.Data());
            return ds;
         }
         fMultiDetExpressions.Add(ds);
         return ds;
      }
   }
   else {
      // no explicit reference to detectors - by default, use detector (1) for Y axis,
      // and the last detector for X axis
      if (axe == "Y" || GetSize() == 1) det = GetDetector(1);
      else det = GetDetector(GetSize());
      sig_type = var;
      det_labels = det->GetLabel();
   }
   ds = det->GetDetectorSignal(sig_type);
   if (!ds) {
      // sig_type is not the name of a known signal: assume it is an expression using known signal names
      if (!det->AddDetectorSignalExpression(sig_type, sig_type)) {
         Error("GetSignalFromGridVar",
               "Problem initialising ID-grid %s coordinate for telescope %s. Request for unknown signal %s for detector %s. Check definition of VAR%s for grid (=%s)",
               axe.Data(), GetName(), sig_type.Data(), det->GetName(), axe.Data(), var.Data());
         ds = nullptr;
      }
      else
         ds = det->GetDetectorSignal(sig_type);
   }
   return ds;
}
 
 
 
 
KVIDGrid* KVIDTelescope::newGrid(bool onlyZ)
{
   TClass* cl = TClass::GetClass(GetDefaultIDGridClass());
   if (!cl || !cl->InheritsFrom("KVIDZAGrid")) cl = TClass::GetClass("KVIDZAGrid");
   KVIDGrid* idgrid = (KVIDGrid*)cl->New();
 
   idgrid->AddIDTelescope(this);
   idgrid->SetOnlyZId(onlyZ);
   idgrid->SetRunList("1-10000");
 
   KVIDCutLine* B_line = (KVIDCutLine*)idgrid->Add("OK", "KVIDCutLine");
   Int_t npoi_bragg = 0;
   B_line->SetName("Bragg_line");
   B_line->SetAcceptedDirection("right");
 
   return idgrid;
}
 
 
 
 
void KVIDTelescope::addLineToGrid(KVIDGrid* idgrid, int zz, int aa, int npoints, std::optional<Double_t> einc_max)
{
 
   //loop over energy
   //first find :
   //  ****E1 = energy at which particle passes 1st detector and starts to enter in the 2nd one****
   //      E2 = energy at which particle passes the 2nd detector
   //then perform npoints calculations between these two energies and use these
   //to construct a KVIDZALine
   //
   // if einc_max is given, limit incident energy of particles
 
   double xfactor = 1.;
 
   KVNucleus part;
 
   KVDetector* det_de = GetDetector(1);
   KVDetector* det_eres = GetDetector(2);
 
   Double_t SeuilE = 0.1;
 
   part.SetZ(zz);
   part.SetA(aa);
 
   Double_t E1, E2;
   //find E1
   //go from SeuilE MeV to det_de->GetEIncOfMaxDeltaE(part.GetZ(),part.GetA()))
   Double_t E1min = SeuilE, E1max = det_de->GetEIncOfMaxDeltaE(zz, aa);
   E1 = (E1min + E1max) / 2.;
 
   while ((E1max - E1min) > SeuilE) {
 
      part.SetEnergy(E1);
      det_de->Clear();
      det_eres->Clear();
 
      det_de->DetectParticle(&part);
      det_eres->DetectParticle(&part);
      if (det_eres->GetEnergy() > SeuilE) {
         //particle got through - decrease energy
         E1max = E1;
         E1 = (E1max + E1min) / 2.;
      }
      else {
         //particle stopped - increase energy
         E1min = E1;
         E1 = (E1max + E1min) / 2.;
      }
   }
 
   //add point to Bragg line
   Double_t dE_B = det_de->GetMaxDeltaE(zz, aa);
   Double_t E_B = det_de->GetEIncOfMaxDeltaE(zz, aa);
   Double_t Eres_B = det_de->GetERes(zz, aa, E_B);
 
   KVIDCutLine* B_line = (KVIDCutLine*)idgrid->GetCut("Bragg_line");
   if (B_line) B_line->SetPoint(B_line->GetN(), Eres_B, dE_B);
 
   //find E2
   //go from E1 MeV to maximum value where the energy loss formula is valid
   Double_t E2min = E1, E2max = det_eres->GetEmaxValid(part.GetZ(), part.GetA());
   E2 = (E2min + E2max) / 2.;
 
   while ((E2max - E2min > SeuilE)) {
 
      part.SetEnergy(E2);
      det_de->Clear();
      det_eres->Clear();
 
      det_de->DetectParticle(&part);
      det_eres->DetectParticle(&part);
      if (part.GetEnergy() > SeuilE) {
         //particle got through - decrease energy
         E2max = E2;
         E2 = (E2max + E2min) / 2.;
      }
      else {
         //particle stopped - increase energy
         E2min = E2;
         E2 = (E2max + E2min) / 2.;
      }
   }
   E2 *= xfactor;
   if ((!strcmp(det_eres->GetType(), "CSI")) && (E2 > 5000)) E2 = 5000;
   if(einc_max) E2 = std::min(E2, einc_max.value());
   //                printf("z=%d a=%d E1=%lf E2=%lf\n",zz,aa,E1,E2);
   KVIDZALine* line = (KVIDZALine*)idgrid->Add("ID", "KVIDZALine");
   if (TMath::Even(zz)) line->SetLineColor(4);
   line->SetZ(zz);
   line->SetA(aa);
 
   Double_t logE1 = TMath::Log(E1);
   Double_t logE2 = TMath::Log(E2);
   Double_t dLog = (logE2 - logE1) / (npoints - 1.);
 
   for (Int_t i = 0; i < npoints; i++) {
      //              Double_t E = E1 + i*(E2-E1)/(npoints-1.);
      Double_t E = TMath::Exp(logE1 + i * dLog);
 
      Double_t Eres = 0.;
      Int_t niter = 0;
      while (Eres < SeuilE && niter <= 20) {
         det_de->Clear();
         det_eres->Clear();
 
         part.SetEnergy(E);
 
         det_de->DetectParticle(&part);
         det_eres->DetectParticle(&part);
 
         Eres = det_eres->GetEnergy();
         E += SeuilE;
         niter += 1;
      }
      if (!(niter > 20)) {
         Double_t dE = det_de->GetEnergy();
         Double_t gEres, gdE;
         line->GetPoint(i - 1, gEres, gdE);
         line->SetPoint(i, Eres, dE);
 
      }
   }
   //printf("sort de boucle points");
 
 
}
 
 
 
 
KVString KVIDTelescope::GetDetectorLabelsForGridCoord(const KVString& axis) const
{
   // Returns a comma-separated list of the labels of the detectors used to determine
   // the "x" or "y" coordinates of the identification grid(s)
   //
   // If there is more than 1 grid and the list is not the same for all grids,
   // prints a warning message.
   //
   // \param[in] axis name of grid axis i.e. "x", "X", "y" or "Y" (case insensitive)
 
   KVString _axis = axis;
   _axis.ToUpper();
 
   if (_axis != "X" && _axis != "Y") {
      Error("GetDetectorLabelsForGridCoord", "Called with illegal axis name '%s'", axis.Data());
      return "";
   }
   KVString lab_list;
 
   for (auto& __gc : fGraphCoords) {
      KVString labs;
      if (_axis == "X") labs = __gc.second.fDetLabelsX;
      else labs = __gc.second.fDetLabelsY;
      if (lab_list.Length() && lab_list != labs) {
         Error("GetDetectorLabelsForGridCoord", "Grids for telescope %s use different detector types for %s-coordinate: "
               "%s and %s", GetName(), _axis.Data(), lab_list.Data(), labs.Data());
         return "";
      }
      lab_list = labs;
   }
   return lab_list;
}
 
 
 
 
 
KVIDGraph* KVIDTelescope::GetIDGrid()
{
   //Return the first in the list of identification grids used by this telescope
   //(this is for backwards compatibility with ID telescopes which had only one grid).
   return (KVIDGraph*)GetListOfIDGrids()->First();
}
 
 
 
 
 
KVIDGraph* KVIDTelescope::GetIDGrid(Int_t index)
{
   //Return pointer to grid using position in list. First grid has index = 1.
   if (index < 1) {
      Error("GetIDGrid(int)", "Index must be >=1!");
      return nullptr;
   }
   return (KVIDGraph*)GetListOfIDGrids()->At(index - 1);
}
 
 
 
 
 
KVIDGraph* KVIDTelescope::GetIDGrid(const Char_t* label)
{
   //Return pointer to grid using "label" to search in list of grids associated
   //to this telescope.
   return (KVIDGraph*)GetListOfIDGrids()->FindObjectByLabel(label);
}
 
 
 
 
 
Double_t KVIDTelescope::GetIDMapX(Option_t*)
{
   AbstractMethod("GetIDMapX");
   return -1.;
}
 
 
 
 
Double_t KVIDTelescope::GetPedestalX(Option_t*)
{
   // Returns the pedestal associated with the 2nd detector of the telescope,
   // optionally depending on the given option string.
   // By default this returns 0, and should be overridden in specific implementations.
 
   return 0.;
}
 
 
 
 
Double_t KVIDTelescope::GetPedestalY(Option_t*)
{
   // Returns the pedestal associated with the 1st detector of the telescope,
   // optionally depending on the given option string.
   // By default this returns 0, and should be overridden in specific implementations.
 
   return 0.;
}
 
 
 
 
Double_t KVIDTelescope::GetIDGridXCoord(KVIDGraph* g, bool& fired) const
{
   // \return value of X coordinate to be used with the given ID grid
   //
   // This corresponds to whatever was given as parameter "VARX" for the grid
   //
   // \param[out] fired set to true or false depending on whether detector signal(s) used for X were fired (read from last raw event)
 
   KVDetectorSignal* ds = fGraphCoords[g].fVarX;
   if (ds) {
      fired = ds->IsFired();
      return ds->GetValue();
   }
   fired = false;
   return -1;
}
 
 
 
 
Double_t KVIDTelescope::GetIDGridYCoord(KVIDGraph* g, bool& fired) const
{
   // \return value of Y coordinate to be used with the given ID grid
   //
   // This corresponds to whatever was given as parameter "VARY" for the grid
   //
   // \param[out] fired set to true or false depending on whether detector signal(s) used for Y were fired (read from last raw event)
 
   KVDetectorSignal* ds = fGraphCoords[g].fVarY;
   if (ds) {
      fired = ds->IsFired();
      return ds->GetValue();
   }
   fired = false;
   return -1;
}
 
 
 
 
 
Double_t KVIDTelescope::GetIDMapY(Option_t*)
{
   AbstractMethod("GetIDMapY");
   return -1.;
}
 
 
 
 
 
void KVIDTelescope::RemoveGrids()
{
   //Remove all identification grids for this ID telescope
   //Grids are not deleted as this is handled by gIDGridManager
   fIDGrids.Clear();
   fGraphCoords.clear();
}
 
 
 
 
 
KVIDTelescope* KVIDTelescope::MakeIDTelescope(const Char_t* uri)
{
   //Static function which will create an instance of the KVIDTelescope-derived class
   //corresponding to 'name'
   //These are defined as 'Plugin' objects in the file $KVROOT/KVFiles/.kvrootrc :
   //~~~~~~~
   // # The KVMultiDetArray::GetIDTelescopes(KVDetector*de, KVDetector*e) method uses these plugins to
   // # create KVIDTelescope instances adapted to the specific array geometry and detector types.
   // # For each pair of detectors we look for a plugin with one of the following names:
   // #    [name_of_dataset].de_detector_type[de detector thickness]-e_detector_type[de detector thickness]
   // # Each characteristic in [] brackets may or may not be present in the name; first we test for names
   // # with these characteristics, then all combinations where one or other of the characteristics is not present.
   // # In addition, we first test all combinations which begin with [name_of_dataset].
   // # The first plugin found in this order will be used.
   // # In addition, if for one of the two detectors there is a plugin called
   // #    [name_of_dataset].de_detector_type[de detector thickness]
   // #    [name_of_dataset].e_detector_type[e detector thickness]
   // # then we add also an instance of this 1-detector identification telescope.
   //~~~~~~~
 
   TPluginHandler* ph;
   //check and load plugin library
   if (!(ph = LoadPlugin("KVIDTelescope", uri)))
      return 0;
 
   //execute constructor/macro for identification telescope
   KVIDTelescope* mda = (KVIDTelescope*) ph->ExecPlugin(0);
   if (mda) {
      //set label of telescope with URI used to find plugin (minus dataset name)
      mda->SetLabelFromURI(uri);
   }
 
   return mda;
}
 
 
 
 
 
void KVIDTelescope::SetLabelFromURI(const Char_t* uri)
{
   //PRIVATE METHOD
   //Sets label of telescope based on URI of plugin describing child class for this telescope
 
   TString _uri(uri);
   if (gDataSet && _uri.BeginsWith(gDataSet->GetName())) _uri.Remove(0, strlen(gDataSet->GetName()) + 1);
   SetLabel(_uri.Data());
}
 
 
 
 
 
Bool_t KVIDTelescope::SetIdentificationParameters(const KVMultiDetArray* multidet)
{
   // Initialise the identification parameters (grids, etc.) of ALL identification telescopes of this
   // kind (label) in the multidetector array. Therefore this method need only be called once, and not
   // called for each telescope. The kind/label (returned by GetLabel) of the telescope corresponds
   // to the URI used to find the plugin class in $KVROOT/KVFiles/.kvrootrc.
   //
   // By default, this method first looks for a file containing the list of all files containing the
   // identification grids. The default name of this file is:
   //
   //~~~
   // [label].IdentificationFiles.dat
   //~~~
   //
   // or any name given by the environment variable
   //
   //~~~
   // [dataset].IdentificationParameterList.[label]:       [filename]
   //~~~
   //
   // If neither file is found, we look for a file with name:
   //
   //~~~
   // IDGrids.[label].dat
   //~~~
   //
   // or any name given by the environment variable
   //
   //~~~
   // [dataset].IdentificationParameterFile.[label]:       [filename]
   //~~~
   //
   // which is assumed to contain all identification grids for this type of telescope.
 
   Bool_t ok = ReadIdentificationParameterFiles(multidet);
   if(!ok) ok = LoadIdentificationParameters(multidet);
 
   if (!ok)
      Warning("SetIdentificationParameters", "No files found. See doc for this method.");
 
   return ok;
}
 
 
 
 
 
Bool_t KVIDTelescope::ReadIdentificationParameterFiles(const KVMultiDetArray* multidet)
{
   // In the case where the identification grids are stored in several files, this method looks
   // for a file with name
   //
   //~~~
   // [label].IdentificationFiles.dat
   //~~~
   //
   // or any name given by the environment variable
   //
   //~~~
   // [dataset].IdentificationParameterList.[label]:       [filename]
   //~~~
   //
   // and if found reads from it the list of files containing the identification grids.
   //
   // \returns kTRUE if file found and read, kFALSE if no file found
 
   TString filename = gDataSet->GetDataSetEnv(Form("IdentificationParameterList.%s", GetLabel()),
                                              Form("%s.IdentificationFiles.dat", GetLabel()));
 
   if(gDataSet->FindDataSetFile(filename))
   {
      KVFileReader fr;
      auto path = gDataSet->GetFullPathToDataSetFile(filename);
      Info("ReadIdentificationParameterFiles", "Using file %s", path.Data());
      fr.OpenFileToRead(path);
      Bool_t ok = kTRUE;
      while (fr.IsOK()) {
         fr.ReadLine(0);
 
         if (fr.GetCurrentLine() != "")
            ok &= LoadIdentificationParameters(multidet, fr.GetCurrentLine().Data());
      }
 
      fr.CloseFile();
      return ok;
   }
 
   return kFALSE;
}
 
 
 
 
 
Bool_t KVIDTelescope::LoadIdentificationParameters(const KVMultiDetArray* multidet, const TString& _filename)
{
   // \param[in] _filename optional name of dataset file containing grids for this telescope type
   // \returns kTRUE if a file is found and read in
   //
   // If no value for _filename is given, this method looks for a dataset file named
   //
   //~~~
   // IDGrids.[label].dat
   //~~~
   //
   // or any name given by the environment variable
   //
   //~~~
   // [dataset].IdentificationParameterFile.[label]:       [filename]
   //~~~
   //
   // which is assumed to contain all identification grids for this type of telescope.
 
   TString filename;
   if(_filename.IsNull())
   {
      filename = gDataSet->GetDataSetEnv(Form("IdentificationParameterFile.%s", GetLabel()),
                                         Form("IDGrids.%s.dat",GetLabel()));
   }
   else
      filename = _filename;
 
   if(gDataSet->FindDataSetFile(filename))
   {
      auto path = gDataSet->GetFullPathToDataSetFile(filename);
      //Read grids from file
      Info("LoadIdentificationParameters", "Using file %s", path.Data());
      multidet->ReadGridsFromAsciiFile(path);
      return kTRUE;
   }
   else if (!_filename.IsNull())
   {
      // if the method was called with a definite filename to look for, we assume that this
      // file is supposed to exist, therefore if it is not found, there is a problem!
      Error("LoadIdentificationParameters",
            "File %s not found. Should be in %s",
            _filename.Data(), gDataSet->GetDataSetDir());
      return kFALSE;
   }
   return kFALSE;
}
 
 
 
 
 
void KVIDTelescope::RemoveIdentificationParameters()
{
   //Remove identification parameters from telescope in such a way that they
   //can subsequently be reset e.g. with a new version.
   //This is used by KVMultiDetArray::UpdateIdentifications.
   //Child classes with specific SetIdentificationParameters methods should
   //also redefine this method in order to remove (destroy) cleanly the objects
   //created in SetIdentificationParameters.
   //
   //This default method takes the list of grids associated to the telescope,
   //and for each one: 1) checks if it is still in the gIDGridManager's list
   //2) if yes, delete the grid and remove it from gIDGridManager
 
   TIter next_grid(GetListOfIDGrids());
   KVIDGrid* grid;
   while ((grid = (KVIDGrid*)next_grid())) {
 
      if (gIDGridManager->GetGrids()->FindObject(grid)) {   //this only works if KVIDTelescope uses TObject:IsEqual method (i.e. compares pointers)
 
         gIDGridManager->DeleteGrid(grid);
 
      }
   }
   //clear list of grids
   fIDGrids.Clear();
   SetHasMassID(kFALSE);
}
 
 
 
//void KVIDTelescope::CalculateParticleEnergy(KVReconstructedNucleus* nuc)
//{
//   // The energy of each particle is calculated as follows:
//   //
//   //      E = dE_1 + dE_2 + ... + dE_N
//   //
//   // dE_1, dE_2, ... = energy losses measured in each detector through which
//   //                          the particle has passed (or stopped, in the case of dE_N).
//   //                         These energy losses are corrected for (Z,A)-dependent effects
//   //                          such as pulse-heigth defect in silicon detectors, losses in
//   //                          windows of gas detectors, etc.
//   //
//   // Whenever possible, the energy loss for fired detectors which are uncalibrated
//   // or not functioning is calculated. In this case the status returned by GetCalibStatus()
//   // will be KVIDTelescope::kCalibStatus_Calculated.
//   // If none of the detectors is calibrated, the particle's energy cannot be calculated &
//   // the status will be KVIDTelescope::kCalibStatus_NoCalibrations.
//   // Otherwise, the status code will be KVIDTelescope::kCalibStatus_OK.
 
//   //status code
//   fCalibStatus = kCalibStatus_NoCalibrations;
 
//   UInt_t z = nuc->GetZ();
//   //uncharged particles
//   if (z == 0) return;
 
//   KVDetector* d1 = GetDetector(1);
//   KVDetector* d2 = (GetSize() > 1 ? GetDetector(2) : 0);
//   Bool_t d1_cal = d1->IsCalibrated();
//   Bool_t d2_cal = (d2 ? d2->IsCalibrated() : kFALSE);
 
//   //no calibrations
//   if (!d1_cal && !d2)
//      return;
//   if ((d1 && d2) && !d1_cal && !d2_cal)
//      return;
 
//   //status code
//   fCalibStatus = kCalibStatus_OK;
 
//   UInt_t a = nuc->GetA();
 
//   // particles stopped in first member of telescope
//   if (nuc->GetStatus() == 3) {
//      if (d1_cal) {
//         nuc->SetEnergy(d1->GetCorrectedEnergy(nuc, -1, kFALSE));  //N.B.: transmission=kFALSE because particle stop in d1
//      }
//      return;
//   }
 
//   Double_t e1, e2, einc;
//   e1 = e2 = einc = 0.0;
 
//   if (!d1_cal) {//1st detector not calibrated - calculate from residual energy in 2nd detector
 
//      //second detector must exist and have all acquisition parameters fired with above-pedestal value
//      if (d2 && d2->Fired("Pall")) e2 = d2->GetCorrectedEnergy(nuc, -1, kFALSE); //N.B.: transmission=kFALSE because particle stop in d2
//      if (e2 <= 0.0) {
//         // zero energy loss in 2nd detector ? can't do anything...
//         fCalibStatus = kCalibStatus_NoCalibrations;
//         return;
//      }
//      //calculate & set energy loss in dE detector
//      //N.B. using e2 for the residual energy after detector 1 means
//      //that we are assuming the particle stops in detector 2.
//      //if this is not true, we will underestimate the energy of the particle.
//      e1 = d1->GetDeltaEFromERes(z, a, e2);
//      if (e1 < 0.0) e1 = 0.0;
//      else {
//         d1->SetEnergyLoss(e1);
//         d1->SetEResAfterDetector(e2);
//         e1 = d1->GetCorrectedEnergy(nuc);
//         //status code
//         fCalibStatus = kCalibStatus_Calculated;
//      }
//   }
//   else {  //1st detector is calibrated too: get corrected energy loss
 
//      e1 = d1->GetCorrectedEnergy(nuc);
 
//   }
 
//   if (d2 && !d2_cal) {//2nd detector not calibrated - calculate from energy loss in 1st detector
 
//      e1 = d1->GetCorrectedEnergy(nuc);
//      if (e1 <= 0.0) {
//         // zero energy loss in 1st detector ? can't do anything...
//         fCalibStatus = kCalibStatus_NoCalibrations;
//         return;
//      }
//      //calculate & set energy loss in 2nd detector
//      e2 = d1->GetEResFromDeltaE(z, a);
//      if (e2 < 0.0) e2 = 0.0;
//      else {
//         e2 = d2->GetDeltaE(z, a, e2);
//         d2->SetEnergyLoss(e2);
//         e2 = d2->GetCorrectedEnergy(nuc);
//         //status code
//         fCalibStatus = kCalibStatus_Calculated;
//      }
//   }
//   else if (d2) {   //2nd detector is calibrated too: get corrected energy loss
 
//      e2 = d2->GetCorrectedEnergy(nuc, -1, kFALSE);//N.B.: transmission=kFALSE because particle assumed to stop in d2
//      // recalculate corrected energy in first stage using info on Eres
//      d1->SetEResAfterDetector(e2);
//      e1 = d1->GetCorrectedEnergy(nuc);
//   }
 
//   //incident energy of particle (before 1st member of telescope)
//   einc = e1 + e2;
 
//   Double_t coherence_tolerance = gEnv->GetValue("KVIDTelescope.CoherencyTolerance", 1.05);
//   if (coherence_tolerance < 1) coherence_tolerance += 1.00;
 
//   //Now we have to work our way up the list of detectors from which the particle was
//   //reconstructed. For each fired & calibrated detector which is only associated with
//   //one particle in the events, we add the corrected measured energy loss
//   //to the particle. For uncalibrated, unfired detectors and detectors through which
//   //more than one particle has passed, we calculate the corrected energy loss and add it
//   //to the particle.
//   int ndets = nuc->GetNumDet();
//   if (ndets > (int)GetSize()) { //particle passed through other detectors before this idtelesocpe
//      //look at detectors not in this id telescope
//      int idet = GetSize();//next detector after delta-e member of IDTelescope (stopping detector = 0)
//      while (idet < ndets) {
 
//         KVDetector* det = nuc->GetDetector(idet);
//         if (det->Fired() && det->IsCalibrated() && det->GetNHits() == 1) {
//            Double_t dE = det->GetEnergy();
//            //in order to check if particle was really the only one to
//            //hit each detector, we calculate the particle's energy loss
//            //from its residual energy. if the measured energy loss is
//            //significantly larger, there may be a second particle.
//            e1 = det->GetDeltaEFromERes(z, a, einc);
//            if (e1 < 0.0) e1 = 0.0;
//            det->SetEResAfterDetector(einc);
//            dE = det->GetCorrectedEnergy(nuc);
//            einc += dE;
//         }
//         else {
//            // Uncalibrated/unfired/multihit detector. Calculate energy loss.
//            //calculate energy of particle before detector from energy after detector
//            e1 = det->GetDeltaEFromERes(z, a, einc);
//            if (e1 < 0.0) e1 = 0.0;
//            if (det->GetNHits() > 1) {
//               //Info("CalculateParticleEnergy",
//               //    "Detector %s was hit by %d particles. Calculated energy loss for particle %f MeV",
//               //    det->GetName(), det->GetNHits(), e1);
//               if (!(det->Fired() && det->IsCalibrated())) {
//                  det->SetEnergyLoss(e1 + det->GetEnergy());// sum up calculated energy losses in uncalibrated detector
//               }
//               //status code
//               fCalibStatus = kCalibStatus_Multihit;
//            }
//            else if (!det->Fired() || !det->IsCalibrated()) {
//               //Info("CalculateParticleEnergy",
//               //    "Detector %s uncalibrated/not fired. Calculated energy loss for particle %f MeV",
//               //    det->GetName(), e1);
//               det->SetEnergyLoss(e1);
//               //status code
//               fCalibStatus = kCalibStatus_Calculated;
//            }
//            det->SetEResAfterDetector(einc);
//            e1 = det->GetCorrectedEnergy(nuc, e1);
//            einc += e1;
//         }
//         idet++;
//      }
//   }
//   //einc is now the energy of the particle before crossing the first detector
//   nuc->SetEnergy(einc);
//}
 
 
 
 
const Char_t* KVIDTelescope::GetDefaultIDGridClass()
{
   // Returns name of default ID grid class for this ID telescope.
   // This is defined in a .kvrootrc configuration file by one of the following:
   // KVIDTelescope.DefaultGrid:
   // KVIDTelescope.DefaultGrid.[type]:
   // where [type] is the type of this identification telescope (which is given
   // by the character string returned by method GetLabel()... sorry :( )
   // If no default grid is defined for the specific type of this telescope,
   // the default defined by KVIDTelescope.DefaultGrid is used.
 
   TString specific;
   specific.Form("KVIDTelescope.DefaultGrid.%s", GetLabel());
   return gEnv->GetValue(specific.Data(), gEnv->GetValue("KVIDTelescope.DefaultGrid", "KVIDGraph"));
}
 
 
 
 
KVIDGrid* KVIDTelescope::CalculateDeltaE_EGrid(const KVNameValueList& AperZ, Int_t npoints, Double_t xfactor)
{
   //Create a dE-E grid (energy loss in detector 1 versus residual energy in detector 2) for a given list of isotopes
   // - AperZ   : list of A for each Z. (ex: "1=1-3,2=3-4 5,3=6-8,4=7 9-12"...)
   // - npoints : number of points in each generated line
   // - xfactor : scales the detector 2 thickness to prolongate the lines
 
   if (GetSize() <= 1) return 0;
   if (!GetDetector(1) || !GetDetector(2))      return 0;
   double thickness = GetDetector(2)->GetThickness();
   GetDetector(2)->SetThickness(thickness * xfactor);
 
   Info("CalculateDeltaE_EGrid", "called with KVNameValueList");
 
   KVIDGrid* idgrid = newGrid(0);
 
   for (auto& par : AperZ) {
      KVString tmp = par.GetName();
      int zz = tmp.Atoi();
      KVNumberList alist = par.GetString();
      for (auto aa : alist) {
         addLineToGrid(idgrid, zz, aa, npoints);
      }
   }
 
   GetDetector(2)->SetThickness(thickness);
 
   return idgrid;
 
}
 
 
 
 
KVIDGrid* KVIDTelescope::CalculateDeltaE_EGrid(const KVNumberList& Zrange, Int_t deltaA, Int_t npoints, Double_t lifetime, UChar_t massformula, Double_t xfactor)
{
   //Create a dE-E grid (energy loss in detector 1 versus residual energy in detector 2) for a given list of isotopes
   // - Zrange  : list of element for which we create lines
   // - deltaA  : number of isotopes generated for each Z around massformula (ex: deltaA=1, Aref-1 Aref Aref+1)
   // - npoints : number of points in each generated line
   // - lifetime: remove isotopes with lifetime lower than this value
   // - xfactor : scales the detector 2 thickness to prolongate the lines
 
   if (GetSize() <= 1) return 0;
   if (!GetDetector(1) || !GetDetector(2))      return 0;
   double thickness = GetDetector(2)->GetThickness();
   GetDetector(2)->SetThickness(thickness * xfactor);
 
   Info("CalculateDeltaE_EGrid", "called with KVNumberList");
 
   KVIDGrid* idgrid = newGrid(0);
   KVNucleus part;
 
   Zrange.Begin();
   while (!Zrange.End()) {
      Int_t zz = Zrange.Next();
      part.SetZ(zz, massformula);
      Int_t aref = part.GetA();
      for (Int_t aa = aref - deltaA; aa <= aref + deltaA; aa += 1) {
         part.SetA(aa);
         if (part.IsKnown() && (part.GetLifeTime() > lifetime)) {
            addLineToGrid(idgrid, zz, aa, npoints);
 
         }
      }
   }
   GetDetector(2)->SetThickness(thickness);
   return idgrid;
}
 
 
 
 
KVIDGrid* KVIDTelescope::CalculateDeltaE_EGrid(const KVNumberList& Zrange, Int_t npoints, Double_t Einc_max, UChar_t massformula)
{
   // Generate ID grid for given list of Z values (1 line per Z), up to the given maximum incident energy
   if (GetSize() <= 1) return 0;
   if (!GetDetector(1) || !GetDetector(2))      return 0;
 
   Info("CalculateDeltaE_EGrid", "called with Einc_max=%f", Einc_max);
 
   KVIDGrid* idgrid = newGrid(0);
   KVNucleus part;
 
   Zrange.Begin();
   while (!Zrange.End()) {
      Int_t zz = Zrange.Next();
      part.SetZ(zz, massformula);
      addLineToGrid(idgrid, zz, part.GetA(), npoints, Einc_max);
   }
   return idgrid;
}
 
 
 
 
KVIDGrid* KVIDTelescope::CalculateDeltaE_EGrid(TH2* haa_zz, Bool_t Zonly, Int_t npoints)
{
   //Create a dE-E grid (energy loss in detector 1 versus residual energy in detector 2) for a given list of isotopes
   // - haa_zz  : lines will be generated for A,Z filled with 1 in this histogram
   // - Zonly   : if true, generate only one line per Z with the <A>(Z) of the histogram
   // - npoints : number of points in each generated line
 
   if (GetSize() <= 1) return 0;
   if (!GetDetector(1) || !GetDetector(2))      return 0;
   double thickness = GetDetector(2)->GetThickness();
 
   Info("CalculateDeltaE_EGrid", "called with TH2");
 
   KVIDGrid* idgrid = newGrid(0);
   KVNucleus part;
 
   for (Int_t nx = 1; nx <= haa_zz->GetNbinsX(); nx += 1) {
 
      Int_t zz = TMath::Nint(haa_zz->GetXaxis()->GetBinCenter(nx));
      KVNumberList nlA;
      Double_t sumA = 0, sum = 0;
      for (Int_t ny = 1; ny <= haa_zz->GetNbinsY(); ny += 1) {
         Double_t stat = haa_zz->GetBinContent(nx, ny);
         if (stat > 0) {
            Double_t val = haa_zz->GetYaxis()->GetBinCenter(ny);
            nlA.Add(TMath::Nint(val));
            sumA += val * stat;
            sum  += stat;
         }
      }
      sumA /= sum;
      Int_t nA = nlA.GetNValues();
      if (nA == 0) {
         Warning("CalculateDeltaE_EGrid", "no count for Z=%d", zz);
      }
      else {
         if (Zonly) {
            nlA.Clear();
            nlA.Add(TMath::Nint(sumA));
         }
         else {
            if (nA == 1) {
               Int_t aref = nlA.Last();
               nlA.Add(aref - 1);
               nlA.Add(aref + 1);
            }
         }
         part.SetZ(zz);
         nlA.Begin();
         while (!nlA.End()) {
            Int_t aa = nlA.Next();
            part.SetA(aa);
            if (part.IsKnown()) {
               addLineToGrid(idgrid, zz, aa, npoints);
            }
         }
      }
   }
   return idgrid;
}
 
 
 
 
Double_t KVIDTelescope::GetMeanDEFromID(Int_t& status, Int_t Z, Int_t A, Double_t Eres)
{
   // Returns the Y-axis value in the 2D identification map containing isotope (Z,A)
   // corresponding to either the given X-axis/Eres value or the current X-axis value given by GetIDGridXCoord()
   // If no mass information is available, just give Z.
   //
   // In this (default) implementation this means scanning the ID grids associated with
   // this telescope until we find an identification line Z or (Z,A), and then interpolating
   // the Y-coordinate for the current X-coordinate value.
   //
   // Status variable can take one of following values:
   //
   //  KVIDTelescope::kMeanDE_OK                    all OK
   //   KVIDTelescope::kMeanDE_XtooSmall      X-coordinate is smaller than smallest X-coordinate of ID line
   //   KVIDTelescope::kMeanDE_XtooLarge     X-coordinate is larger than largest X-coordinate of ID line
   //   KVIDTelescope::kMeanDE_NoIdentifie    No identifier found for Z or (Z,A)
 
   status = kMeanDE_OK;
   // loop over grids
   TIter next(GetListOfIDGrids());
   KVIDGrid* grid;
   KVIDLine* idline = 0;
   while ((grid = (KVIDGrid*)next())) {
      idline = (KVIDLine*)grid->GetIdentifier(Z, A);
      if (idline) break;
   }
   if (!idline) {
      status = kMeanDE_NoIdentifier;
      return -1.;
   }
   Double_t x, x1, y1, x2, y2; bool f;
   x = (Eres < 0 ? GetIDGridXCoord(grid,f) : Eres);
   idline->GetEndPoint(x2, y2);
   if (x > x2) {
      status = kMeanDE_XtooLarge;
      return -1;
   }
   idline->GetStartPoint(x1, y1);
   if (x < x1) {
      status = kMeanDE_XtooSmall;
      return -1.;
   }
   return idline->Eval(x);
}
 
 
 
 
 
Bool_t KVIDTelescope::CheckTheoreticalIdentificationThreshold(KVNucleus* ION, Double_t EINC)
{
   // Return kTRUE if energy of ION is > minimum incident energy required for identification
   // This theoretical limit is defined here to be the incident energy for which the
   // dE in the first detector of a dE-E telescope is maximum.
   // If EINC>0 it is assumed to be the energy of the ion just before the first detector
   // (case where ion would have to pass other detectors before reaching this telescope).
   //
   // If this is not a dE-E telescope, we return kTRUE by default.
 
   if (GetSize() < 2) return kTRUE;
 
   KVDetector* dEdet = GetDetector(1);
   Double_t emin = dEdet->GetEIncOfMaxDeltaE(ION->GetZ(), ION->GetA());
   if (EINC > 0.0) return (EINC > emin);
   return (ION->GetEnergy() > emin);
}
 
 
 
 
void KVIDTelescope::SetIdentificationStatus(KVIdentificationResult* IDR, const KVNucleus* n)
{
   // For filtering simulations
   //
   // \param IDR identification result object for this telescope
   // \param n the simulated particle currently being considered
   //
   // Set the IDR->Zident and IDR->Aident status of the identification of n.
   // In principle this depends on whether this telescope provides mass
   // identification or not, but this may depend on the particle's energy.
   //
   // In order to enable mass identification for certain telescopes without a dedicated
   // implementation (e.g. for simulating array response), put the following lines
   // in your .kvrootrc:
   //
   //   [dataset].[telescope label].MassID:   yes
   //
   // If you want to limit mass identification to certain values of Z and/or A,
   // add the following line:
   //
   //   [dataset].[telescope label].MassID.Validity:    [expression]
   //
   // where [expression] is some valid C++ boolean expression involving Z and/or A,
   // for example
   //
   //   [dataset].[telescope label].MassID.Validity:  (Z>3)&&(A<20)
   //
   // Then this expression will be tested here in order to determine particle
   // identification status
 
   KVNucleus test(IDR->Z, IDR->A);
   if (!test.IsKnown()) { // reject weird identifications as pile-up (in between lines)
      IDR->Zident = false;
      IDR->IDquality = 4;
      IDR->IDOK = false;
      return;
   }
   IDR->Zident = true;
   if (!HasMassID()) {
      IDR->Aident = false;
   }
   else {
      if (fMassIDValidity) IDR->Aident = fMassIDValidity->Test(n); // test expression for mass ID validity
      else IDR->Aident = true; // no expression set; all nuclei are identified in mass
   }
}
 
 
 
 
void KVIDTelescope::OpenIdentificationBilan(const TString& path)
{
   // Open IdentificationBilan.dat file with given path
 
   if (fgIdentificationBilan) delete fgIdentificationBilan;
   fgIdentificationBilan = new TEnv(path);
}
 
 
 
 
void KVIDTelescope::CheckIdentificationBilan(const TString& system)
{
   // Set status of ID Telescope for given system
   if (!(fgIdentificationBilan->GetValue(Form("%s.%s", system.Data(), GetName()), kTRUE))) ResetBit(kReadyForID);
}