KVINDRAForwardGroupReconstructor.cpp

//Created by KVClassFactory on Tue Feb 27 11:35:10 2018
//Author: John Frankland,,,
 
#include "KVINDRAForwardGroupReconstructor.h"
#include "KVINDRACodeMask.h"
#include "KVIDGChIoSi.h"
#include "KVIDGCsI.h"
 
#include <KVCsI.h>
#include <KVSilicon.h>
 
ClassImp(KVINDRAForwardGroupReconstructor)
 
 
 
bool KVINDRAForwardGroupReconstructor::DoCoherencyAnalysis(KVReconstructedNucleus& PART)
{
   // Coherency analysis for forward rings 1-9 of INDRA
 
   bool ok = false;
   if (PART.GetStoppingDetector()->IsType("CSI")) {
      // particles stopping in CsI detectors
      // check coherency of CsI and Si-CsI identifications
      ok = CoherencySiCsI(PART);
      if(theChio)
      {
         // we check that the ChIo contribution is sane:
         // if no other particles hit this group, the Z given by the ChIoSi
         // must be <= the Z found from Si-CsI or CsI identification
         PART.SetParameter("UseFullChIoEnergyForCalib", CoherencyChIoSiCsI(PART));
      }
   }
   else {
      // particle stopped in Si (=> ChIo-Si, Si-PSA) or ChIo (=> Zmin)
      ok = PART.IsIdentified() && identifying_telescope;
      if(theChio)
         PART.SetParameter("UseFullChIoEnergyForCalib", !(theChio->GetNHits() > 1));
   }
 
   return ok;
}
 
 
 
 
void KVINDRAForwardGroupReconstructor::DoCalibration(KVReconstructedNucleus* PART)
{
   // Special calibration for particles in rings 1 to 9
   // We set the energy calibration code for the particle here:
   //
   //~~~{.cpp}
   //   KVINDRA::ECodes::NO_CALIBRATION_ATTEMPTED = 0, ///< particle stopped in detectors with no available calibration
   //   KVINDRA::ECodes::NORMAL_CALIBRATION = 1,       ///< normal well-calibrated particle with no problems
   //   KVINDRA::ECodes::SOME_ENERGY_LOSSES_CALCULATED = 2, ///< particle calibration OK, with some detector energies calculated
   //   KVINDRA::ECodes::WARNING_CSI_MAX_ENERGY = 3,    ///< particle calibration OK, although apparent energy would mean punching through the CsI
   //   KVINDRA::ECodes::BAD_CALIBRATION = 15          ///< calibration attempted but bad result (negative energies etc.)
   //~~~
   //
   // The contributions from ChIo, Si, and CsI are stored in particle parameters
   // INDRA.ECHIO, INDRA.ESI and INDRA.ECSI
   // If the contribution is calculated rather than measured (see below), it is stored as a negative value.
   // NOTE: in no case do we ever calculate an energy for uncalibrated detector using measured energy
   //       loss in the detector placed in front (i.e. calculate ECsI from deSi, or calculate ESi
   //       from deChIo) as this gives wildly varying (mostly false) results.
   //
   // For neutrons (IDCODE=1): if CsI is calibrated, we use the CsI light response to calculate
   //                          the equivalent energy for a proton.
   // For particles stopping in the CsI detector (IDCODE=2 or 3):
   //    - for "8Be" we use half the CsI light output to calculate 4He kinetic energy which is then doubled
   //      (we assume two 4He of identical energy), and ECODE=kECode2 (calculated)
   //    - if calibrated CsI energy is bizarre (<=0), we cannot calibrate the particle (ECODE=kECode15)
   //    - if the Silicon detector is not calibrated, or if coherency analysis indicates a pile-up
   //      in the silicon or other inconsistency, the Silicon energy is calculated from the calibrated
   //      CsI energy (ECODE=kECode2)
   //    - if the Ionisation Chamber is not calibrated, or if coherency analysis indicates a pile-up
   //      in the ChIo or other inconsistency, this contribution is calculated using the total
   //      cumulated energy measured in the CsI and/or Silicon detectors (ECODE=kECode2)
   //
   // For particles stopping in the Silicon detector (IDCODE=3, 4 or 5):
   //    - if no calibration exists for the Silicon, we cannot calibrate the particle (ECODE=kECode0)
   //    - if the Ionisation Chamber is not calibrated, or if coherency analysis indicates a pile-up
   //      in the ChIo or other inconsistency, this contribution is calculated using the total
   //      cumulated energy measured in the CsI and/or Silicon detectors (ECODE=kECode2)
 
   bool fUseFullChIoEnergyForCalib = PART->GetParameters()->GetBoolValue("UseFullChIoEnergyForCalib");
   bool fPileup = PART->GetParameters()->GetBoolValue("Pileup");
   bool fCoherent = PART->GetParameters()->GetBoolValue("Coherent");
 
   if (PART->GetIDCode() == KVINDRA::IDCodes::ID_NEUTRON) { // neutrons
      DoNeutronCalibration(PART);
      return;
   }
 
   Bool_t stopped_in_silicon = kTRUE;
 
   auto si = GetSi(PART);
   auto csi = GetCsI(PART);
   if (csi) {
      stopped_in_silicon = kFALSE;
      if (csi->IsCalibrated(Form("Z=%d,A=%d",PART->GetZ(),PART->GetA())) && csi->GetDetectorSignalValue("TotLight") > 0) {
         /* CSI ENERGY CALIBRATION */
         if (PART->GetIDCode() == KVINDRA::IDCodes::ID_CSI_PSA && PART->IsIsotope(4, 8)) {
            fECsI = DoBeryllium8Calibration(PART);
         }
         else
            fECsI = csi->GetCorrectedEnergy(PART, -1., kFALSE);
      }
      else if ((PART->GetIDCode()!=KVINDRA::IDCodes::ID_GAMMA) && si && !fPileup && fCoherent && si->IsCalibrated()) {
         // use silicon energy to calculate CsI energy
         fESi = si->GetCorrectedEnergy(PART);
         if (fESi <= 0) {
            // problem with silicon calibration
            SetBadCalibrationStatus(PART);
            return;
         }
         fECsI = si->GetEResFromDeltaE(PART->GetZ(), PART->GetA());
         // sanity check: very small DE can give very large Eres
         auto de_max = csi->GetMaxDeltaE(PART->GetZ(), PART->GetA());
         if (fECsI > de_max)
            SetCalibrationStatus(*PART, KVINDRA::WARNING_CSI_MAX_ENERGY);
         else
            SetCalibrationStatus(*PART, KVINDRA::SOME_ENERGY_LOSSES_CALCULATED);
         fECsI = -fECsI;
      }
      else {
         SetNoCalibrationStatus(PART);// no CsI calibration - no calibration performed
         return;
      }
      if (fECsI <= 0 && (PART->GetECode() != KVINDRA::ECodes::SOME_ENERGY_LOSSES_CALCULATED)) { // DoBeryllium8Calibration returns fECsI<0 with ECode=SOME_ENERGY_LOSSES_CALCULATED if OK
         SetBadCalibrationStatus(PART);// problem with CsI energy - no calibration
         return;
      }
   }
   if (si && (PART->GetIDCode()!=KVINDRA::IDCodes::ID_GAMMA)) {
      /* SILICIUM ENERGY CONTRIBUTION */
      // if fPileup = kTRUE, the Silicon energy appears to include a contribution from another particle
      //     therefore we have to estimate the silicon energy for this particle using the CsI energy
      // if fCoherent = kFALSE, the Silicon energy is too small to be consistent with the CsI identification,
      //     therefore we have to estimate the silicon energy for this particle using the CsI energy
 
      if (PART->GetIDCode() == KVINDRA::IDCodes::ID_STOPPED_IN_FIRST_STAGE && stopped_in_silicon) {
         // particles stopped in Silicon detectors with no identification except estimation of Zmin (case of INDRA without Ionization Chambers):
         // if silicon is calibrated, give measured energy in silicon as particle energy
         if (si->IsCalibrated())
            PART->SetEnergy(si->GetEnergy());
         else
            SetNoCalibrationStatus(PART);
         return;
      }
 
      if (!fPileup && fCoherent && si->IsCalibrated()) {
         // all is apparently well. use full energy deposited in silicon.
         Bool_t si_transmission = kTRUE;
         if (stopped_in_silicon) {
            si_transmission = kFALSE;
         }
         else {
            si->SetEResAfterDetector(TMath::Abs(fECsI));
         }
         fESi = si->GetCorrectedEnergy(PART, -1., si_transmission);
         if (fESi <= 0) {
            // problem with silicon calibration
            print_part = true;
            if (!stopped_in_silicon && (TMath::Abs(fECsI) > 0.0)) {
               if (!CalculateSiliconDEFromResidualEnergy(PART, si)) return;
            }
            else {
               // stopped in Si with bad Si calibration - no good
               SetBadCalibrationStatus(PART);
               return;
            }
         }
      }
      else {
         if (!stopped_in_silicon && (TMath::Abs(fECsI) > 0.0)) {
            if (!CalculateSiliconDEFromResidualEnergy(PART, si)) return;
         }
         else {
            // cannot calibrate particle stopping in silicon in this case
            SetNoCalibrationStatus(PART);
            return;
         }
      }
 
   }
   if (theChio && (PART->GetIDCode()!=KVINDRA::IDCodes::ID_GAMMA)) {
      /* CHIO ENERGY CONTRIBUTION */
      // if fUseFullChIoEnergyForCalib = kFALSE, we have to estimate the ChIo energy for this particle
      Double_t ERES = TMath::Abs(fESi) + TMath::Abs(fECsI);
      if (fUseFullChIoEnergyForCalib && theChio->IsCalibrated()) {
         // all is apparently well
         if (!PART->GetStoppingDetector()->IsType("CI")) {
            theChio->SetEResAfterDetector(ERES);
            fEChIo = theChio->GetCorrectedEnergy(PART);
         }
         else {
            // particle stopped in chio, not in transmission
            fEChIo = theChio->GetCorrectedEnergy(PART, -1., kFALSE);
         }
         if (fEChIo < 0.) {
            // bad chio calibration
            if (!PART->GetStoppingDetector()->IsType("CI") && (ERES > 0.0)) {
               if (!CalculateChIoDEFromResidualEnergy(PART, ERES)) return;
            }
         }
      }
      else {
         if (!PART->GetStoppingDetector()->IsType("CI") && (ERES > 0.0)) {
            if (!CalculateChIoDEFromResidualEnergy(PART, ERES)) return;
         }
         else {
            // particle stopped in ChIo, no calibration available
            SetNoCalibrationStatus(PART);
            return;
         }
      }
   }
   PART->SetEnergy(TMath::Abs(fECsI) + TMath::Abs(fESi) + TMath::Abs(fEChIo));
}
 
 
 
 
 
void KVINDRAForwardGroupReconstructor::DoNeutronCalibration(KVReconstructedNucleus* PART)
{
   // use energy of CsI calculated using the Z & A of the CsI identification
   // to calculate the energy deposited by the neutron
   KVIdentificationResult* IDcsi = PART->GetIdentificationResult(1);
   KVNucleus tmp(IDcsi->Z, IDcsi->A);
   auto csi = GetCsI(PART);
   if (csi && csi->IsCalibrated(Form("Z=%d,A=%d",IDcsi->Z,IDcsi->A))) {
      fECsI = csi->GetCorrectedEnergy(&tmp, -1., kFALSE);
      PART->SetEnergy(fECsI);
      SetCalibrationStatus(*PART, KVINDRA::ECodes::SOME_ENERGY_LOSSES_CALCULATED); // not a real energy measure
   }
   else
      SetNoCalibrationStatus(PART);
}
 
 
 
 
 
Bool_t KVINDRAForwardGroupReconstructor::CoherencyChIoSiCsI(KVReconstructedNucleus& PART)
{
   // Called by DoCoherency() for particles stopping in CsI detectors on rings 2-9,
   // which have a ChIo in front of them.
   //
   // we check that the ChIo contribution is sane:
   // if no other particles hit this group, the Z given by the ChIoSi
   // must be <= the Z found from Si-CsI or CsI-RL identification
   //
   // in this case the measured energy loss of the ChIo can be solely attributed to this particle
   // and we return kTRUE;
 
   // ChIo was hit by more than one particle in group
   if (theChio && theChio->GetNHits() > 1) {
      return kFALSE;
   }
 
   KVIdentificationResult* IDchiosi = PART.GetIdentificationResult(3);
   if (!IDchiosi->IDattempted) {
      // no ChIo-Si identification ? assume coherency ?
      return kTRUE;
   }
 
   // if we have a successful ChIo-Si id with a Z > Z given by CsI or Si-CsI id,
   // then we consider that there is a pile-up in the ChIo
   // note that we consider Z_ChIoSi = 1 to mean the ChIo is in the pedestal i.e. nothing seen
   // we also require the chio-si identification to be above the punch-through line
   if (IDchiosi->IDOK && IDchiosi->Z > 1 && IDchiosi->IDquality != KVIDGChIoSi::k_BelowPunchThrough && IDchiosi->Z > partID.Z) {
      return kFALSE;
   }
 
   return kTRUE;
}
 
 
 
 
Bool_t KVINDRAForwardGroupReconstructor::CoherencySiCsI(KVReconstructedNucleus& PART)
{
   // Called by DoCoherency() for particles stopping in CsI detectors on rings 1-9,
   // which have a Silicon detector just in front of them.
   //
   // coherent = kFALSE if CsI and Si-CsI identifications are not coherent.
   //                 the CsI identification is kept, either the Si signal
   //                  was not good (particle hitting dead zone), or it actually corresponds
   //                  to two particles reaching the CsI at the same time
   // pileup = kTRUE means that the (charged) particle identified in CsI-R/L is correct,
   //              and there is probably a second particle which stopped in the silicon
   //              detector at the same time, to be identified in ChIo-Si/Si-PSA after
   //              subtraction of the Silicon contribution
 
   KVIdentificationResult* IDcsi = PART.GetIdentificationResult(1);
   KVIdentificationResult* IDsicsi = PART.GetIdentificationResult(2);
   KVIDTelescope* idt_csi = (KVIDTelescope*)PART.GetReconstructionTrajectory()->GetIDTelescopes()->FindObjectByType(IDcsi->GetType());
   KVIDTelescope* idt_sicsi = (KVIDTelescope*)PART.GetReconstructionTrajectory()->GetIDTelescopes()->FindObjectByType(IDsicsi->GetType());
   // in case we have PSA identification in Si detector
   KVIdentificationResult* IDsi = PART.GetIdentificationResult("SI");
   KVIDTelescope* idt_si = (KVIDTelescope*)PART.GetReconstructionTrajectory()->GetIDTelescopes()->FindObjectByType("SI");
   auto the_CsI = idt_csi->GetDetector(1);
   bool fCoherent = kTRUE;
   bool fPileup = kFALSE;
   auto nid = PART.GetNumberOfIdentificationResults();
 
   // Unsuccessful/no CsI id attempt with successful id further up trajectory ?
   // Then use first successful identification result
   if (!IDcsi->IDOK) {
      if(nid>1)
      {
         for(int i_id=2; i_id<=nid; ++i_id)
         {
            auto next_idr = PART.GetIdentificationResult(i_id);
            if(next_idr && next_idr->IDattempted && next_idr->IDOK)
            {
               auto next_idt = (KVIDTelescope*)PART.GetReconstructionTrajectory()->GetIDTelescopes()->FindObjectByType(next_idr->GetIDType());
               if(next_idt)
               {
                  if(next_idt->HasDetector(the_CsI))
                  {
                     // CsI is still the stopping detector, we just change the ID
                     partID = *next_idr;
                     identifying_telescope = next_idt;
                     PART.SetParameter("Coherent", fCoherent);
                     PART.SetParameter("Pileup", fPileup);
                     PART.SetParameter("No CsI id", Form("Replaced with %s id", partID.GetIDType()));
                     return kTRUE;
                  }
                  else
                  {
                     // ID achieved without CSI, change stopping detector & recon trajectory
                     auto new_stop_det = next_idt->GetDetector(next_idt->GetSize());
                     auto traj = (KVGeoDNTrajectory*)new_stop_det->GetNode()->GetTrajectories()->First();
                     auto new_recon_traj = (const KVReconNucTrajectory*)GetGroup()->GetTrajectoryForReconstruction(traj, new_stop_det->GetNode());
                     PART.ModifyReconstructionTrajectory(new_recon_traj);
                     partID = *next_idr;
                     identifying_telescope = next_idt;
                     PART.SetParameter("Coherent", fCoherent);
                     PART.SetParameter("Pileup", fPileup);
                     PART.SetParameter("No CsI id", Form("Replaced with %s id, stopping detector changed to %s", partID.GetIDType(), new_stop_det->GetName()));
                     return kTRUE;
                  }
               }
            }
         }
      }
      // no identification possible
      if(PART.IsIdentified())
         PART.SetIsUnidentified();
      return kFALSE;
   }
 
   // check coherency of CsI-R/L and Si-CsI identifications
   if (IDcsi->IDOK) {
      partID = *IDcsi;
      identifying_telescope = idt_csi;
 
      /************ CHECKS FOR 'GAMMA' IDENTIFICATION IN CSI ID TELESCOPE ************/
      if(IDcsi->IDquality==KVIDGCsI::kICODE10)
      {
         if(nid>1)
         {
            for(int i_id=2; i_id<=nid; ++i_id)
            {
               auto next_idr = PART.GetIdentificationResult(i_id);
               if(next_idr && next_idr->IDattempted && next_idr->IDOK)
               {
                  auto next_idt = (KVIDTelescope*)PART.GetReconstructionTrajectory()->GetIDTelescopes()->FindObjectByType(next_idr->GetIDType());
                  if(next_idt)
                  {
                     if(next_idt->HasDetector(the_CsI))
                     {
                        // We check the next id telescope along the trajectory which includes the CsI detector as its second member
                        // If this telescope has a successful (charged) particle identification, we keep that one instead.
                        // -- check for incorrect gamma identification in CsI, with good charged particle id in DE-CsI
                        //    this may occur for low energy Z=1,2, etc. which are difficult to distinguish from gamma line
                        partID = *next_idr;
                        identifying_telescope = next_idt;
                        PART.SetParameter("Coherent", fCoherent);
                        PART.SetParameter("Pileup", fPileup);
                        PART.SetParameter("Bad gamma id", Form("switched to %s id: Z=%d A=%d",next_idr->GetIDType(),next_idr->Z,next_idr->A));
                        return kTRUE;
                     }
                     else
                     {
                        // a gamma identified in the CsI will hide any charged particles which stop in a detector further
                        // up the reconstruction trajectory. if an identification which does NOT include the CsI was successful,
                        // we add a coherency particle stopping in the 'last' detector of the corresponding idtelescope.
                        auto new_stop_det = next_idt->GetDetector(next_idt->GetSize());
                        coherency_particles.push_back(
                                 {
                                    &PART, new_stop_det->GetNode(), (KVGeoDNTrajectory*)the_CsI->GetNode()->GetTrajectories()->First(), next_idt, i_id, nid
                                 }
                                 );
                        fPileup=true;
                        PART.SetParameter("Coherent", fCoherent);
                        PART.SetParameter("Pileup", fPileup);
                        PART.SetParameter("Gamma pileup", Form("added new particle identified in %s", next_idt->GetName()));
                        return kTRUE;
                     }
                  }
               }
            }
         }
       }
      // Neutrons have no energy loss in Si detector (pedestal), have a reaction in the CsI and create a Z=1
      // or Z=2 which is identified in CsI R-L,  while they show up in Si-CsI maps as a horizontal
      // band around the Si pedestal for low energies (energies where proton dE is significantly larger than
      // the pedestal).
      if ((IDcsi->Z == 1 || IDcsi->Z == 2) && IDsicsi->IDattempted) {
         if (IDsicsi->deltaEpedestal == KVIdentificationResult::deltaEpedestal_YES) {
            // explicit treatment of 'neutron-like' particles with a cut in Si-CsI id grid
            partID = *IDsicsi;
            identifying_telescope = idt_sicsi;
            partID.IDOK = kTRUE;
            partID.Zident = kTRUE;
            partID.Aident = kTRUE;
            partID.Z = 0;
            partID.A = 1;
            partID.IDcode = KVINDRA::IDCodes::ID_NEUTRON; // general code for neutrons
            PART.SetParameter("Coherent", fCoherent);
            PART.SetParameter("Pileup", fPileup);
            return kTRUE;
         }
      }
 
      // We check the coherency of the mass and charge identifications
      //
      // If the CsI identification is not isotopic i.e. IDcsi->Aident=false and Si-CsI identification *is* isotopic,
      // we prefer the Si-CsI identification always.
      //
      // Otherwise, if CsI identification *is* isotopic:
      // If a successful Si-CsI identification is available we check:
      //   if Si-CsI has isotopic identification:
      //       - for Z=1, at low energies Si-CsI identification is more sure than PSA in CsI, as all isotopes merge
      //         together and can even be confused with gammas at low energy in the CsI.
      //         Therefore as long as Si-CsI ID is good for Z=1, we prefer Si-CsI
      //       - for Z>1, must have same Z, A within +/-1 unit
      //                                      if Z or A smaller => incoherency : pile-up of two particles in CsI ?
      //                                               or bad signal from Si detector (dead zones) ?
      //                                      if Z or A larger, CsI identification is good,
      //                                      assume another particle stopped in Si (identify in ChIo-Si)
      //   if Si-CsI gives just Z - must have same Z
      //                                      if Z smaller => incoherency : pile-up of two particles in CsI ?
      //                                               or bad signal from Si detector (dead zones) ?
      //                                      if Z larger => CsI identification is good,
      //                                      assume another particle stopped in Si (identify in ChIo-Si)
      //
      // If CsI identification gives code KVIDGCsI::kICode1 ou KVIDGCsI::kICode3 and the
      // Si-CsI gives the same Z but A = Acsi + 1, we use the Si-CsI identification.
      //
      // If CsI identification gives code KVIDGCsI::kICode2 ou KVIDGCsI::kICode3 and the
      // Si-CsI gives the same Z but A = Acsi - 1, we use the Si-CsI identification.
      //
      // N.B. if CsI-R/L identification gives "8Be" (2 alpha particles) then there are two correct possibilities:
      //     1) Si-CsI identification gives 7Li => CsI identification is correct ("8Be")
      //     2) Si-CsI identification gives 8He => the particle is 8He (falls on same R/L line as 2*alpha)
      // Seeing the very low statistics for 8He compared to 8Be/2alpha, we assume that if Si-CsI id
      // gives something above 8He it is either incoherent (below 7Li) or 8Be + something else in ChIo-Si
      // (above 7Li).
      Int_t Zref = IDcsi->Z;
      Int_t Aref = IDcsi->A;
      if (IDsicsi->IDOK) {
 
         if (IDsicsi->Aident) { // Si-CsI provides mass identification
 
            if(!IDcsi->Aident)
            {
               // CsI identification is approximative (Z only), Si-CsI has full Z & A identification
               // => use Si-CsI identification, don't ask questions about 'coherency'
               partID = *IDsicsi;
               identifying_telescope = idt_sicsi;
               PART.SetParameter("Coherent", fCoherent);
               PART.SetParameter("Pileup", fPileup);
               return kTRUE;
            }
 
            // for Z=1: at low energies Si-CsI identification is more sure than PSA in CsI, as all isotopes merge
            // together and can even be confused with gammas at low energy in the CsI.
            // Therefore as long as Si-CsI ID is good for Z=1, we prefer Si-CsI
            if(IDcsi->Z==1 && IDsicsi->Z==1)
            {
               partID = *IDsicsi;
               identifying_telescope = idt_sicsi;
               PART.SetParameter("Coherent", fCoherent);
               PART.SetParameter("Pileup", fPileup);
               return kTRUE;
            }
            if (IDcsi->Z == 4 && IDcsi->A == 8) {
               // traitement special 8Be
               // if sicsi => 7Li, it is 8Be (2alpha)
               // if sicsi => 8He, it is 8He
               if (IDsicsi->Z < 2 || (IDsicsi->Z == 2 && IDsicsi->A < 7)) {
                  fCoherent = kFALSE;
                  IDsicsi->SetComment("CsI-R/L & Si-CsI identifications not coherent");
                  PART.SetParameter("Coherent", fCoherent);
                  PART.SetParameter("Pileup", fPileup);
                  return kTRUE;
               }
               else if (IDsicsi->Z == 2 && IDsicsi->A > 6 && IDsicsi->A < 10) {
                  // accept helium-7,8,9 as 8He
                  partID = *IDsicsi;
                  identifying_telescope = idt_sicsi;
                  PART.SetParameter("Coherent", fCoherent);
                  PART.SetParameter("Pileup", fPileup);
                  return kTRUE;
               }
               else if ((IDsicsi->Z == 2 && IDsicsi->A > 9) || (IDsicsi->Z == 3 && IDsicsi->A < 6)) {
                  fCoherent = kFALSE;
                  IDsicsi->SetComment("CsI-R/L & Si-CsI identifications not coherent");
                  PART.SetParameter("Coherent", fCoherent);
                  PART.SetParameter("Pileup", fPileup);
                  return kTRUE;
               }
               else if (IDsicsi->Z == 3 && IDsicsi->A > 5 && IDsicsi->A < 9) {
                  // accept lithium-6,7,8 as 7Li
                  PART.SetParameter("Coherent", fCoherent);
                  PART.SetParameter("Pileup", fPileup);
                  return kTRUE;
               }
               else if ((IDsicsi->Z == 3 && IDsicsi->A > 8) || IDsicsi->Z > 3) {
                  fPileup = kTRUE;
                  IDsicsi->SetComment("Second particle stopping in Si, identification ChIo-Si required");
                  PART.SetParameter("Coherent", fCoherent);
                  PART.SetParameter("Pileup", fPileup);
                  return kTRUE;
               }
            }
            // if CsI says A could be bigger and Si-CsI gives same Z and A+1, use Si-CsI
            if ((IDsicsi->Z == Zref) && (IDsicsi->A == (Aref + 1))
                  && (IDcsi->IDquality == KVIDGCsI::kICODE1 || IDcsi->IDquality == KVIDGCsI::kICODE3)) {
               partID = *IDsicsi;
               identifying_telescope = idt_sicsi;
               PART.SetParameter("Coherent", fCoherent);
               PART.SetParameter("Pileup", fPileup);
               return kTRUE;
            }
            // if CsI says A could be smaller and Si-CsI gives same Z and A-1, use Si-CsI
            if ((IDsicsi->Z == Zref) && (IDsicsi->A == (Aref - 1))
                  && (IDcsi->IDquality == KVIDGCsI::kICODE2 || IDcsi->IDquality == KVIDGCsI::kICODE3)) {
               partID = *IDsicsi;
               identifying_telescope = idt_sicsi;
               PART.SetParameter("Coherent", fCoherent);
               PART.SetParameter("Pileup", fPileup);
               return kTRUE;
            }
            // Si-CsI gives unambiguous mass identification, CsI is less sure, and Si-CsI gives same Z, smaller A (=>no pileup): use Si-CsI
            if((IDsicsi->Z == Zref) && (IDsicsi->A < Aref)
               && (IDsicsi->IDquality == KVIDZAGrid::kICODE0 && IDcsi->IDquality > KVIDZAGrid::kICODE0))
              {
                partID = *IDsicsi;
                identifying_telescope = idt_sicsi;
                PART.SetParameter("Coherent", fCoherent);
                PART.SetParameter("Pileup", fPileup);
                return kTRUE;
              }
            // everything else - Z must be same, A +/- 1 unit
            if (IDsicsi->Z == Zref && TMath::Abs(IDsicsi->A - Aref) < 2) {
               PART.SetParameter("Coherent", fCoherent);
               PART.SetParameter("Pileup", fPileup);
               return kTRUE;
            }
            else if (IDsicsi->Z < Zref || (IDsicsi->Z == Zref && IDsicsi->A < (Aref - 1))) {
               fCoherent = kFALSE;
               IDsicsi->SetComment("CsI-R/L & Si-CsI identifications not coherent");
               PART.SetParameter("Coherent", fCoherent);
               PART.SetParameter("Pileup", fPileup);
               return kTRUE;
            }
            else  if (IDsicsi->Z > Zref || (IDsicsi->Z == Zref && IDsicsi->A > (Aref + 1))) {
               fPileup = kTRUE;
               IDsicsi->SetComment("Second particle stopping in Si");
               PART.SetParameter("Coherent", fCoherent);
               PART.SetParameter("Pileup", fPileup);
               return kTRUE;
            }
         }
         else {   // only Z identification from Si-CsI
            if (IDcsi->Z == 4 && IDcsi->A == 8) {
               // traitement special 8Be
               // we ask for Z to be equal 3 in SiCsI, but with no mass identification
               // we do not try for 8He identification
               if (IDsicsi->Z < 3) {
                  fCoherent = kFALSE;
                  IDsicsi->SetComment("CsI-R/L & Si-CsI identifications not coherent");
                  PART.SetParameter("Coherent", fCoherent);
                  PART.SetParameter("Pileup", fPileup);
                  return kTRUE;
               }
               else if (IDsicsi->Z == 3) {
                  PART.SetParameter("Coherent", fCoherent);
                  PART.SetParameter("Pileup", fPileup);
                  return kTRUE;
               }
               else {
                  fPileup = kTRUE;
                  IDsicsi->SetComment("Second particle stopping in Si, identification ChIo-Si required");
                  PART.SetParameter("Coherent", fCoherent);
                  PART.SetParameter("Pileup", fPileup);
                  return kTRUE;
               }
            }
            // everything else
            if (IDsicsi->Z == Zref) {
               PART.SetParameter("Coherent", fCoherent);
               PART.SetParameter("Pileup", fPileup);
               return kTRUE;
            }
            else if (IDsicsi->Z < Zref) {
               fCoherent = kFALSE;
               IDsicsi->SetComment("CsI-R/L & Si-CsI identifications not coherent");
               PART.SetParameter("Coherent", fCoherent);
               PART.SetParameter("Pileup", fPileup);
               return kTRUE;
            }
            else {
               fPileup = kTRUE;
 
               IDsicsi->SetComment("Second particle stopping in Si, identification ChIo-Si required");
               PART.SetParameter("Coherent", fCoherent);
               PART.SetParameter("Pileup", fPileup);
               return kTRUE;
            }
         }
      }
      else if (IDsicsi->IDattempted && IDsicsi->IDquality == KVIDZAGrid::kICODE8)
      {
         // Si-CsI id failed because totally out of grid i.e. no Si signal fired, but we have a good CsI id:
         // use the CsI id in this case
         partID = *IDcsi;
         identifying_telescope = idt_csi;
         PART.SetParameter("Coherent", fCoherent);
         PART.SetParameter("Pileup", fPileup);
         return kTRUE;
      }
      else if (IDsi && IDsi->IDOK)
      {
         // check for particle stopping in Si detector
         if(IDsi->Z > Zref)
         {
            // pile-up: a second charged particle stopped in Si
            auto new_stop_det = idt_si->GetDetector(1);
            coherency_particles.push_back(
                     {
                        &PART, new_stop_det->GetNode(), (KVGeoDNTrajectory*)the_CsI->GetNode()->GetTrajectories()->First(),
                        idt_si, 3, nid
                     }
                     );
            fPileup=true;
            PART.SetParameter("Coherent", fCoherent);
            PART.SetParameter("Pileup", fPileup);
            PART.SetParameter("Ion pileup", Form("added new particle identified in %s", idt_si->GetName()));
            return kTRUE;
         }
      }
      else if (IDsicsi->IDattempted)
      {
         if(IDsicsi->IDquality == KVIDZAGrid::kICODE4 || IDsicsi->IDquality == KVIDZAGrid::kICODE5)
         {
            // "noise" in between Si-CsI identification lines => pile-up with good particle
            // in CsI plus another particle stopped in the Silicon
            if (IDsicsi->Z < Zref) fCoherent = kFALSE;
            else fPileup = kTRUE;
         }
         else if (IDsicsi->IDquality == KVIDZAGrid::kICODE7)
         {
            // Point to identify in Si-CsI matrix is above the last line in the grid
            // Z given to identification result is Zmin = lastZ + 1
            // if this is > Z from CsI identification, pileup in Si
            if (IDsicsi->Z > Zref) fPileup = kTRUE;
         }
      }
      PART.SetParameter("Coherent", fCoherent);
      PART.SetParameter("Pileup", fPileup);
      return kTRUE;
   }
   PART.SetParameter("Coherent", fCoherent);
   PART.SetParameter("Pileup", fPileup);
   return kFALSE;
}
 
 
 
 
Bool_t KVINDRAForwardGroupReconstructor::CalculateSiliconDEFromResidualEnergy(KVReconstructedNucleus* n, KVDetector* si)
{
   // calculate fESi from fECsI
   // returns kTRUE if OK
   //
   // this sets directly the value of fESi
   //
   // if result is not good, bad calibration status is set for particle and we return false: in this case fESi=0
   //
   // if result is good, particle calibration status set to SOME_ENERGY_LOSSES_CALCULATED and fESi<0
 
   Double_t e0 = si->GetDeltaEFromERes(n->GetZ(), n->GetA(), TMath::Abs(fECsI));
   si->SetEResAfterDetector(TMath::Abs(fECsI));
   fESi = si->GetCorrectedEnergy(n, e0);
   if (fESi <= 0) {
      // can't calculate fESi from CsI energy - bad
      SetBadCalibrationStatus(n);
      fESi = 0;
      return kFALSE;
   }
   else {
      // calculated energy: negative
      fESi = -TMath::Abs(fESi);
      SetCalibrationStatus(*n, KVINDRA::ECodes::SOME_ENERGY_LOSSES_CALCULATED);
   }
   return kTRUE;
}