KVFilterGroupReconstructor.cpp

#include "KVFilterGroupReconstructor.h"
#include "KVIDZAGrid.h"
#include "KVMultiDetArray.h"
#include "KVNucleusEvent.h"
 
std::unordered_map<std::string,bool> KVFilterGroupReconstructor::audit_id_extensions;
 
 
 
KVReconNucTrajectory* KVFilterGroupReconstructor::get_recon_traj_for_particle(const KVGeoDNTrajectory* traj, const KVGeoDetectorNode* node)
{
   // We look for a simulated particle stopped in the same node with the same trajectory.
   // We use the actual trajectory of the simulated particle to get the right reconstruction trajectory
 
   //Info("get_recon_traj","traj=%s node=%s",traj->GetPathString().Data(), node->GetName());
   KVReconNucTrajectory* Rtraj{nullptr};
   for (auto& n : EventIterator(fSimEvent.get())) {
      TString stop = n.GetParameters()->GetTStringValue("STOPPING DETECTOR");
      //Info("get_recon_traj", "stop=%s", stop.Data());
      if (stop == node->GetName()) {
         Rtraj = (KVReconNucTrajectory*)GetGroup()->GetTrajectoryForReconstruction(traj, node);
         if (Rtraj) {
            //Info("get_recon_traj","Rtraj=%s TRAJECTORY=%s",Rtraj->GetPathString().Data(),n.GetParameters()->GetStringValue("TRAJECTORY"));
            auto traj_string = n.GetParameters()->GetTStringValue("TRAJECTORY");
            auto traj_in_grp = (KVGeoDNTrajectory*)GetGroup()->GetTrajectories()->FindObject(traj_string);
            if ((Rtraj->GetPathString() == traj_string) || (traj_in_grp && Rtraj->ContainsTrajectory(traj_in_grp))) {
               // correspondence recon <-> simu
               //Info("get_recon_traj","sim particle added to part_correspond");
               part_correspond[current_nuc_recon].Add(&n);
               // copy simulated particle parameters to reconstructed particle
               n.GetParameters()->Copy(*current_nuc_recon->GetParameters());
            }
         }
      }
   }
   if (!Rtraj) {
      Info("get_recon_traj_for_particle", "%p : noRtraj for this: %s %s", this, node->GetName(), traj->GetPathString().Data());
      if (!fSimEvent->GetMult()) {
         Info("get_recon_traj_for_particle", "%p : sim event is empty", this);
      }
   }
   return Rtraj;
}
 
 
 
 
void KVFilterGroupReconstructor::identify_particle(KVIDTelescope* idt, KVIdentificationResult* IDR, KVReconstructedNucleus& nuc)
{
   // check associated simulated particle passes identification threshold
   //
   // if more than one simulated particle stopped in the same detector, we add the Z, A and E of each particle together
   // and merge the lists of parameters into one.
 
   int multihit;
   if ((multihit = part_correspond[&nuc].GetEntries()) > 1) {
      KVNucleus sum_part;
      TIter next(&part_correspond[&nuc]);
      KVNucleus* N;
      KVNameValueList q;
      // the parameter lists are added to q because the KVNucleus operator+ does
      // not retain any parameters from either of its arguments
      while ((N = (KVNucleus*)next())) {
         sum_part += *N;
         q.Merge(*N->GetParameters());
      }
      sum_part.GetParameters()->Merge(q);
      auto sim_part = (KVNucleus*)part_correspond[&nuc].First();
      // replace properties of first nucleus in list with summed properties
      part_correspond[&nuc].Clear();
      sum_part.Copy(*sim_part);
      part_correspond[&nuc].Add(sim_part);
      sim_part->SetParameter("PILEUP", Form("%d particles in %s", multihit, nuc.GetStoppingDetector()->GetName()));
      // check 'string' parameters (TRAJECTORY,STOPPING DETECTOR,DETECTED)
      // 'merging' string parameters forms a comma-separated list of individual values
      // these should be replaced with a single value - unless they are not all the same!?!
      for (auto& param : *sim_part->GetParameters()) {
         KVString unique_value;
         if (param.IsString()) {
            KVString sparam(param.GetString());
            sparam.Begin(",");
            auto next = sparam.Next();
            if (unique_value.IsNull())
               unique_value = next;
            else {
               if (next != unique_value) {
                  KVError::Warning(this, "identify_particle", "different values of %s for %d-particle pileup in %s:%s, %s, ...",
                          param.GetName(), multihit, nuc.GetStoppingDetector()->GetName(), unique_value.Data(), next.Data());
               }
            }
            param.Set(unique_value);
         }
      }
      nuc.GetParameters()->Concatenate(*sim_part->GetParameters());
   }
   auto sim_part = (KVNucleus*)part_correspond[&nuc].First();
 
   IDR->SetIDType(idt->GetType());
   IDR->IDattempted = true;
 
   // why does this happen???
   if (!sim_part) {
      KVError::Warning(this, "identify_particle", "sim_part=nullptr; multihit=%d", multihit);
      KVError::Warning(this, "identify_particle", "sim_event contains:");
      fSimEvent->Print();
      KVError::Warning(this, "identify_particle", "recon_nuc stopping detector is: %s", nuc.GetStoppingDetector()->GetName());
      nuc.Print();
      KVError::Error(this, "identify_particle", "no sim_part for recon_nuc!");
      IDR->IDOK=false;
      return;
   }
 
   check_identification_threshold(sim_part, idt, IDR);
 
   // particles which punch-through ("DETECTED"="PUNCH THROUGH")
   //   for DE-E identification (identifying telescope has 2 members) particles which punch-through the
   //   the 2nd detector will not be identified correctly - probably rejected by identification grid.
   //    => if "RESIDUAL ENERGY" parameter is >2% of total energy of particle, reject as unidentified
   //
   //   for PSA identification in silicon (identifying telescope has 1 member of type "Si" or "SI")
   //   particles which punch-through the the detector will not be identified correctly
   //   - probably rejected by identification grid.
   //    => if "RESIDUAL ENERGY" parameter is >2% of total energy of particle, reject as unidentified
   //
   //   light particles which punch through CsI detectors are estimated to be identifiable,
   //   and will have a lower reconstructed energy than the real value
   if (sim_part->GetParameters()->IsValue("DETECTED", "PUNCH THROUGH"))
   {
      if(idt->GetSize()==2 || (idt->GetSize()==1 && idt->GetDetector(1)->IsType("Si",true)))
      {
         if(sim_part->GetParameters()->GetDoubleValue("RESIDUAL ENERGY")>0.02*sim_part->GetE())
         {
            IDR->IDOK = false;
            IDR->IDcode = GetGroup()->GetParentStructure<KVMultiDetArray>()->GetBadIDCode();
            IDR->IDquality = KVIDZAGrid::kICODE8;
            IDR->SetComment("punch through particle, identification not possible/good");
         }
      }
   }
}
 
 
 
 
void KVFilterGroupReconstructor::check_identification_threshold(KVNucleus* nuc, KVIDTelescope* idt, KVIdentificationResult* IDR)
{
   // Check that identification is possible for the given nucleus with its Z, A and (total) kinetic energy.
   // The nucleus in question may be from the simulation (KVSimNucleus) or from the reconstruction (KVReconstructedNucleus).
   // The KVIdentificationResult (identification status of the nucleus) may be changed by this method:
   // a particle previously thought identified may become unidentified here, as the reconstructed energy is
   // lower than the simulated energy, and therefore below the experimental threshold, although the "real"
   // energy of the simulated nucleus permitted its identification. If we don't do this, we end up with the
   // inconsistent final result of identified particles with energies below the experimental identification threshold.
 
   if (fDataQAudit) {
      // use experimental thresholds etc. to decide what can be identified
      int A = nuc->GetA();
      auto status = fDataQAudit->CanIdentify(get_id_telescope_name_for_audit(idt), nuc->GetZ(), A, nuc->GetEnergy());
 
      // if no Z as high as the one of this particle is included in the audit, and if it has been
      // authorized for this kind of telescope, we will use the old method of theoretical threshold
      // to determine whether it can be identified or not (see below)
      if(!((status == KVDataQualityAudit::IDStatus::EStatus::identification_not_possible_for_this_Z)
            && can_extend_identification_beyond_audit(idt->GetLabel())))
      {
         std::stringstream comments;
         comments << status;
         if (status) {
            IDR->IDOK = true;
            IDR->Z = nuc->GetZ();
            IDR->Zident = true;
            IDR->A = A; // may have been changed in CanIdentify()
            IDR->IDcode = idt->GetIDCode();
            IDR->Aident = false;
            if (status == KVDataQualityAudit::IDStatus::EStatus::full_isotopic_identification
                || status == KVDataQualityAudit::IDStatus::EStatus::Z_identified_but_not_A)
               IDR->IDquality = KVIDZAGrid::kICODE0;
            else if (status == KVDataQualityAudit::IDStatus::EStatus::partial_isotopic_identification) {
               IDR->IDquality = KVIDZAGrid::kICODE1;
               IDR->Aident = true;
            }
            else if (status == KVDataQualityAudit::IDStatus::EStatus::A_never_measured_but_between_Amin_and_Amax) {
               IDR->IDquality = KVIDZAGrid::kICODE3;
               IDR->Aident = true;
            }
            if (status == KVDataQualityAudit::IDStatus::EStatus::full_isotopic_identification)
               IDR->Aident = true;
         }
         else {
            IDR->IDOK = false;
            if (status == KVDataQualityAudit::IDStatus::EStatus::A_greater_than_measured_Amax)
               IDR->IDquality = KVIDZAGrid::kICODE5;
            else if (status == KVDataQualityAudit::IDStatus::EStatus::A_less_than_measured_Amin)
               IDR->IDquality = KVIDZAGrid::kICODE4;
            else if (status == KVDataQualityAudit::IDStatus::EStatus::weird_result_1
                     || status == KVDataQualityAudit::IDStatus::EStatus::weird_result_2
                     || status == KVDataQualityAudit::IDStatus::EStatus::weird_result_3
                     || status == KVDataQualityAudit::IDStatus::EStatus::weird_result_4
                     )
               KVError::Warning(this, "identify_particle",
                       "%s : %s Z=%d A=%d E=%f", comments.str().c_str(),
                       idt->GetName(), nuc->GetZ(), nuc->GetA(), nuc->GetEnergy());
         }
         IDR->SetComment(comments.str().c_str());
         return;
      }
   }
 
   // either no data quality audit is used, or no Z as high as the one of this particle is included in the audit,
   // and it is authorized for this kind of telescope to try to extend the identification beyond the limit of the audit
   if (idt->CanIdentify(nuc->GetZ(), nuc->GetA())) {
      IDR->Z = nuc->GetZ();
      IDR->A = nuc->GetA();
      if (idt->CheckTheoreticalIdentificationThreshold(nuc)) {
         IDR->IDOK = true;
         IDR->IDcode = idt->GetIDCode();
         IDR->IDquality = KVIDZAGrid::kICODE0;
         // set mass identification status depending on telescope, Z, A & energy of simulated particle
         idt->SetIdentificationStatus(IDR, nuc);
      }
      else {
         IDR->IDOK = false;
         IDR->IDcode = GetGroup()->GetParentStructure<KVMultiDetArray>()->GetIDCodeForParticlesStoppingInFirstStageOfTelescopes();
         IDR->IDquality = KVIDZAGrid::kICODE8;
         IDR->SetComment("below threshold for identification");
      }
   }
   else {
      IDR->IDOK = false;
   }
}
 
 
 
 
void KVFilterGroupReconstructor::PerformSecondaryAnalysis()
{
   // After first round of identification in group, try to identify remaining particles
 
   int num_ident_0, num_ident;
   int num_unident_0, num_unident;
   do {
      AnalyseParticles();
      num_ident = num_ident_0 = GetNIdentifiedInGroup();
      num_unident_0 = GetNUnidentifiedInGroup();
      Reconstruct();
      num_unident = GetNUnidentifiedInGroup();
      if (num_unident > num_unident_0) {
         Identify();
         num_ident = GetNIdentifiedInGroup();
         if (num_ident > num_ident_0) Calibrate();
      }
   }
   while (num_ident > num_ident_0); // continue until no more new identifications occur
}
 
 
 
 
 
 
void KVFilterGroupReconstructor::ReconstructParticle(KVReconstructedNucleus* part, const KVGeoDNTrajectory* traj, const KVGeoDetectorNode* node)
{
   current_nuc_recon = part;
 
   try {
      KVGroupReconstructor::ReconstructParticle(part, traj, node);
   }  catch (const std::exception& e) {
      KVError::Error(this, "ReconstructParticle", "Caught exception: %s", e.what());
      throw;
   }
 
   part->GetReconstructionTrajectory()->IterateFrom();
   while (auto Nd = part->GetReconstructionTrajectory()->GetNextNode()) {
      // energy_loss will contain the sum of energy losses (dE) measured in the *active* layer of each detector
      // the contribution to the energy of each particle may be greater than this after correcting for
      // dead layers etc.
      energy_loss[Nd->GetName()] = Nd->GetDetector()->GetEnergyLoss();
      ++number_uncalibrated[Nd->GetName()];
   }
}
 
 
 
 
void KVFilterGroupReconstructor::IdentifyParticle(KVReconstructedNucleus& PART)
{
   KVGroupReconstructor::IdentifyParticle(PART);
 
   if (PART.IsIdentified()) {
      PART.GetReconstructionTrajectory()->IterateFrom();
      KVGeoDetectorNode* node;
      while ((node = PART.GetReconstructionTrajectory()->GetNextNode())) {
         --number_unidentified[node->GetName()];
      }
   }
}
 
 
 
 
void KVFilterGroupReconstructor::CalibrateParticle(KVReconstructedNucleus* PART)
{
   // Calculate and set the energy of a (previously identified) reconstructed particle,
   // including an estimate of the energy loss in the target.
 
   PART->SetIsCalibrated();
   PART->SetECode(GetGroup()->GetParentStructure<KVMultiDetArray>()->GetNormalCalibrationCode()); // => general code for calibration with no problems, no calculated energy losses
   double Einc = 0;
   PART->GetReconstructionTrajectory()->IterateFrom();
   while (auto node = PART->GetReconstructionTrajectory()->GetNextNode()) {
      double edet, dE;
      auto det = node->GetDetector();
      //Info("Calib", "det=%s", det->GetName());
      //Info("Calib", "number_uncalibrated=%d, energy_loss=%f", number_uncalibrated[det->GetName()], energy_loss[det->GetName()]);
      // deal with pileup in detectors
      if (number_uncalibrated[det->GetName()] > 1 || (det != PART->GetStoppingDetector() && number_unidentified[det->GetName()] > 0)) {
         //Info("Calib", "number_uncalib[%s]=%d, Einc=%f", det->GetName(), number_uncalibrated[det->GetName()], Einc);
         // there are other uncalibrated particles which hit this detector
         // the contribution for this particle must be calculated from the residual energy (Einc)
         // deposited in all detectors so far (if any)
         if (Einc > 0) {
            // calculate dE in active layer for particle
            dE = det->GetDeltaEFromERes(PART->GetZ(), PART->GetA(), Einc);
            if (dE == 0) {
               //Info("Calib", "dE from Eres gives 0! dead");
               PART->SetECode(GetGroup()->GetParentStructure<KVMultiDetArray>()->GetNoCalibrationCode());
               PART->SetIsUncalibrated();
               break;
            }
            det->SetEResAfterDetector(Einc);
            edet = det->GetCorrectedEnergy(PART, dE);
            //Info("Calib", "calculated edet=%f from dE=%f", edet, dE);
            // negative parameter for calculated contribution
            PART->SetParameter(Form("%s.E%s", GetGroup()->GetParentStructure<KVMultiDetArray>()->GetName(), det->GetLabel()), -edet);
            PART->SetECode(GetGroup()->GetParentStructure<KVMultiDetArray>()->GetCalculatedCalibrationCode()); // calculated
         }
         else {
            //Info("Calib", "Einc=0, dead");
            // nothing to do here
            PART->SetECode(GetGroup()->GetParentStructure<KVMultiDetArray>()->GetNoCalibrationCode());
            PART->SetIsUncalibrated();
            break;
         }
      }
      else if (det->GetNHits() > 1) {
         // other particles hit this detector, but their contributions have already been calculated
         // and subtracted. use remaining energy calculated in detector.
         //Info("Calib", "det->nhits>1");
         dE = energy_loss[det->GetName()];
         det->SetEResAfterDetector(Einc);
         edet = det->GetCorrectedEnergy(PART, dE);
         // negative parameter for calculated contribution
         PART->SetParameter(Form("%s.E%s", GetGroup()->GetParentStructure<KVMultiDetArray>()->GetName(), det->GetLabel()), -edet);
         PART->SetECode(GetGroup()->GetParentStructure<KVMultiDetArray>()->GetCalculatedCalibrationCode()); // calculated
      }
      else {
         dE = det->GetEnergyLoss();
         det->SetEResAfterDetector(Einc);
         edet = det->GetCorrectedEnergy(PART, -1., (Einc > 0));
         //Info("Calib", "Simple normal dE=%f edet=%f", dE, edet);
         PART->SetParameter(Form("%s.E%s", GetGroup()->GetParentStructure<KVMultiDetArray>()->GetName(), node->GetDetector()->GetLabel()), edet);
      }
      Einc += edet;
      --number_uncalibrated[det->GetName()];
      energy_loss[det->GetName()] -= dE;
      det->SetEnergyLoss(energy_loss[det->GetName()]);
      //Info("Calib", "After: number_uncalibrated=%d, energy_loss=%f", number_uncalibrated[det->GetName()], energy_loss[det->GetName()]);
      //Info("Calib", "After: number_unidentified=%d", number_unidentified[det->GetName()]);
   }
   if (PART->GetECode() > 0) {
      PART->SetEnergy(Einc);
      // set particle momentum from telescope dimensions (random)
      PART->GetAnglesFromReconstructionTrajectory();
      //Info("CalibrateParticle","Calibrated particle with %f MeV",Einc);
      //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);
         //Info("CalibrateParticle","Target energy loss %f MeV",E_targ);
         PART->SetTargetEnergyLoss(E_targ);
      }
      Double_t E_tot = PART->GetEnergy() + E_targ;
      PART->SetEnergy(E_tot);
      //Info("CalibrateParticle","Total energy now %f MeV",E_tot);
 
      // to avoid inconsistencies in reconstructed data, now that we have 'calibrated' the particle,
      // we re-check that its identification is consistent with its apparent energy if the latter
      // is lower than the energy of the simulated particle
//      if (E_tot < PART->GetParameters()->GetDoubleValue("SIM:ENERGY")) {
//         //Warning("CalibrateParticle", "Apparent energy lower than real energy. Initial idendtification:");
//         auto idt = PART->GetIdentifyingTelescope();
//         auto IDR = PART->GetIdentificationResult(idt);
//         //std::cout << std::endl; IDR->Print(); std::cout << std::endl;
//         check_identification_threshold(PART, idt, IDR);
//         //Warning("CalibrateParticle", "Revised idendtification:");
//         //std::cout << std::endl; IDR->Print(); std::cout << std::endl;
//         if (IDR->IDOK) {
//            // Initial identification is still valid, but fine details may have changed...
//            PART->SetIdentification(IDR, idt);
//            // in case particle mass number changes, which changes the total (and kinetic) energy:
//            // we keep the same energy as calculated from measured losses in detectors,
//            // to stay coherent with the latter
//            PART->SetEnergy(E_tot);
//         }
//         else {
//            // initial identification is no longer valid. look for another?
//            // set first successful identification as particle identification
//            // as long as the id telescope concerned contains the stopping detector!
//            const KVSeqCollection* idt_list = PART->GetReconstructionTrajectory()->GetIDTelescopes();
//            TIter next(idt_list);
//            Int_t id_no = 1;
//            Bool_t ok = kFALSE;
//            KVIdentificationResult* pid = PART->GetIdentificationResult(id_no);
//            idt = (KVIDTelescope*)next();
//            while (idt && idt->HasDetector(PART->GetStoppingDetector())) {
//               if (pid->IDattempted && pid->IDOK) {
//                  ok = kTRUE;
//                  partID = *pid;
//                  identifying_telescope = idt;
//                  break;
//               }
//               ++id_no;
//               pid = PART->GetIdentificationResult(id_no);
//               idt = (KVIDTelescope*)next();
//            }
//            if (ok) { // particle was already identified
//               PART->SetIdentification(&partID, identifying_telescope);
//               // in case particle mass number changes, which changes the total (and kinetic) energy:
//               // we keep the same energy as calculated from measured losses in detectors,
//               // to stay coherent with the latter
//               PART->SetEnergy(E_tot);
//            }
//            else {
//               PART->SetIsUnidentified();
//               PART->SetIDCode(GetGroup()->GetParentStructure<KVMultiDetArray>()->GetBadIDCode());
//            }
//         }
//      }
 
   }
}
 
 
ClassImp(KVFilterGroupReconstructor)