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 "KVTelescope.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);
}
 
 
 
 
 
TGraph* KVIDTelescope::MakeIDLine(KVNucleus* nuc, Double_t Emin,
                                  Double_t Emax, Double_t Estep)
{
   //For a given nucleus, generate a TGraph representing the line in the deltaE-E
   //plane of the telescope which can be associated with nuclei of this (A,Z) with total
   //incident energies between the two limits.
   //NOTE: if there are other absorbers/detectors placed before this telescope,
   //the energy losses of the particle in these will be taken into account.
   //If the step in energy is not given, it is taken equal to 100 equal steps between min and max.
   //The TGraph should be deleted by the user after use.
 
 
   if (!Estep)
      Estep = (Emax - Emin) / 100.;
 
   Int_t nsteps = 1 + (Int_t)((Emax - Emin) / Estep);
 
   if (nsteps < 1)
      return 0;
 
   Double_t* y = new Double_t[nsteps];
   Double_t* x = new Double_t[nsteps];
   Int_t step = 0;
 
   //get list of all detectors through which particle must pass in order to reach
   //2nd member of ID Telescope
   TList* detectors =
      GetDetector(2)->GetAlignedDetectors(KVGroup::kForwards);
   //detectors->ls();
   TIter next_det(detectors);
   //cout << "nsteps =" << nsteps << endl;
 
   for (Double_t E = Emin; E <= Emax; E += Estep) {
      //Set energy of nucleus
      nuc->SetEnergy(E);
      //cout << "Einc=" << E << endl;
 
      //Calculate energy loss in each member and stock in arrays x & y
      //first member
      KVDetector* det = 0;
      x[step] = y[step] = -1;
      while ((det = (KVDetector*) next_det())) {
         //det->Print();
         Double_t eloss = det->GetELostByParticle(nuc);
         if (det == GetDetector(1))
            y[step] = eloss;
         else if (det == GetDetector(2))
            x[step] = eloss;
         Double_t E1 = nuc->GetEnergy() - eloss;
         nuc->SetEnergy(E1);
         //cout << "Eloss=" << eloss << endl;
         //cout << "Enuc=" << nuc->GetEnergy() << endl;
         if (E1 < 1.e-3) break;
      }
 
      //cout << "step = " << step << " x = " << x[step] << " y = " << y[step] << endl;
 
      //make sure that some energy is lost in each member
      //otherwise miss a step and reduce number of points in graph
      if (x[step] > 0 && y[step] > 0) {
         step++;
      }
      else {
         nsteps--;
      }
 
      //cout << "nsteps =" << nsteps << endl;
      //reset iterator ready for next loop on detectors
      next_det.Reset();
   }
   TGraph* tmp = 0;
   if (nsteps > 1)
      tmp = new TGraph(nsteps, x, y);
   delete[]x;
   delete[]y;
   return tmp;
}
 
 
 
 
 
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;
      GetIDGridCoords(e, de, grid, x, y);
      idr->SetGridName(grid->GetName());
      if (grid->IsIdentifiable(e, de, &idr->Rejecting_Cut)) {
         grid->Identify(e, de, idr);
         if (idr->IDOK) break; // stop on first successful identification
      }
      else {
         // particle rejected by cut in grid. 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)
{
 
   //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
 
   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;
   //                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) 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
 
   KVDetectorSignal* ds = fGraphCoords[g].fVarX;
   if (ds) return ds->GetValue();
   return -1;
}
 
 
 
 
Double_t KVIDTelescope::GetIDGridYCoord(KVIDGraph* g) 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
 
   KVDetectorSignal* ds = fGraphCoords[g].fVarY;
   if (ds) return ds->GetValue();
   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 looks for the file with name given by the environment variable
   //
   // [dataset name].IdentificationParameterList.[telescope label]:       [filename]
   //
   // which is assumed to contain the list of files containing the identification grids.
   //
   // If not such envionment variable is found, the method looks for another one:
   //
   // [dataset name].IdentificationParameterFile.[telescope label]:       [filename]
   //
   // which is assumed to contain identification grids.
 
   TString filename = gDataSet->GetDataSetEnv(Form("IdentificationParameterList.%s", GetLabel()));
 
   if (filename != "") {
      ReadIdentificationParameterFiles(filename.Data(), multidet);
   }
   else {
      filename = gDataSet->GetDataSetEnv(Form("IdentificationParameterFile.%s", GetLabel()));
 
      if (filename == "") {
         Warning("SetIdentificationParameters",
                 "No filename defined. Should be given by %s.IdentificationParameterFile.%s or %s.IdentificationParameterFile.%s",
                 gDataSet->GetName(), GetLabel(), gDataSet->GetName(), GetLabel());
         return kFALSE;
      }
 
      LoadIdentificationParameters(filename, multidet);
   }
   return kTRUE;
}
 
 
 
 
 
void KVIDTelescope::ReadIdentificationParameterFiles(const Char_t* filename, const KVMultiDetArray* multidet)
{
   // In the case where the identification grids are stored in several files, this method parse
   // the file found with the following environment variable:
   //
   // [dataset name].IdentificationParameterList.[telescope label]:       [filename]
   //
   // which contains the list of files containing the identification grids.
 
   KVFileReader fr;
   fr.OpenFileToRead(gDataSet->GetFullPathToDataSetFile(filename));
 
   while (fr.IsOK()) {
      fr.ReadLine(0);
 
      if (fr.GetCurrentLine() != "") LoadIdentificationParameters(fr.GetCurrentLine().Data(), multidet);
   }
 
   fr.CloseFile();
}
 
 
 
 
 
void KVIDTelescope::LoadIdentificationParameters(const Char_t* filename, const KVMultiDetArray* multidet)
{
   // This method add to the gIDGridManager list the identification grids.
 
   TString path;
 
   if ((path = gDataSet->GetFullPathToDataSetFile(filename)) == "") {
      Error("LoadIdentificationParameters",
            "File %s not found. Should be in %s",
            filename, gDataSet->GetDataSetDir());
      return;
   }
   //
   //Read grids from file
   Info("LoadIdentificationParameters", "Using file %s", path.Data());
   multidet->ReadGridsFromAsciiFile(path);
}
 
 
 
 
 
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(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;
   x = (Eres < 0 ? GetIDGridXCoord(grid) : 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(KVReconstructedNucleus* n)
{
   // For filtering simulations
   // Set the n->IsZMeasured() and n->IsAMeasured() status of the particle
   // In principle this depends on whether this telescope provides mass
   // identification or not, but this may depend on the particle's energy.
   // If A was not measured, it will be replaced with a value calculated
   // from whatever mass formula is used for the particle.
   //
   // 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
 
   n->SetZMeasured();
   if (!HasMassID()) {
      n->SetAMeasured(kFALSE);
      // beware - changing particle's mass changes its KE (momentum is conserved)
      double e = n->GetE();
      n->SetZ(n->GetZ());// use mass formula for A
      n->SetE(e);
   }
   else {
      if (fMassIDValidity) n->SetAMeasured(fMassIDValidity->Test(n)); // test expression for mass ID validity
      else n->SetAMeasured();   // no expression set; all nuclei are identified in mass
      if (!n->IsAMeasured()) {
         // beware - changing particle's mass changes its KE (momentum is conserved)
         double e = n->GetE();
         n->SetZ(n->GetZ()); // use mass formula for A
         n->SetE(e);
      }
   }
}
 
 
 
 
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);
}