KVDetector.cpp

#include "KVDetector.h"
#include "TROOT.h"
#include "KVGroup.h"
#include "KVCalibrator.h"
#include "TPluginManager.h"
#include "TClass.h"
/*** geometry ***/
#include "TGeoVolume.h"
#include "TGeoManager.h"
#include "TGeoMatrix.h"
#include "TGeoBBox.h"
#include "TGeoArb8.h"
#include "TGeoTube.h"
#include "KVCalibratedSignal.h"
#include "KVDetectorSignalTrace.h"
#include "KVDetectorSignalExpression.h"
#include "KVZDependentCalibratedSignal.h"
#include "KVMaterialStack.h"
#include <TGeoPhysicalNode.h>
#include <TGraph.h>
 
using namespace std;
 
ClassImp(KVDetector)
 
Int_t KVDetector::fDetCounter = 0;
 
 
 
void KVDetector::init()
{
   //default initialisations
   fCalibrators = nullptr;
   fParticles.SetCleanup();
   fActiveLayer = nullptr;
   fIdentP = fUnidentP = 0;
   fELossF = fEResF = fRangeF = nullptr;
   fEResforEinc = -1.;
   fSimMode = kFALSE;
   fPresent = kTRUE;
   fDetecting = kTRUE;
   fParentStrucList.SetCleanup();
   fSingleLayer = kTRUE;
   fNode.SetDetector(this);
   // detector owns any signals which are defined for it
   fDetSignals.SetOwner();
   // adding a new signal with the same name as an existing one
   // will delete the existing signal and replace it
   fDetSignals.ReplaceObjects();
}
 
 
 
 
KVDetector::KVDetector()
{
//default ctor
   init();
   fDetCounter++;
   SetName(Form("Det_%d", fDetCounter));
}
 
 
 
 
KVDetector::KVDetector(const Char_t* type,
                       const Float_t thick): KVMaterial()
{
   // Create a new detector of a given material and thickness in centimetres (default value = 0.0)
 
   init();
   SetType("DET");
   fDetCounter++;
   SetName(Form("Det_%d", fDetCounter));
   AddAbsorber(new KVMaterial(type, thick));
}
 
 
 
 
KVDetector::KVDetector(const Char_t* gas, const Double_t thick, const Double_t pressure, const Double_t temperature)
{
   // Create gaseous dteector with given type, linear thickness in cm, pressure in Torr, and temperature in degrees C (default value 19°C).
   //
   // \note this just defines some gas, the 'detector' has no windows!
 
   init();
   SetType("DET");
   fDetCounter++;
   SetName(Form("Det_%d", fDetCounter));
   AddAbsorber(new KVMaterial(gas, thick, pressure, temperature));
}
 
 
 
 
KVDetector::KVDetector(const KVDetector& obj) : KVMaterial(), KVPosition()
{
//copy ctor
   init();
   fDetCounter++;
   SetName(Form("Det_%d", fDetCounter));
#if ROOT_VERSION_CODE >= ROOT_VERSION(3,4,0)
   obj.Copy(*this);
#else
   ((KVDetector&) obj).Copy(*this);
#endif
}
 
 
#if ROOT_VERSION_CODE >= ROOT_VERSION(3,4,0)
 
 
void KVDetector::Copy(TObject& obj) const
#else
void KVDetector::Copy(TObject& obj)
#endif
{
   //copy 'this' to 'obj'
   //The structure of the detector is copied, with new cloned objects for
   //each absorber layer. The active layer is set in the new detector.
 
   TIter next(&fAbsorbers);
   KVMaterial* mat;
   while ((mat = (KVMaterial*) next())) {
      ((KVDetector&) obj).AddAbsorber((KVMaterial*) mat->Clone());
   }
   //set active layer
   Int_t in_actif = fAbsorbers.IndexOf(fActiveLayer);
   ((KVDetector&) obj).SetActiveLayer(((KVDetector&)obj).GetAbsorber(in_actif));
}
 
 
 
 
 
KVDetector::~KVDetector()
{
   SafeDelete(fCalibrators);
   SafeDelete(fELossF);
   SafeDelete(fEResF);
   SafeDelete(fRangeF);
}
 
 
 
 
void KVDetector::SetMaterial(const Char_t* type)
{
   // Set material of active layer.
   // If no absorbers have been added to the detector, create and add
   // one (active layer by default)
 
   if (!GetActiveLayer())
      AddAbsorber(new KVMaterial(type));
   else
      GetActiveLayer()->SetMaterial(type);
}
 
 
 
 
void KVDetector::DetectParticle(KVNucleus* kvp, TVector3* norm)
{
   //Calculate the energy loss of a charged particle traversing the detector,
   //the particle is slowed down, it is added to the list of all particles hitting the
   //detector.
   //
   //The apparent energy loss of the particle in the active layer of the
   //detector is set.
   //
   //If the optional argument 'norm' is given, it is supposed to be a vector
   //normal to the detector, oriented from the origin towards the detector.
   //
   //In this case the effective thicknesses of the detector's absorbers 'seen' by the particle
   //depending on its direction of motion is used for the calculation.
   //
   //Note that if the detector is part of an array geometry, GetEntranceWindow().GetNormalIntoShape()
   // returns a TVector3 unit vector which is the normal to the detector.
 
   if (kvp->GetKE() <= 0.)
      return;
 
   AddHit(kvp);                 //add nucleus to list of particles hitting detector in the event
   //set flag to say that particle has been slowed down
   kvp->SetIsDetected();
   //If this is the first absorber that the particle crosses, we set a "reminder" of its
   //initial energy
   if (!kvp->GetPInitial())
      kvp->SetE0();
 
   Double_t eloss(0), dE(0);
   if (norm)
   {
       KVMaterialStack stack(this);
       stack.SetIncidenceAngle(TMath::ACos(norm->Unit().Dot(kvp->GetMomentum().Unit())*TMath::RadToDeg()));
       eloss = stack.GetTotalDeltaE(kvp->GetZ(), kvp->GetA(), kvp->GetEnergy());
       dE = stack.GetDeltaE(kvp->GetZ(), kvp->GetA(), kvp->GetEnergy());
   }
   else
   {
       eloss = GetTotalDeltaE(kvp->GetZ(), kvp->GetA(), kvp->GetEnergy());
       dE = GetDeltaE(kvp->GetZ(), kvp->GetA(), kvp->GetEnergy());
   }
   auto epart = kvp->GetEnergy() - eloss;
   if (epart < 1e-3) {
      //printf("%s, pb d arrondi on met l energie de la particule a 0\n",GetName());
      epart = 0.0;
   }
   kvp->SetEnergy(epart);
   auto eloss_old = GetEnergyLoss();
   SetEnergyLoss(eloss_old + dE);
}
 
 
 
 
 
Double_t KVDetector::GetELostByParticle(KVNucleus* kvp, TVector3* norm)
{
   //Calculate the total energy loss of a charged particle traversing the detector.
   //
   //This does not affect the "stored" energy loss value of the detector nor the energy of the particle.
   //
   //If the optional argument 'norm' is given, it is supposed to be a vector
   //normal to the detector, oriented from the origin towards the detector.
   //
   //In this case the effective thicknesses of the detector's absorbers 'seen' by the particle
   //depending on its direction of motion is used for the calculation.
 
   if (norm)
   {
       KVMaterialStack stack(this);
       stack.SetIncidenceAngle(TMath::ACos(norm->Unit().Dot(kvp->GetMomentum().Unit())*TMath::RadToDeg()));
       return stack.GetTotalDeltaE(kvp->GetZ(), kvp->GetA(), kvp->GetEnergy());
   }
   return GetTotalDeltaE(kvp->GetZ(), kvp->GetA(), kvp->GetEnergy());
}
 
 
 
 
 
Double_t KVDetector::GetParticleEIncFromERes(KVNucleus* kvp, TVector3* norm)
{
   //Calculate the energy of particle 'kvn' before its passage through the detector,
   //based on the current kinetic energy, Z & A of nucleus 'kvn', supposed to be
   //after passing through the detector.
   //
   //If the optional argument 'norm' is given, it is supposed to be a vector
   //normal to the detector, oriented from the origin towards the detector.
   //
   //In this case the effective thicknesses of the detector's absorbers 'seen' by the particle
   //depending on its direction of motion is used for the calculation.
 
   Double_t Einc = 0.0;
   //make 'clone' of particle
   KVNucleus clone_part(kvp->GetZ(), kvp->GetA());
   clone_part.SetMomentum(kvp->GetMomentum());
   //composite detector - calculate losses in all layers
   KVMaterial* abs;
   TIter next(&fAbsorbers, kIterBackward); //work through list backwards
   while ((abs = (KVMaterial*) next())) {
 
      //calculate energy of particle before current absorber
      Einc = abs->GetParticleEIncFromERes(&clone_part, norm);
 
      //set new energy of particle
      clone_part.SetKE(Einc);
   }
   return Einc;
}
 
 
 
 
 
void KVDetector::Print(Option_t* opt) const
{
   //Print info on this detector
   //if option="data" the energy loss and raw data are displayed
 
   if (!strcmp(opt, "data")) {
      cout << Form("%10s -- ", GetName());
      // print raw data signals
      KVDetectorSignal* sig;
      TIter it(&GetListOfDetectorSignals());
      while ((sig = (KVDetectorSignal*)it())) {
         if (sig->IsRaw()) {
                std::cout << sig->GetName() << "=" << sig->GetValue() << " | ";
         }
      }
      // print calibrated data
      it.Reset();
      while ((sig = (KVDetectorSignal*)it())) {
         if (!sig->IsRaw() && !sig->GetValueNeedsExtraParameters()) {
            std::cout << sig->GetName() << "=" << sig->GetValue() << " | ";
         }
      }
      cout << " ";
      if (BelongsToUnidentifiedParticle())
         cout << "(Belongs to an unidentified particle)";
      cout << endl;
   }
   else if (!strcmp(opt, "all")) {
      //Give full details of detector
      //
      TString option("    ");
      cout << option << ClassName() << " : " << ((KVDetector*) this)->
           GetName() << endl;
      //composite detector - print all layers
      KVMaterial* abs;
      TIter next(&fAbsorbers);
      while ((abs = (KVMaterial*) next())) {
         if (GetActiveLayer() == abs)
            cout << " ### ACTIVE LAYER ### " << endl;
         cout << option << option;
         abs->Print();
         if (GetActiveLayer() == abs)
            cout << " #################### " << endl;
      }
      if (fParticles.GetEntries()) {
         cout << option << " --- Detected particles:" << endl;
         fParticles.Print();
      }
   }
   else {
      //just print name
      cout << ClassName() << ": " << ((KVDetector*) this)->
           GetName() << endl;
   }
}
 
 
 
 
Bool_t KVDetector::AddCalibrator(KVCalibrator* cal, const KVNameValueList& opts)
{
   // Associate a calibration with this detector.
   //
   // This will add a new signal to the list of the detector's signals.
   //
   // Also sets calibrator's name to `[detname]_[caltype]` where `caltype` is the type of the KVCalibration object.
   //
   // \param[in] cal pointer to KVCalibrator object (must be on heap, i.e. created with new: detector handles deletion)
   // \param[in] opts
   // \parblock
   // can be used to pass any extra parameters/options needed by the calibrator.
   // For example, if it contains a parameter `ZRange`:
   //
   //~~~~~~~~~~~~~~~
   // ZRange=1-10
   //~~~~~~~~~~~~~~~
   //
   // then the calibrator will be handled by a KVZDependentCalibratedSignal (handles several
   // calibrators which provide the same output signal, each one is used for a specific range
   // of atomic numbers)
   // \endparblock
   //
   // \returns kFALSE in case of problems (non-existent input signal for calibrator, output signal not defined for calibrator),
   // otherwise kTRUE
 
   if (!cal) return kFALSE;
   // look for input signal
   KVDetectorSignal* in  = GetDetectorSignal(cal->GetInputSignalType());
   // input signal must exist
   if (!in) {
      KVError::Error(this, "AddCalibrator", "Detector:%s has no signal %s, required as input by calibrator %s to provide output %s",
            GetName(), cal->GetInputSignalType().Data(), cal->GetType(), cal->GetOutputSignalType().Data());
      return kFALSE;
   }
   if (!fCalibrators)
      fCalibrators = new KVList();
 
   cal->SetDetector(this);
   fCalibrators->Add(cal);
   cal->SetName(Form("%s_%s", GetName(), cal->GetType()));
 
   if (cal->GetOutputSignalType() == "") {
      KVError::Warning(this, "AddCalibrator", "%s : output signal not defined for calibrator %s. No output signal created.",
              GetName(), cal->GetType());
   }
   else {
      KVDetectorSignal* new_cal_sig(nullptr);
      if (opts.HasParameter("ZRange")) {
         // If 'ZRange' parameter is given we need to find the KVZDependentCalibratedSignal
         // and add a new signal to it.
         KVDetectorSignal* sig = GetDetectorSignal(cal->GetOutputSignalType());
         if (!sig) {
            new_cal_sig = sig = new KVZDependentCalibratedSignal(in, cal->GetOutputSignalType());
         }
         dynamic_cast<KVZDependentCalibratedSignal*>(sig)->AddSignal(new KVCalibratedSignal(in, cal), opts.GetStringValue("ZRange"));
      }
      else {
         new_cal_sig = new KVCalibratedSignal(in, cal);
      }
      if (new_cal_sig) fDetSignals.Add(new_cal_sig);
   }
   return kTRUE;
}
 
 
 
 
Bool_t KVDetector::ReplaceCalibrator(const Char_t* type, KVCalibrator* cal, const KVNameValueList& opts)
{
   // Replace calibrator of given type with the given calibrator object
   // The calibrator object should not be shared with any other detectors: it now belongs
   // to this detector, which will delete it when necessary.
   // If an exising calibrator with the same type is already defined, it will be
   // deleted and removed from the detector's calibrator list
   //
   // Returns kFALSE in case of problems.
   //
   // The (optional) KVNameValueList argument can be used to pass any extra parameters/options.
 
   KVCalibrator* old_cal = GetCalibrator(type);
   if (old_cal) {
      fCalibrators->Remove(old_cal);
      remove_signal_for_calibrator(old_cal);
      delete old_cal;
   }
   return AddCalibrator(cal, opts);
}
 
 
 
 
Bool_t KVDetector::IsCalibrated(const KVNameValueList& params) const
{
   // A detector is considered to be calibrated if it has
   // a signal "Energy" available and if depending on the supplied parameters
   // this signal can be calculated
 
   if(IsSimMode())
      return IsDetecting();
 
   KVCalibratedSignal* e_sig = dynamic_cast<KVCalibratedSignal*>(GetDetectorSignal("Energy"));
   return (e_sig && e_sig->IsAvailableFor(params) && IsOK());
}
 
 
 
 
 
void KVDetector::Clear(Option_t* opt)
{
   //Set energy loss(es) etc. to zero
   //
   //If opt="N":
   //   + we do not reset acquisition parameters/raw detector signals
   //   + in SimMode we do not reset energy losses (this is called before reconstruction happens)
 
   bool option_N = (strncmp(opt, "N", 1) == 0);
   bool sim_mode_option_N = option_N && IsSimMode();
 
   if (!sim_mode_option_N) KVMaterial::Clear(opt);
   SetAnalysed(kFALSE);
   fIdentP = fUnidentP = 0;
   ResetBit(kIdentifiedParticle);
   ResetBit(kUnidentifiedParticle);
   if (!option_N) {
      TIter it(&GetListOfDetectorSignals());
      KVDetectorSignal* ds;
      while ((ds = (KVDetectorSignal*)it())) {
         if (ds->IsRaw() && !ds->IsExpression()) ds->Reset();
      }
   }
   ClearHits();
   if (!sim_mode_option_N) {
      //reset all layers in detector
      KVMaterial* mat;
      TIter next(&fAbsorbers);
      while ((mat = (KVMaterial*) next())) {
         mat->Clear();
      }
   }
   fEResforEinc = -1.;
}
 
 
 
 
void KVDetector::AddAbsorber(KVMaterial* mat)
{
   // Add a layer of absorber material to the detector
   // By default, the first layer added is set as the "Active" layer.
   // Call SetActiveLayer to change this.
   fAbsorbers.Add(mat);
   if (!TestBit(kActiveSet))
      SetActiveLayer(mat);
   if (fAbsorbers.GetSize() > 1) fSingleLayer = kFALSE;
}
 
 
 
 
KVMaterial* KVDetector::GetAbsorber(Int_t i) const
{
   //Returns pointer to the i-th absorber in the detector (i=0 first absorber, i=1 second, etc.)
 
   if (!fAbsorbers.GetEntries()) {
      KVError::Error(this, "GetAbsorber", "No absorbers defined for detector");
      return nullptr;
   }
   if (i < 0 || i >= fAbsorbers.GetEntries()) {
      KVError::Error(this, "GetAbsorber", "i=%d but valid values are 0-%d [number of absorbers in list]", i,
            fAbsorbers.GetEntries());
      return nullptr;
   }
   return (KVMaterial*) fAbsorbers.At(i);
}
 
 
 
 
void KVDetector::RemoveAllAbsorbers()
{
   // Completely reset the KVDetector as if it had just been created by a call to the default constructor:
   //   + removes all absorber layers
   //   + removes all calibrators/detector signals
 
   Clear();
   RemoveCalibrators();
   SafeDelete(fCalibrators);
   fAbsorbers.Clear();
   SafeDelete(fELossF);
   SafeDelete(fEResF);
   SafeDelete(fRangeF);
   fCalibrators = nullptr;
   fActiveLayer = nullptr;
   ResetBit(kActiveSet);
   fIdentP = fUnidentP = 0;
   fELossF = fEResF = fRangeF = nullptr;
   fEResforEinc = -1.;
   fSimMode = kFALSE;
   fPresent = kTRUE;
   fDetecting = kTRUE;
   fSingleLayer = kTRUE;
}
 
 
 
 
void KVDetector::remove_signal_for_calibrator(KVCalibrator* K)
{
   // Used when a calibrator object is removed or replaced
   // We remove and delete the corresponding output signal from the list of detector signals
 
   if (K->GetOutputSignalType() != "") {
      KVDetectorSignal* ds = GetDetectorSignal(K->GetOutputSignalType());
      if (ds) {
         fDetSignals.Remove(ds);
         delete ds;
      }
   }
}
 
 
 
 
void KVDetector::RemoveCalibrators()
{
   // Removes all calibrations associated to this detector: in other words, we delete all
   // the KVCalibrator objects in list fCalibrators.
   //
   // We also destroy all signals provided by these calibrators
 
   if (fCalibrators) {
      KVCalibrator* K;
      TIter it(fCalibrators);
      while ((K = (KVCalibrator*)it())) {
         remove_signal_for_calibrator(K);
      }
      fCalibrators->Delete();
   }
}
 
 
 
 
 
Double_t KVDetector::GetCorrectedEnergy(KVNucleus* nuc, Double_t e, Bool_t transmission)
{
   // \returns the calculated *total* energy loss in the detector for a given nucleus including inactive absorber layers.
   // 
   // \param[in] e [MeV] energy loss of nucleus in active layer - if not given, we use current value
   // \param[in] transmission
   //               + if kTRUE (default): the particle is assumed to emerge with a non-zero residual energy \f$E_{res}\f$ after the detector
   //               + if kFALSE: the particle is assumed to stop in the detector.
   //
   // \warning if transmission=kTRUE, and if the residual energy after the detector
   //   is known (i.e. measured in a detector placed after this one), you should
   //   first call method SetEResAfterDetector() with the value of the residual energy
   //   before calling this method. Otherwise, especially for heavy ions, the
   //   correction may be false for particles which are just above the punch-through energy.
   //
   // #### Angle of incidence and effective thickness
   // It is possible that the measured energy loss in the active layer is greater than the maximum
   //   theoretical value calculated from range tables. Apart from the uncertainty in the actual measurement
   //   (plus the uncertainty of the detector calibration), for any detector covering a finite solid angle
   //   the angles of incidence of particles hitting the detector will have a distribution, with very
   //   few (even none) arriving in a direction perpendicular to the detector's surface, i.e. in the direction of
   // the axis along which the detector's nominal thickness is defined.
   //
   //  In this case, the effective thickness of the detector is increased by a factor \f$1/\cos\psi\f$, where \f$\psi\f$
   //   is the angle between the particle's trajectory and the normal to the detector. We take this into account when
   //   trying to decide if the measured energy loss is consistent with the nucleus in question, by using a
   //   KVMaterialStack.
   //
   //   KVMaterialStack can calculate the minimum value of \f$\psi_{\Delta E}\f$ that would be compatible
   //   with a given energy loss (i.e. that for which the maximum possible energy loss is that which was measured): see KVMaterialStack::GetMinimumIncidentAngleForDEMax().
   //   We can then check if such an angle is compatible with the geometry of the detector - GetMaximumMisalignmentAngle()
   //   gives the largest possible value, \f$\psi_{max}\f$, for particles impinging on the detector.
   //
   //   If a minimum angle of incidence \f$\psi_{\Delta E}\f$ which is consistent with the geometry can be found,
   //   it is stored in the parameters of the reconstructed nucleus with name `GetCorrectedEnergy.IncidenceAngle`
   //   i.e. for any reconstructed nucleus you can do
   //
   //~~~{.cpp}
   // KVReconstructedNucleus& rnuc; // reference to a reconstructed nucleus
   // if( rnuc.GetParameters()->HasParameter("GetCorrectedEnergy.IncidenceAngle") )
   // {
   //    auto psi = rnuc.GetParameters()->GetDoubleValue("GetCorrectedEnergy.IncidenceAngle");
   // }
   //~~~
   //
   //  If \f$\psi_{\Delta E}>\psi_{max}\f$, we use the (apparently impossible) incidence angle in order to calculate the
   // corrected energy, and we add the following warning parameters to the particle, each of which begins with the
   // name of the detector (unlike the previous one):
   //
   //~~~
   // [det].GetCorrectedEnergy.Warning = 1
   // [det].GetCorrectedEnergy.MeasuredDE = [value]
   // [det].GetCorrectedEnergy.MaxDE = [value]
   // [det].GetCorrectedEnergy.Transmission = 0 or 1
   // [det].GetCorrectedEnergy.ERES = [value]
   // [det].GetCorrectedEnergy.RequiredIncidenceAngle = [value]
   // [det].GetCorrectedEnergy.MaximumIncidenceAngle = [value]
   //~~~
 
   Int_t z = nuc->GetZ();
   Int_t a = nuc->GetA();
 
   if (e < 0.) e = GetEnergy();
   if (e <= 0) {
      SetEResAfterDetector(-1.);
      return 0;
   }
 
   auto get_stack_correct_energy = [&](KVMaterialStack& stack, double maxDE, bool save_angle_if_good = false)
   {
       auto inc_angle = stack.GetMinimumIncidentAngleForDEMax(z, a, e);
       stack.SetIncidenceAngle(inc_angle);
       // check incidence angle is consistent with geometry
       if(inc_angle > GetEntranceWindow().GetMaximumMisalignmentAngle())
       {
           nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.Warning", GetName()), 1);
           nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.MeasuredDE", GetName()), e);
           nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.MaxDE", GetName()), maxDE);
           nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.Transmission", GetName()), (Int_t)transmission);
           nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.ERES", GetName()), GetEResAfterDetector());
           nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.RequiredIncidenceAngle", GetName()), inc_angle);
           nuc->GetParameters()->SetValue(Form("%s.GetCorrectedEnergy.MaximumIncidenceAngle", GetName()), GetEntranceWindow().GetMaximumMisalignmentAngle());
       }
       else
       {
           if(save_angle_if_good)
           {
               // if the incidence angle is consistent with geometry, we keep it for use with other detectors on the particle's
               // trajectory - in principle it should be the same for all
               nuc->GetParameters()->SetValue("GetCorrectedEnergy.IncidenceAngle", inc_angle);
           }
       }
   };
   // check if calling this function for a previous detector in the particle's trajectory
   // led to using an effective incident angle for the particle: if so, use it here
   if (nuc->GetParameters()->HasDoubleParameter("GetCorrectedEnergy.IncidenceAngle")) {
       KVMaterialStack stack(this);
       stack.SetIncidenceAngle(nuc->GetParameters()->GetDoubleValue("GetCorrectedEnergy.IncidenceAngle"));
       // check that apparent energy loss in detector is compatible with a & z
       Double_t maxDE = stack.GetMaxDeltaE(z, a);
       if (e > maxDE)
           get_stack_correct_energy(stack,maxDE);
       return get_corrected_energy(&stack, nuc, e, transmission);
   }
   // check that apparent energy loss in detector is compatible with a & z
   Double_t maxDE = GetMaxDeltaE(z, a);
   if (e > maxDE) {
       KVMaterialStack stack(this);
       get_stack_correct_energy(stack,maxDE,true);
       return get_corrected_energy(&stack, nuc, e, transmission);
   }
   return get_corrected_energy(this, nuc, e, transmission);
}
 
 
 
 
 
Int_t KVDetector::FindZmin(Double_t ELOSS, Char_t mass_formula)
{
   //For particles which stop in the first stage of an identification telescope,
   //we can at least estimate a minimum Z value based on the energy lost in this
   //detector.
   //
   //This is based on the KVMaterial::GetMaxDeltaE method, giving the maximum
   //energy loss in the active layer of the detector for a given nucleus (A,Z).
   //
   //The "Zmin" is the Z of the particle which gives a maximum energy loss just greater
   //than that measured in the detector. Particles with Z<Zmin could not lose as much
   //energy and so are excluded.
   //
   //If ELOSS is not given, we use the current value of GetEnergy()
   //Use 'mass_formula' to change the formula used to calculate the A of the nucleus
   //from its Z. Default is valley of stability value. (see KVNucleus::GetAFromZ).
   //
   //If the value of ELOSS or GetEnergy() is <=0 we return Zmin=0
 
   ELOSS = (ELOSS < 0. ? GetEnergy() : ELOSS);
   if (ELOSS <= 0) return 0;
 
   UInt_t z = 40;
   UInt_t zmin, zmax;
   zmin = 1;
   zmax = 92;
   KVNucleus particle(zmin);
 
   if (mass_formula > -1)
      particle.SetMassFormula((UChar_t)mass_formula);
 
   if (GetMaxDeltaE(particle.GetZ(), particle.GetA()) > ELOSS)
      return zmin; // Z=1 solution
 
   do {
 
      particle.SetZ(z);
 
      auto difference = GetMaxDeltaE(z, particle.GetA()) - ELOSS;
      //if difference < 0 the z is too small
      if (difference < 0.0) {
 
         zmin = z;
         z += (UInt_t)((zmax - z) / 2 + 0.5);
 
      }
      else {
 
         zmax = z;
         z -= (UInt_t)((z - zmin) / 2 + 0.5);
 
      }
   }
   while (zmax > (zmin + 1));
 
   return zmax;
}
 
 
 
 
 
Double_t KVDetector::ELossActive(Double_t* x, Double_t* par)
{
   // Calculates energy loss (in MeV) in active layer of detector, taking into account preceding layers
   //
   // Arguments are:
   //    x[0] is incident energy in MeV
   // Parameters are:
   //   par[0]   Z of ion
   //   par[1]   A of ion
 
   Double_t e = x[0];
   TIter next(&fAbsorbers);
   KVMaterial* mat;
   mat = (KVMaterial*)next();
   while (fActiveLayer != mat) {
      // calculate energy losses in absorbers before active layer
      e = mat->GetERes(par[0], par[1], e);     //residual energy after layer
      if (e <= 0.)
         return 0.;          // return 0 if particle stops in layers before active layer
      mat = (KVMaterial*)next();
   }
   //calculate energy loss in active layer
   return mat->GetDeltaE(par[0], par[1], e);
}
 
 
 
 
 
Double_t KVDetector::RangeDet(Double_t* x, Double_t* par)
{
   // Calculates range (in centimetres) of ions in detector as a function of incident energy (in MeV),
   // taking into account all layers of the detector.
   //
   // Arguments are:
   //   x[0]  =  incident energy in MeV
   // Parameters are:
   //   par[0] = Z of ion
   //   par[1] = A of ion
 
   Double_t Einc = x[0];
   Int_t Z = (Int_t)par[0];
   Int_t A = (Int_t)par[1];
   Double_t range = 0.;
   TIter next(&fAbsorbers);
   KVMaterial* mat = (KVMaterial*)next();
   if (!mat) return 0.;
   do {
      // range in this layer
      Double_t this_range = mat->GetLinearRange(Z, A, Einc);
      KVMaterial* next_mat = (KVMaterial*)next();
      if (this_range > mat->GetThickness()) {
         // particle traverses layer.
         if (next_mat)
            range += mat->GetThickness();
         else // if this is last layer, the range continues to increase beyond size of detector
            range += this_range;
         // calculate incident energy for next layer (i.e. residual energy after this layer)
         Einc = mat->GetERes(Z, A, Einc);
      }
      else {
         // particle stops in layer
         range += this_range;
         return range;
      }
      mat = next_mat;
   }
   while (mat);
   // particle traverses detector
   return range;
}
 
 
 
 
 
Double_t KVDetector::EResDet(Double_t* x, Double_t* par)
{
   // Calculates residual energy (in MeV) of particle after traversing all layers of detector.
   // Returned value is -1000 if particle stops in one of the layers of the detector.
   //
   // Arguments are:
   //    x[0] is incident energy in MeV
   // Parameters are:
   //   par[0]   Z of ion
   //   par[1]   A of ion
 
   Double_t e = x[0];
   TIter next(&fAbsorbers);
   KVMaterial* mat;
   while ((mat = (KVMaterial*)next())) {
      Double_t eres = mat->GetERes(par[0], par[1], e);     //residual energy after layer
      if (eres <= 0.)
         return -1000.;          // return -1000 if particle stops in layers before active layer
      e = eres;
   }
   return e;
}
 
 
 
 
 
KVDetector* KVDetector::MakeDetector(const Char_t* name, Float_t thickness)
{
   //Static function which will create an instance of the KVDetector-derived class
   //corresponding to 'name'
   //These are defined as 'Plugin' objects in the file $KVROOT/KVFiles/.kvrootrc :
   //   [name_of_dataset].detector_type
   //   detector_type
   // To use the dataset-dependent plugin, call this method with
   //  name = "[name_of_dataset].detector_type"
   // If not, the default plugin will be used
   //first we check if there is a special plugin for the DataSet
   //if not we take the default one
   //
   //'thickness' is passed as argument to the constructor for the detector plugin
 
   //check and load plugin library
   TPluginHandler* ph = nullptr;
   KVString nom(name);
   if (!nom.Contains(".") && !(ph = LoadPlugin("KVDetector", name))) return nullptr;
   if (nom.Contains(".")) {
      // name format like [dataset].[det_type]
      // in case dataset name contains "." we parse string to find detector type: assumed after the last "."
      nom.RBegin(".");
      KVString det_type = nom.RNext();
      if (!(ph = LoadPlugin("KVDetector", name))) {
         if (!(ph = LoadPlugin("KVDetector", det_type))) {
            return nullptr;
         }
      }
   }
 
   //execute constructor/macro for detector
   return (KVDetector*) ph->ExecPlugin(1, thickness);
}
 
 
 
 
Double_t KVDetector::GetEntranceWindowSurfaceArea()
{
   // Return surface area of first layer of detector in cm2.
 
   if (!GetEntranceWindow().GetShape()->InheritsFrom("TGeoArb8")) {
      // simple shape area
      return GetEntranceWindow().GetShape()->GetFacetArea(1);
   }
   // Monte Carlo calculation for TGeoArb8 shapes
   return GetEntranceWindow().GetSurfaceArea();
}
 
 
 
 
TF1* KVDetector::GetEResFunction(Int_t Z, Int_t A)
{
   // Return pointer toTF1 giving residual energy after detector as function of incident energy,
   // for a given nucleus (Z,A).
   // The TF1::fNpx parameter is taken from environment variable KVDetector.ResidualEnergy.Npx
 
   if (!fEResF) {
      fEResF = new TF1(Form("KVDetector:%s:ERes", GetName()), this, &KVDetector::EResDet,
                       0., 1.e+04, 2, "KVDetector", "EResDet");
      fEResF->SetNpx(gEnv->GetValue("KVDetector.ResidualEnergy.Npx", 20));
   }
   fEResF->SetParameters((Double_t)Z, (Double_t)A);
   fEResF->SetRange(0., GetSmallestEmaxValid(Z, A));
   fEResF->SetTitle(Form("Residual energy [MeV] after detector %s for Z=%d A=%d", GetName(), Z, A));
 
   return fEResF;
}
 
 
 
 
TF1* KVDetector::GetRangeFunction(Int_t Z, Int_t A)
{
   // Return pointer toTF1 giving range (in centimetres) in detector as function of incident energy,
   // for a given nucleus (Z,A).
   // The TF1::fNpx parameter is taken from environment variable KVDetector.Range.Npx
 
   if (!fRangeF) {
      fRangeF = new TF1(Form("KVDetector:%s:Range", GetName()), this, &KVDetector::RangeDet,
                        0., 1.e+04, 2, "KVDetector", "RangeDet");
      fRangeF->SetNpx(gEnv->GetValue("KVDetector.Range.Npx", 20));
   }
   fRangeF->SetParameters((Double_t)Z, (Double_t)A);
   fRangeF->SetRange(0., GetSmallestEmaxValid(Z, A));
   fRangeF->SetTitle(Form("Range [cm] in detector %s for Z=%d A=%d", GetName(), Z, A));
 
   return fRangeF;
}
 
 
 
 
TF1* KVDetector::GetELossFunction(Int_t Z, Int_t A)
{
   // Return pointer to TF1 giving energy loss in active layer of detector as function of incident energy,
   // for a given nucleus (Z,A).
   // The TF1::fNpx parameter is taken from environment variable KVDetector.EnergyLoss.Npx
 
   if (!fELossF) {
      fELossF = new TF1(Form("KVDetector:%s:ELossActive", GetName()), this, &KVDetector::ELossActive,
                        0., 1.e+04, 2, "KVDetector", "ELossActive");
      fELossF->SetNpx(gEnv->GetValue("KVDetector.EnergyLoss.Npx", 20));
   }
   fELossF->SetParameters((Double_t)Z, (Double_t)A);
   fELossF->SetRange(0., GetSmallestEmaxValid(Z, A));
   fELossF->SetTitle(Form("Energy loss [MeV] in detector %s for Z=%d A=%d", GetName(), Z, A));
   return fELossF;
}
 
 
 
 
Double_t KVDetector::GetEIncOfMaxDeltaE(Int_t Z, Int_t A)
{
   // Returns incident energy corresponding to maximum energy loss in the
   // active layer of the detector, for a given nucleus.
 
   return GetELossFunction(Z, A)->GetMaximumX();
}
 
 
 
 
Double_t KVDetector::GetMaxDeltaE(Int_t Z, Int_t A)
{
   // Returns maximum energy loss in the
   // active layer of the detector, for a given nucleus.
 
   return GetELossFunction(Z, A)->GetMaximum();
}
 
 
 
 
Double_t KVDetector::GetDeltaE(Int_t Z, Int_t A, Double_t Einc, Double_t)
{
   // Returns energy loss of given nucleus in the active layer of the detector.
 
   // optimization for single-layer detectors
   if (fSingleLayer) {
      return fActiveLayer->GetDeltaE(Z, A, Einc);
   }
   return GetELossFunction(Z, A)->Eval(Einc);
}
 
 
 
 
Double_t KVDetector::GetTotalDeltaE(Int_t Z, Int_t A, Double_t Einc)
{
   // Returns calculated total energy loss of ion in ALL layers of the detector.
   // This is just (Einc - GetERes(Z,A,Einc))
 
   return Einc - GetERes(Z, A, Einc);
}
 
 
 
 
Double_t KVDetector::GetERes(Int_t Z, Int_t A, Double_t Einc, Double_t)
{
   // Returns residual energy of given nucleus after the detector.
   // Returns 0 if Einc<=0
 
   if (Einc <= 0.) return 0.;
   Double_t eres = GetEResFunction(Z, A)->Eval(Einc);
   // Eres function returns -1000 when particle stops in detector,
   // in order for function inversion (GetEIncFromEres) to work
   if (eres < 0.) eres = 0.;
   return eres;
}
 
 
 
 
Double_t KVDetector::GetIncidentEnergy(Int_t Z, Int_t A, Double_t delta_e, enum SolType type)
{
   // Returns incident energy corresponding to energy loss delta_e in active layer of detector for a given nucleus.
   // If delta_e is not given, the current energy loss in the active layer is used.
   //
   // By default the solution corresponding to the highest incident energy is returned
   // This is the solution found for Einc greater than the maximum of the dE(Einc) curve.
   // If you want the low energy solution set SolType = KVIonRangeTable::kEmin.
   //
   // \warning Calculating the incident energy of a particle using only the dE in a detector
   // is ambiguous, as in general (and especially for very heavy ions) the maximum of the dE
   // curve occurs for Einc greater than the punch-through energy, therefore it is not always
   // true to assume that if the particle does not stop in the detector the required solution
   // is that for type=KVIonRangeTable::kEmax. For a range of energies between punch-through
   // and dE_max, the required solution is still that for type=KVMaterial::kEmin.
   // If the residual energy of the particle is unknown, there is no way to know which is the
   // correct solution.
   //
   // \warning If the given energy loss in the active layer is greater than the maximum theoretical dE
   // for given Z & A, (dE > GetMaxDeltaE(Z,A)) then we return a NEGATIVE incident energy
   // corresponding to the maximum, GetEIncOfMaxDeltaE(Z,A)
 
   if (Z < 1) return 0.;
 
   Double_t DE = (delta_e > 0 ? delta_e : GetEnergyLoss());
 
   // If the given energy loss in the active layer is greater than the maximum theoretical dE
   // for given Z & A, (dE > GetMaxDeltaE(Z,A)) then we return a NEGATIVE incident energy
   // corresponding to the maximum, GetEIncOfMaxDeltaE(Z,A)
   if (DE > GetMaxDeltaE(Z, A)) return -GetEIncOfMaxDeltaE(Z, A);
 
   TF1* dE = GetELossFunction(Z, A);
   Double_t e1, e2;
   dE->GetRange(e1, e2);
   switch (type) {
      case kEmin:
         e2 = GetEIncOfMaxDeltaE(Z, A);
         break;
      case kEmax:
         e1 = GetEIncOfMaxDeltaE(Z, A);
         break;
   }
   KVBase::GetX_status status;
   Double_t EINC = ProtectedGetX(dE, DE, status, e1, e2);
   if (status == KVBase::GetX_status::above_maximum || status == KVBase::GetX_status::below_minimum) {
//      Warning("GetIncidentEnergy",
//              "In %s : Called for Z=%d A=%d with dE=%f sol_type %d",
//              GetName(), Z, A, DE, type);
//      Warning("GetIncidentEnergy",
//              "Max DeltaE in this case is %f MeV",
//              GetMaxDeltaE(Z, A));
//      Warning("GetIncidentEnergy",
//              "Between Einc limits [%f,%f] no solution found",
//              e1,e2);
//      Warning("GetIncidentEnergy",
//              "Returned value is -Einc for %s dE, %f",
//              (status>0 ? "maximum" : "minimum"), EINC);
      return -EINC;
   }
   return EINC;
}
 
 
 
 
Double_t KVDetector::GetDeltaEFromERes(Int_t Z, Int_t A, Double_t Eres)
{
   // Overrides KVMaterial::GetDeltaEFromERes
   //
   // Calculate energy loss in active layer of detGetAlignedDetector for nucleus (Z,A)
   // having a residual kinetic energy Eres (MeV)
 
   if (Z < 1 || Eres <= 0.) return 0.;
   Double_t Einc = GetIncidentEnergyFromERes(Z, A, Eres);
   if (Einc <= 0.) return 0.;
   return GetELossFunction(Z, A)->Eval(Einc);
}
 
 
 
 
Double_t KVDetector::GetIncidentEnergyFromERes(Int_t Z, Int_t A, Double_t Eres)
{
   // Overrides KVMaterial::GetIncidentEnergyFromERes
   //
   // Calculate incident energy of nucleus from residual energy.
   //
   // Returns -1 if Eres is out of defined range of values
 
   if (Z < 1 || Eres <= 0.) return 0.;
   //return GetEResFunction(Z, A)->GetX(Eres);
   KVBase::GetX_status status;
   Double_t einc = KVBase::ProtectedGetX(GetEResFunction(Z, A), Eres, status);
   if (status == KVBase::GetX_status::above_maximum || status == KVBase::GetX_status::below_minimum){// problem with inversion - value out of defined range of function
      return -1;
   }
   return einc;
}
 
 
 
 
Double_t KVDetector::GetSmallestEmaxValid(Int_t Z, Int_t A) const
{
   // Returns the smallest maximum energy for which range tables are valid
   // for all absorbers in the detector, and given ion (Z,A)
 
   Double_t maxmin = -1.;
   TIter next(&fAbsorbers);
   KVMaterial* abs;
   while ((abs = (KVMaterial*)next())) {
      if (maxmin < 0.) maxmin = abs->GetEmaxValid(Z, A);
      else {
         if (abs->GetEmaxValid(Z, A) < maxmin) maxmin = abs->GetEmaxValid(Z, A);
      }
   }
   return maxmin;
}
 
 
 
 
void KVDetector::ReadDefinitionFromFile(const Char_t* envrc)
{
   // Create detector from text file in 'TEnv' format.
   //
   // Example:
   // ========
   //
   // Layer:  Gold
   // Gold.Material:   Au
   // Gold.AreaDensity:   200.*KVUnits::ug
   // +Layer:  Gas1
   // Gas1.Material:   C3F8
   // Gas1.Thickness:   5.*KVUnits::cm
   // Gas1.Pressure:   50.*KVUnits::mbar
   // Gas1.Active:    yes
   // +Layer:  Si1
   // Si1.Material:   Si
   // Si1.Thickness:   300.*KVUnits::um
 
   TEnv fEnvFile(envrc);
 
   KVString layers(fEnvFile.GetValue("Layer", ""));
   layers.Begin(" ");
   while (!layers.End()) {
      KVString layer = layers.Next();
      KVString mat = fEnvFile.GetValue(Form("%s.Material", layer.Data()), "");
      KVString tS = fEnvFile.GetValue(Form("%s.Thickness", layer.Data()), "");
      KVString pS = fEnvFile.GetValue(Form("%s.Pressure", layer.Data()), "");
      KVString dS = fEnvFile.GetValue(Form("%s.AreaDensity", layer.Data()), "");
      Double_t thick, dens, press;
      thick = dens = press = 0.;
      KVMaterial* M = 0;
      if (pS != "" && tS != "") {
         press = (Double_t)gROOT->ProcessLineFast(Form("%s*1.e+12", pS.Data()));
         press /= 1.e+12;
         thick = (Double_t)gROOT->ProcessLineFast(Form("%s*1.e+12", tS.Data()));
         thick /= 1.e+12;
         M = new KVMaterial(mat.Data(), thick, press);
      }
      else if (tS != "") {
         thick = (Double_t)gROOT->ProcessLineFast(Form("%s*1.e+12", tS.Data()));
         thick /= 1.e+12;
         M = new KVMaterial(mat.Data(), thick);
      }
      else if (dS != "") {
         dens = (Double_t)gROOT->ProcessLineFast(Form("%s*1.e+12", dS.Data()));
         dens /= 1.e+12;
         M = new KVMaterial(dens, mat.Data());
      }
      if (M) {
         AddAbsorber(M);
         if (fEnvFile.GetValue(Form("%s.Active", layer.Data()), kFALSE)) SetActiveLayer(M);
      }
   }
}
 
 
 
 
Double_t KVDetector::GetRange(Int_t Z, Int_t A, Double_t Einc)
{
   // WARNING: SAME AS KVDetector::GetLinearRange
   // Only linear range in centimetres is calculated for detectors!
   return GetLinearRange(Z, A, Einc);
}
 
 
 
 
Double_t KVDetector::GetLinearRange(Int_t Z, Int_t A, Double_t Einc)
{
   // Returns range of ion in centimetres in this detector,
   // taking into account all layers.
   // Note that for Einc > punch through energy, this range is no longer correct
   // (but still > total thickness of detector).
   return GetRangeFunction(Z, A)->Eval(Einc);
}
 
 
 
 
Double_t KVDetector::GetPunchThroughEnergy(Int_t Z, Int_t A)
{
   // Returns energy (in MeV) necessary for ion (Z,A) to punch through all
   // layers of this detector
 
   if (fSingleLayer) {
      // Optimize calculation time for single-layer detector
      return fActiveLayer->GetPunchThroughEnergy(Z, A);
   }
   //return GetRangeFunction(Z, A)->GetX(GetTotalThicknessInCM());
   KVBase::GetX_status status;
   return ProtectedGetX(GetRangeFunction(Z, A),GetTotalThicknessInCM(),status);
}
 
 
 
 
 
KVDrawable<TGraph> KVDetector::DrawPunchThroughEnergyVsZ(Int_t massform)
{
   // Creates and fills a KVDrawable<TGraph> with the punch through energy in MeV vs. Z for the given detector,
   // for Z=1-92. The mass of each nucleus is calculated according to the given mass formula
   // (see KVNucleus).
 
   TGraph* punch = new TGraph(92);
   punch->SetName(Form("KVDetpunchthrough_%s_mass%d", GetName(), massform));
   punch->SetTitle(Form("Simple Punch-through %s (MeV) (mass formula %d)", GetName(), massform));
   KVNucleus nuc;
   nuc.SetMassFormula(massform);
   for (int Z = 1; Z <= 92; Z++) {
      nuc.SetZ(Z);
      punch->SetPoint(Z - 1, Z, GetPunchThroughEnergy(nuc.GetZ(), nuc.GetA()));
   }
   return punch;
}
 
 
 
 
KVDrawable<TGraph> KVDetector::DrawPunchThroughEsurAVsZ(Int_t massform)
{
   // Creates and fills a KVDrawable<TGraph> with the punch through energy in MeV/nucleon vs. Z for the given detector,
   // for Z=1-92. The mass of each nucleus is calculated according to the given mass formula
   // (see KVNucleus).
 
   TGraph* punch = new TGraph(92);
   punch->SetName(Form("KVDetpunchthroughEsurA_%s_mass%d", GetName(), massform));
   punch->SetTitle(Form("Simple Punch-through %s (AMeV) (mass formula %d)", GetName(), massform));
   KVNucleus nuc;
   nuc.SetMassFormula(massform);
   for (int Z = 1; Z <= 92; Z++) {
      nuc.SetZ(Z);
      punch->SetPoint(Z - 1, Z, GetPunchThroughEnergy(nuc.GetZ(), nuc.GetA()) / nuc.GetA());
   }
   return punch;
}
 
 
 
 
 
KVGroup* KVDetector::GetGroup() const
{
   return (KVGroup*)GetParentStructure("GROUP");
}
 
 
 
 
UInt_t KVDetector::GetGroupNumber()
{
   return (GetGroup() ? GetGroup()->GetNumber() : 0);
}
 
 
 
 
void KVDetector::AddParentStructure(KVGeoStrucElement* elem)
{
   fParentStrucList.Add(elem);
}
 
 
 
 
void KVDetector::RemoveParentStructure(KVGeoStrucElement* elem)
{
   fParentStrucList.Remove(elem);
}
 
 
 
 
KVGeoStrucElement* KVDetector::GetParentStructure(const Char_t* type, const Char_t* name) const
{
   // Get parent geometry structure element of given type.
   // Give unique name of structure if more than one element of same type is possible.
   KVGeoStrucElement* el = 0;
   if (strcmp(name, "")) {
      KVSeqCollection* strucs = fParentStrucList.GetSubListWithType(type);
      el = (KVGeoStrucElement*)strucs->FindObject(name);
      delete strucs;
   }
   else
      el = (KVGeoStrucElement*)fParentStrucList.FindObjectByType(type);
   return el;
}
 
 
 
 
void KVDetector::SetActiveLayerMatrix(const TGeoHMatrix* m)
{
   // Set ROOT geometry global matrix transformation to coordinate frame of active layer volume
   SetMatrix(m);
}
 
 
 
 
void KVDetector::SetActiveLayerShape(TGeoBBox* s)
{
   // Set ROOT geometry shape of active layer volume
   SetShape(s);
}
 
 
 
 
void KVDetector::SetEntranceWindowMatrix(const TGeoHMatrix* m)
{
   // Set ROOT geometry global matrix transformation to coordinate frame of entrance window
   fEWPosition.SetMatrix(m);
}
 
 
 
 
void KVDetector::SetEntranceWindowShape(TGeoBBox* s)
{
   // Set ROOT geometry shape of entrance window
   fEWPosition.SetShape(s);
}
 
 
 
 
void KVDetector::SetThickness(Double_t thick)
{
   // Overrides KVMaterial::SetThickness
   //
   // If ROOT geometry is defined, we modify the DZ thickness of the volume representing
   // this detector in accordance with the new thickness.
   //
   // This is only implemented for single-layer detectors with the following shapes:
   //  - TGeoBBox (rectangular box)
   //  - TGeoTube (tube)
   //  - TGeoArb8 (arbitrary trapezoid with less than 8 vertices standing on two parallel planes perpendicular to Z axis)
 
   if (ROOTGeo() && fSingleLayer) {
      TGeoPhysicalNode* pn = (TGeoPhysicalNode*)gGeoManager->GetListOfPhysicalNodes()->FindObject(GetNode()->GetFullPathToNode());
      if (!pn) pn = gGeoManager->MakePhysicalNode(GetNode()->GetFullPathToNode());
      TGeoBBox* shape = (TGeoBBox*)pn->GetShape();
      TGeoShape* newshape = nullptr;
      // bad kludge - is there no better way to clone a shape and change its dZ?
      if (shape->IsA() == TGeoBBox::Class()) {
         newshape = new TGeoBBox(shape->GetDX(), shape->GetDY(), 0.5 * thick);
      }
      else if (shape->IsA() == TGeoTube::Class()) {
         TGeoTube* oldtube = static_cast<TGeoTube*>(shape);
         newshape = new TGeoTube(oldtube->GetRmin(), oldtube->GetRmax(), 0.5 * thick);
      }
      else if (shape->IsA() == TGeoArb8::Class()) {
         TGeoArb8* oldtube = static_cast<TGeoArb8*>(shape);
         auto vert = oldtube->GetVertices();
         newshape = new TGeoArb8(0.5 * thick, vert);
      }
      else {
         KVError::Error(this, "SetThickness", "No implementation for %s (%s)", shape->IsA()->GetName(), GetName());
      }
      if (newshape) {
         pn->Align(nullptr, newshape);
         gGeoManager->RefreshPhysicalNodes(kFALSE);
      }
   }
   KVMaterial::SetThickness(thick);
}
 
 
 
 
Bool_t KVDetector::HasSameStructureAs(const KVDetector* other) const
{
   // Return kTRUE if the two detectors have the same internal structure, i.e.
   //  - the same number of absorber layers
   //  - in the same order
   //  - with the same material & thickness
 
   int nabs = GetNumberOfAbsorberLayers();
   if (other->GetNumberOfAbsorberLayers() != nabs) return kFALSE;
   bool same = true;
   for (int iabs = 0; iabs < nabs; ++iabs) {
      KVMaterial* this_abs = GetAbsorber(iabs);
      KVMaterial* other_abs = other->GetAbsorber(iabs);
      if (!this_abs->IsType(other_abs->GetType())
            || this_abs->GetMass() != other_abs->GetMass()
            || this_abs->GetThickness() != other_abs->GetThickness())
         same = false;
   }
   return same;
}
 
 
 
 
KVDetectorSignal* KVDetector::AddDetectorSignal(const KVString& type)
{
   // Add a new signal to the list of detector's signals.
   // \param[in] type define the name of the signal to add
   // \returns pointer to the new signal object
   // \note do not `delete`{.cpp} the signal object: the detector handles deletion
   //
   // If the type is "SignalTrace", we add a special KVDetectorSignalTrace object
 
   if(type=="SignalTrace")
      return AddDetectorSignalTrace();
   auto signal = new KVDetectorSignal(type, this);
   fDetSignals.Add(signal);
   return signal;
}
 
 
 
 
KVDetectorSignalTrace* KVDetector::AddDetectorSignalTrace()
{
   // Add a new signal trace (with type "SignalTrace") to the list of detector's signals.
   //
   // \returns pointer to the new signal object
   // \note do not `delete`{.cpp} the signal object: the detector handles deletion
 
   auto signal = new KVDetectorSignalTrace(this);
   fDetSignals.Add(signal);
   return signal;
}
 
 
 
 
Bool_t KVDetector::AddDetectorSignalExpression(const KVString& type, const KVString& _expr)
{
   // Add a new KVDetectorSignalExpression to this detector
   //
   // \param[in] type the name/type of the new signal
   // \param[in] _expr mathematical expression using any of the known signals of the detector
   //
   // \note If the expression is not valid, no signal will be created and method returns kFALSE.
 
   KVDetectorSignalExpression* ds = new KVDetectorSignalExpression(type, _expr, this);
   if (!ds->IsValid()) {
      delete ds;
      ds = nullptr;
   }
   else
      AddDetectorSignal(ds);
   return (ds != nullptr);
}
 
 
 
 
void KVDetector::SetPressure(Double_t P)
{
   // \param[in] P pressure in [Torr]
   //
   // For a gaseous detector, set/change the pressure of the active gas layer.
   //
   // For ROOT geometries, we change the medium/material of the corresponding node in the geometry
   // in order to reflect the change in pressure (gases are represented by different media/materials
   // for each pressure/temperature) so that it will be taken into account for example when filtering simulated data.
 
   KVMaterial::SetPressure(P);
   if (ROOTGeo()) {
      // find node in geometry
      TGeoPhysicalNode* pn = (TGeoPhysicalNode*)gGeoManager->GetListOfPhysicalNodes()->FindObject(GetNode()->GetFullPathToNode());
      if (!pn) pn = gGeoManager->MakePhysicalNode(GetNode()->GetFullPathToNode());
      // use new medium reflecting change in pressure
      pn->GetVolume()->SetMedium(GetActiveLayer()->GetGeoMedium());
   }
}
 
 
 
 
void KVDetector::SetTemperature(Double_t T)
{
   // \param[in] T pressure in [deg. C]
   //
   // For a gaseous detector, set/change the temperature of the active gas layer.
   //
   // For ROOT geometries, we change the medium/material of the corresponding node in the geometry
   // in order to reflect the change in temperature (gases are represented by different media/materials
   // for each pressure/temperature) so that it will be taken into account for example when filtering simulated data.
 
   KVMaterial::SetTemperature(T);
   if (ROOTGeo()) {
      // find node in geometry
      TGeoPhysicalNode* pn = (TGeoPhysicalNode*)gGeoManager->GetListOfPhysicalNodes()->FindObject(GetNode()->GetFullPathToNode());
      if (!pn) pn = gGeoManager->MakePhysicalNode(GetNode()->GetFullPathToNode());
      // use new medium reflecting change in pressure
      pn->GetVolume()->SetMedium(GetActiveLayer()->GetGeoMedium());
   }
}