KVINDRAGroupReconstructor.cpp

//Created by KVClassFactory on Wed Feb 21 13:42:47 2018
//Author: John Frankland,,,
 
#include "KVINDRAGroupReconstructor.h"
#include "KVCalibratedSignal.h"
#include <KVIDGCsI.h>
#include <KVIDINDRACsI.h>
#include <KVINDRADetector.h>
 
ClassImp(KVINDRAGroupReconstructor)
 
TString KVINDRAGroupReconstructor::CSI_ID_TYPE = "CSI";
 
 
 
void KVINDRAGroupReconstructor::Identify()
{
   KVGroupReconstructor::Identify();
   for (auto& rnuc : ReconEventIterator(GetEventFragment()))
   {
      if (!rnuc.IsIdentified()) rnuc.SetIDCode(KVINDRA::IDCodes::NO_IDENTIFICATION); // unidentifiable particle
   }
}
 
 
 
 
void KVINDRAGroupReconstructor::IdentifyParticle(KVReconstructedNucleus& PART)
{
   // INDRA-specific particle identification.
   // Here we attribute the Veda6-style general identification codes depending on the
   // result of KVReconstructedNucleus::Identify and the subcodes from the different
   // identification algorithms:
   // If the particle's mass A was NOT measured, we make sure that it is calculated
   // from the measured Z using the mass formula defined by default
   //
   //IDENTIFIED PARTICLES
   //Identified particles with ID code = 2 with subcodes 4 & 5
   //(masse hors limite superieure/inferieure) are relabelled
   //with kIDCode10 (identification entre les lignes CsI)
   //
   //UNIDENTIFIED PARTICLES
   //Unidentified particles receive the general ID code for non-identified particles (kIDCode14)
   //EXCEPT if their identification in CsI gave subcodes 6 or 7
   //(Zmin) then they are relabelled "Identified" with IDcode = 9 (ident. incomplete dans CsI ou Phoswich (Z.min))
   //Their "identifying" telescope is set to the CsI ID telescope
 
   KVGroupReconstructor::IdentifyParticle(PART);
   // if a successful identification has occurred, identifying_telescope & partID are now set
   // and the particle has IsIdentified()=true
 
   // INDRA coherency treatment
   PART.SetParameter("Coherent", kTRUE);
   PART.SetParameter("Pileup", kFALSE);
   if(theChio)
      PART.SetParameter("UseFullChIoEnergyForCalib", kTRUE);
   Bool_t ok = DoCoherencyAnalysis(PART);
 
   if (ok) { // identification may change here due to coherency analysis
      PART.SetIsIdentified();
      PART.SetIdentification(&partID, identifying_telescope);
   } // if not ok, do we need to unset any previously identified particle?
 
   if (PART.IsIdentified()) {
 
      /******* IDENTIFIED PARTICLES *******/
      if (partID.IsType(CSI_ID_TYPE)) {     /**** CSI IDENTIFICATION ****/
 
         //Identified particles with ID code = 2 with subcodes 4 & 5
         //(masse hors limite superieure/inferieure) are relabelled
         //with kIDCode10 (identification entre les lignes CsI)
 
         Int_t grid_code = partID.IDquality;
         if (grid_code == KVIDGCsI::kICODE4 || grid_code == KVIDGCsI::kICODE5) {
            partID.IDcode = KVINDRA::IDCodes::ID_CSI_MASS_OUT_OF_RANGE;
            PART.SetIdentification(&partID, identifying_telescope);
         }
 
         // gamma particles identified as such...
         if(PART.GetZ()==0)
         {
            partID.IDcode = KVINDRA::IDCodes::ID_GAMMA;
            partID.Aident = true;
            partID.A = 0;
            PART.SetIdentification(&partID, identifying_telescope);
         }
 
      }
 
   }
   else {
 
      /******* UNIDENTIFIED PARTICLES *******/
 
      /*** general ID code for non-identified particles ***/
      PART.SetIDCode(KVINDRA::IDCodes::NO_IDENTIFICATION);
      auto csirl = id_by_type.find(CSI_ID_TYPE.Data());
      if (csirl != id_by_type.end()) {
         //Particles remaining unidentified are checked:
         //  - if in CsI identification they were in a zone marked 'nomassid|Z=...', we label them as identified
         //    in Z only (no mass identification), with the usual CsI IDCODE. This was introduced for E818 to handle
         //    the effects of defects in CsI crystals leading to poor discrimination of (p,d,t) above a certain energy/light
         //  - if their identification in CsI gave subcodes 6 or 7 (Zmin) then they are relabelled "Identified" with IDcode = 9
         //    (ident. incomplete dans CsI ou Phoswich (Z.min))
         //  - subcodes 4 & 5 (masse hors limite superieure/inferieure) are relabelled with kIDCode10
         //    (identification entre les lignes CsI)
         //In each case their "identifying" telescope is set to the CsI ID telescope
         KVIdentificationResult* idr = csirl->second;
         if (idr->IDattempted) {
            if (idr->IdentifyingGridHasFlagWhichBegins("nomassid|Z=")) {
               KVString flag = idr->IdentifyingGridGetFlagWhichBegins("nomassid|Z=");
               flag.Begin("=");
               flag.Next(); // split e.g. "nomassid|Z=1" into "nomassid|Z=" and "1"
               auto z = flag.Next().Atoi();
               PART.SetIsIdentified();
               idr->IDcode = KVINDRA::IDCodes::ID_CSI_PSA;
               idr->Z = z;
               idr->Zident = true;
               idr->Aident = false;
               idr->SetComment(Form("Z=%d identified without mass discrimination (due to CsI crystal defects)", z));
               partID = *idr;
               identifying_telescope = (KVIDTelescope*)PART.GetReconstructionTrajectory()->GetIDTelescopes()->FindObjectByType(CSI_ID_TYPE);
               PART.SetIdentification(&partID, identifying_telescope);
            }
            else if (idr->IDquality == KVIDGCsI::kICODE6 || idr->IDquality == KVIDGCsI::kICODE7) {
               PART.SetIsIdentified();
               idr->IDcode = KVINDRA::IDCodes::ID_CSI_FRAGMENT;
               partID = *idr;
               identifying_telescope = (KVIDTelescope*)PART.GetReconstructionTrajectory()->GetIDTelescopes()->FindObjectByType(CSI_ID_TYPE);
               PART.SetIdentification(&partID, identifying_telescope);
            }
            else if (idr->IDquality == KVIDGCsI::kICODE4 || idr->IDquality == KVIDGCsI::kICODE5) {
               PART.SetIsIdentified();
               idr->IDcode = KVINDRA::IDCodes::ID_CSI_MASS_OUT_OF_RANGE;
               partID = *idr;
               identifying_telescope = (KVIDTelescope*)PART.GetReconstructionTrajectory()->GetIDTelescopes()->FindObjectByType(CSI_ID_TYPE);
               PART.SetIdentification(&partID, identifying_telescope);
            }
         }
      }
 
   }
}
 
 
 
 
Bool_t KVINDRAGroupReconstructor::CalculateChIoDEFromResidualEnergy(KVReconstructedNucleus* n, Double_t ERES)
{
   // calculate fEChIo from residual energy
   //
   // returns kFALSE if it doesn't work, and sets particle bad calibration status
   //
   // returns kTRUE if it works, and sets calib status to SOME_ENERGY_LOSSES_CALCULATED
 
   Double_t e0 = theChio->GetDeltaEFromERes(n->GetZ(), n->GetA(), ERES);
   theChio->SetEResAfterDetector(ERES);
   fEChIo = theChio->GetCorrectedEnergy(n, e0);
   if (fEChIo <= 0) {
      SetBadCalibrationStatus(n);
      fEChIo = 0;
      return kFALSE;
   }
   fEChIo = -TMath::Abs(fEChIo);
   n->SetECode(KVINDRA::ECodes::SOME_ENERGY_LOSSES_CALCULATED);
   return kTRUE;
}
 
 
 
 
void KVINDRAGroupReconstructor::CalibrateCoherencyParticle(KVReconstructedNucleus* nuc)
{
   CalibrateParticle(nuc);
}
 
 
 
 
void KVINDRAGroupReconstructor::CalibrateParticle(KVReconstructedNucleus* PART)
{
   // Calculate and set the energy of a (previously identified) reconstructed particle
   //
   // This is only possible for correctly identified particles.
   // We exclude IDCODE9 particles (Zmin in CsI-RL)
 
   fEChIo = fESi = fECsI = 0;
 
   print_part = false;
 
   SetNoCalibrationStatus(PART);
 
   if (PART->GetIDCode() != KVINDRA::IDCodes::ID_CSI_FRAGMENT && PART->GetIDCode() != KVINDRA::IDCodes::ID_CSI_MASS_OUT_OF_RANGE) {
      // this status may be modified depending on what happens in DoCalibration
      SetCalibrationStatus(*PART, KVINDRA::ECodes::NORMAL_CALIBRATION);
      DoCalibration(PART);
   }
 
   if(theChio)
      PART->SetParameter("INDRA.ECHIO", fEChIo);
   PART->SetParameter("INDRA.ESI", fESi);
   PART->SetParameter("INDRA.ECSI", fECsI);
   //add correction for target energy loss - moving charged particles only!
   Double_t E_targ = 0.;
   if (PART->GetZ() && PART->GetEnergy() > 0) {
      E_targ = GetTargetEnergyLossCorrection(PART);
      PART->SetTargetEnergyLoss(E_targ);
   }
   Double_t E_tot = PART->GetEnergy() + E_targ;
   PART->SetEnergy(E_tot);
   // set particle momentum from telescope dimensions (random)
   PART->GetAnglesFromReconstructionTrajectory();
   CheckCsIEnergy(PART);
 
   //if(print_part) PART->Print();
}
 
 
 
 
double KVINDRAGroupReconstructor::DoBeryllium8Calibration(KVReconstructedNucleus* n)
{
   // Beryllium-8 = 2 alpha particles of same energy
   // We halve the total light output of the CsI to calculate the energy of 1 alpha
   // Then multiply resulting energy by 2
   // Note: fECsI is -ve, because energy is calculated not measured
 
   auto csi = GetCsI(n);
   // energy is usually calculated from "TotLight", but it could be "ScaledTotLight", or something else...
   auto input_signal = ((KVCalibratedSignal*)csi->GetDetectorSignal("Energy"))->GetInputSignal();
   Double_t half_light = input_signal->GetValue(Form("Z=%d,A=%d", n->GetZ(), n->GetA())) * 0.5;
   KVNucleus tmp(2, 4);
   double ecsi = 2.*csi->GetCorrectedEnergy(&tmp, half_light, kFALSE);
   if (ecsi > 0) {
      SetCalibrationStatus(*n, KVINDRA::ECodes::SOME_ENERGY_LOSSES_CALCULATED);
      // calculated energy returned as negative value
      return -ecsi;
   }
   else SetBadCalibrationStatus(n);
   return 0;
}
 
 
 
 
void KVINDRAGroupReconstructor::CheckCsIEnergy(KVReconstructedNucleus* n)
{
   // Check calculated CsI energy loss of particle:
   //   - If the particle was identified by the CsI but only in Z (case of particles falling in zones labelled "nomassid|Z=...")
   //     then we give it the (bad) calibration code KVINDRA::ECodes::WARNING_CSI_NO_MASS_ID
   //   - If it is greater than the maximum theoretical energy loss (depending on the length of CsI, the Z & A of the particle)
   //     we set the energy calibration code to KVINDRA::ECodes::WARNING_CSI_MAX_ENERGY (=3, historical VEDA code for particles with E_csi > E_max_csi)
 
   if (n->GetIDCode() == KVINDRA::IDCodes::ID_CSI_PSA && !n->IsAMeasured()) {
      SetCalibrationStatus(*n, KVINDRA::ECodes::WARNING_CSI_NO_MASS_ID);
      return;
   }
   KVDetector* csi = GetCsI(n);
   if (csi && (n->GetZ() > 0) && (n->GetZ() < 3) && (csi->GetDetectorSignalValue("Energy", Form("Z=%d,A=%d", n->GetZ(), n->GetA())) > csi->GetMaxDeltaE(n->GetZ(), n->GetA()))) {
      SetCalibrationStatus(*n, KVINDRA::ECodes::WARNING_CSI_MAX_ENERGY);
   }
}
 
 
 
 
void KVINDRAGroupReconstructor::AddCoherencyParticles()
{
   // Called to add any nuclei not included in the initial reconstruction, but "revealed"
   // by consistency checks between identifications and calibrations of other nuclei
   //
   // These particles have a (string) parameter "COHERENCY"
 
   for (auto& part : coherency_particles) {
 
      auto rnuc = GetEventFragment()->AddParticle();
 
      // reconstruction
      auto Rtraj = (const KVReconNucTrajectory*)GetGroup()->GetTrajectoryForReconstruction(part.stopping_trajectory, part.stopping_detector_node);
 
      rnuc->SetReconstructionTrajectory(Rtraj);
      rnuc->SetParameter("ARRAY", "INDRA");
      rnuc->SetParameter("COHERENCY",
                         "Particle added to event after consistency checks between identifications and calibrations of other nuclei");
      // for coherency particles in pileup with a charged particle (original particle has "Ion pileup" parameter),
      // when possible, we subtract energy losses of original particle from detectors on this particle's trajectory
      // this is only done if all detectors are calibrated and have invertible calibrations
      bool ok = false;
      if(part.original_particle->GetParameters()->HasParameter("Ion pileup"))
      {
         ok = true;
         Rtraj->IterateFrom();
         KVGeoDetectorNode* N;
         while( (N = Rtraj->GetNextNode() ))
         {
            auto d = N->GetDetector();
            if(!d->IsCalibrated())
            {
               ok = false;
               break;
            }
         }
         if(ok)
         {
            Rtraj->IterateFrom();
            while( (N = Rtraj->GetNextNode() ))
            {
               auto d = N->GetDetector();
               TString eparam = Form("INDRA.E%s",d->GetType());
               if(part.original_particle->GetParameters()->HasParameter(eparam))
               {
                  auto de = TMath::Abs(part.original_particle->GetParameters()->GetDoubleValue(eparam));
                  auto cal = static_cast<KVCalibratedSignal*>(d->GetDetectorSignal("Energy"));
                  auto e_orig = cal->GetValue();
                  if(cal->GetParentSignal() && de < e_orig)
                  {
                     // store initial parameter values in original particle
                     part.original_particle->SetParameter(Form("%s.%s",d->GetType(),cal->GetParentSignal()->GetType()), cal->GetParentSignal()->GetValue());
                     ok = cal->SetRawValueFromCalibrated(e_orig - de);
                  }
                  else
                     ok = false;
                  if(!ok)
                     break;
               }
            }
         }
      }
      if(ok) // for Ion pileup, if all went well, now (re)identify the particle
      {
         IdentifyParticle(*rnuc);
      }
      else // just reuse the identifications from the original particle
      {
         Int_t idnumber = 1;
         for (int i = part.first_id_result_to_copy; i <= part.max_id_result_index; ++i) {
            auto IDR = rnuc->GetIdentificationResult(idnumber++);
            // copy existing identification results from "parent" particle
            part.original_particle->GetIdentificationResult(i)->Copy(*IDR);
         }
         rnuc->SetIsIdentified();
         rnuc->SetIdentification(rnuc->GetIdentificationResult(1), part.identifying_telescope);
         rnuc->SetParameter("Pileup",false);
         rnuc->SetParameter("Coherent",true);
      }
      CalibrateCoherencyParticle(rnuc);
   }
}
 
 
/*
   The following is a copy of the old
void KVINDRAReconEvent::SecondaryAnalyseGroup(KVGroup* grp)
   method from branch 1.10
 
void KVINDRAGroupReconstructor:PerformSecondaryAnalysis()
{
   // Perform identifications and calibrations of particles not included
   // in first round (methods IdentifyEvent() and CalibrateEvent()).
   //
   // Here we treat particles with GetStatus()==KVReconstructedNucleus::kStatusOKafterSub
   // after subtracting the energy losses of all previously calibrated particles in group from the
   // measured energy losses in the detectors they crossed.
 
   // loop over al identified & calibrated particles in group and subtract calculated
   // energy losses from all detectors
   KVINDRAReconNuc* nuc;
   TList sixparts;
   TIter parts(grp->GetParticles());
   while ((nuc = (KVINDRAReconNuc*)parts())) {
      if (nuc->IsIdentified() && nuc->IsCalibrated()) {
         nuc->SubtractEnergyFromAllDetectors();
         // reconstruct particles from pile-up in silicon detectors revealed by coherency CsIR/L - SiCsI
         if (nuc->IsSiPileup() && nuc->GetSi()->GetEnergy() > 0.1) {
            KVINDRAReconNuc* SIX = AddParticle();
            SIX->Reconstruct(nuc->GetSi());
            sixparts.Add(SIX);
         }
         // reconstruct particles from pile-up in si75 detectors revealed by coherency
         if (nuc->IsSi75Pileup()) {
            KVINDRAReconNuc* SIX = AddParticle();
            SIX->Reconstruct(nuc->GetSi75());
            sixparts.Add(SIX);
         }
         // reconstruct particles from pile-up in sili detectors revealed by coherency
         if (nuc->IsSiLiPileup()) {
            KVINDRAReconNuc* SIX = AddParticle();
            SIX->Reconstruct(nuc->GetSiLi());
            sixparts.Add(SIX);
         }
 
         // reconstruct particles from pile-up in ChIo detectors revealed by coherency CsIR/L - ChIoCsI
         if (nuc->IsChIoPileup() && nuc->GetChIo()->GetEnergy() > 1.0) {
            KVINDRAReconNuc* SIX = AddParticle();
            SIX->Reconstruct(nuc->GetChIo());
            sixparts.Add(SIX);
         }
      }
   }
   // reanalyse group
   KVReconstructedNucleus::AnalyseParticlesInGroup(grp);
 
   Int_t nident = 0; //number of particles identified in each step
   if (sixparts.GetEntries()) { // identify any particles added by coherency CsIR/L - SiCsI
      KVINDRAReconNuc* SIX;
      TIter nextsix(&sixparts);
      while ((SIX = (KVINDRAReconNuc*)nextsix())) {
         if (SIX->GetStatus() == KVReconstructedNucleus::kStatusOK) {
            SIX->Identify();
            if (SIX->IsIdentified()) {
               nident++;
               if (SIX->GetCodes().TestIDCode(kIDCode5)) SIX->SetIDCode(kIDCode7);
               else SIX->SetIDCode(kIDCode6);
               SIX->Calibrate();
               if (SIX->IsCalibrated()) SIX->SubtractEnergyFromAllDetectors();
            } else {
               // failure of ChIo-Si identification: particle stopped in ChIo ?
               // estimation of Z (minimum) from energy loss (if detector is calibrated)
               UInt_t zmin = ((KVDetector*)SIX->GetDetectorList()->Last())->FindZmin(-1., SIX->GetMassFormula());
               if (zmin) {
                  SIX->SetZ(zmin);
                  SIX->SetIsIdentified();
                  SIX->SetIDCode(kIDCode7);
                  // "Identifying" telescope is taken from list of ID telescopes
                  // to which stopping detector belongs
                  SIX->SetIdentifyingTelescope((KVIDTelescope*)SIX->GetStoppingDetector()->GetIDTelescopes()->Last());
                  SIX->Calibrate();
               }
            }
         }
      }
   }
   if (nident) { // newly-identified particles may change status of others in group
      // reanalyse group
      KVReconstructedNucleus::AnalyseParticlesInGroup(grp);
      nident = 0;
   }
 
   TIter parts2(grp->GetParticles()); // list may have changed if we have added particles
   // identify & calibrate any remaining particles with status=KVReconstructedNucleus::kStatusOK
   while ((nuc = (KVINDRAReconNuc*)parts2())) {
      if (!nuc->IsIdentified() && nuc->GetStatus() == KVReconstructedNucleus::kStatusOK) {
         nuc->ResetNSegDet();
         nuc->Identify();
         if (nuc->IsIdentified()) {
            nident++;
            nuc->Calibrate();
            if (nuc->IsCalibrated()) nuc->SubtractEnergyFromAllDetectors();
         }
      }
   }
   if (nident) { // newly-identified particles may change status of others in group
      // reanalyse group
      KVReconstructedNucleus::AnalyseParticlesInGroup(grp);
      nident = 0;
   }
 
   // any kStatusOKafterShare particles ?
   TList shareChIo;
   parts2.Reset();
   while ((nuc = (KVINDRAReconNuc*)parts2())) {
      if (!nuc->IsIdentified() && nuc->GetStatus() == KVReconstructedNucleus::kStatusOKafterShare) {
         shareChIo.Add(nuc);
      }
   }
   Int_t nshares = shareChIo.GetEntries();
   if (nshares) {
      KVChIo* theChIo = ((KVINDRAReconNuc*)shareChIo.At(0))->GetChIo();
      if (theChIo && nshares > 1) {
         // divide chio energy equally
         Double_t Eshare = theChIo->GetEnergyLoss() / nshares;
         theChIo->SetEnergyLoss(Eshare);
         // modify PG and GG of ChIo according to new energy loss
         Double_t volts = theChIo->GetVoltsFromEnergy(Eshare);
         Double_t GG = theChIo->GetCanalGGFromVolts(volts);
         Double_t PG = theChIo->GetCanalPGFromVolts(volts);
         theChIo->GetACQParam("PG")->SetData(TMath::Min(4095, (Int_t)PG));
         theChIo->GetACQParam("GG")->SetData(TMath::Min(4095, (Int_t)GG));
      }
      // now try to identify
      TIter nextSh(&shareChIo);
      while ((nuc = (KVINDRAReconNuc*)nextSh())) {
         nuc->SetNSegDet(10);
         nuc->Identify();
         if (nuc->IsIdentified()) {
            nuc->SetIDCode(kIDCode8);
            nuc->Calibrate();
         }
      }
   }
 
   // any remaining stopped in first stage particles ?
   parts2.Reset();
   while ((nuc = (KVINDRAReconNuc*)parts2())) {
      if (!nuc->IsIdentified() && nuc->GetStatus() == KVReconstructedNucleus::kStatusStopFirstStage) {
         // estimation of Z (minimum) from energy loss (if detector is calibrated)
         UInt_t zmin = ((KVDetector*)nuc->GetDetectorList()->Last())->FindZmin(-1., nuc->GetMassFormula());
         if (zmin) {
            nuc->SetZ(zmin);
            nuc->SetIsIdentified();
            nuc->SetIDCode(kIDCode5);
            // "Identifying" telescope is taken from list of ID telescopes
            // to which stopping detector belongs
            nuc->SetIdentifyingTelescope((KVIDTelescope*)nuc->GetStoppingDetector()->GetIDTelescopes()->Last());
            nuc->Calibrate();
         }
      }
   }
}
*/