KVDetectionSimulator.cpp

//Created by KVClassFactory on Sat Oct 10 09:37:42 2015
//Author: John Frankland,,,
 
#include "KVDetectionSimulator.h"
#include "KVGeoNavigator.h"
 
ClassImp(KVDetectionSimulator)
 
 
 
KVDetectionSimulator::KVDetectionSimulator(KVMultiDetArray* a, Double_t e_cut_off) :
   KVBase(Form("DetectionSimulator_%s", a->GetName()),
          Form("Simulate detection of particles or events in detector array %s", a->GetTitle())),
   fArray(a), fCalcTargELoss(kTRUE)
{
   // Initialise a detection simulator
   //
   // The detector array is put into simulation mode, and the minimum cut-off energy
   // for propagation of particles is set (1 keV is default)
   SetArray(a);
}
 
 
 
 
void KVDetectionSimulator::DetectEvent(KVEvent* event, const Char_t* detection_frame)
{
   //Simulate detection of event by multidetector array.
   //
   // optional argument detection_frame(="" by default) can be used to give name of
   // inertial reference frame (defined for all particles of 'event') to be used.
   // e.g. if the simulated event's default reference frame is the centre of mass frame, before calling this method
   // you should create the 'laboratory' or 'detector' frame with KVEvent::SetFrame(...), and then give the
   // name of the 'LAB' or 'DET' frame as 3rd argument here.
   //
   //For each particle in the event we calculate first its energy loss in the target (if the target has been defined, see KVMultiDetArray::SetTarget).
   //By default these energy losses are calculated from a point half-way along the beam-direction through the target (taking into account the orientation
   //of the target), if you want random depths for each event call GetTarget()->SetRandomized() before using DetectEvent().
   //
   //If the particle escapes the target then we look for the group in the array that it will hit. If there is one, then the detection of this particle by the
   //different members of the group is simulated.
   //
   //The detectors concerned have their fEloss members set to the energy lost by the particle when it crosses them.
   //
   //Give tags to the simulated particles via KVNucleus::AddGroup() method
   //Two general tags :
   //   - DETECTED : cross at least one active layer of one detector
   //   - UNDETECTED : go through dead zone or stopped in target
   //We add also different sub group :
   //   - For UNDETECTED particles : "NO HIT", "NEUTRON", "DEAD ZONE", "STOPPED IN TARGET" and "THRESHOLD", the last one concerns particles
   //     which go through the first detection stage of the multidetector array but stopped in an absorber (ie an inactive layer)
   //   - For DETECTED particles :
   //       *  "PUNCH THROUGH" corrresponds to particle which cross all the materials in front of it
   //         (high energy particle punch through)
   //       *  "INCOMPLETE"  corresponds to particles which stopped in the first detection stage of the multidetector
   //         in a detector which can not give alone a clear identification,
   //
 
   fDetectionFrame = detection_frame;
 
   if (get_array_navigator()->IsTracking()) {
      // clear any tracks created by last event
      gGeoManager->ClearTracks();
      get_array_navigator()->ResetTrackID();
   }
 
   // Reset detectors in array hit by any previous events
   ClearHitGroups();
 
   // if any particle traverses >1 group of the array, we need to add copies of the particle
   // to the simulated event, 1 for each group crossed (with the corresponding informations
   // for reconstruction procedures)
   KVUnownedList multi_group_parts;
 
   // iterate through the particles of the event
   for (auto& part : EventIterator(event)) {
      // reference to particle in requested detection frame
      auto part_to_detect = (KVNucleus*)part.GetFrame(detection_frame, kFALSE);
 
      // store initial energy of particle in detection frame
      part_to_detect->SetE0();
      part.SetParameter("SIM:Z", part.GetZ());
      part.SetParameter("SIM:A", part.GetA());
      part.SetParameter("SIM:ENERGY", part_to_detect->GetE());
      part.SetParameter("SIM:THETA", part_to_detect->GetTheta());
      part.SetParameter("SIM:PHI", part_to_detect->GetPhi());
 
      // neutral particles & those with less than the cut-off energy are not detected
      if ((part.GetZ() == 0) && !get_array_navigator()->IsTracking()) {
         // neutrons are included in tracking, if active
         part.SetParameter("UNDETECTED", "NEUTRON");
      }
      else if (part.GetZ() && !get_array_navigator()->CheckIonForRangeTable(part.GetZ(), part.GetA())) {
         // ignore charged particles which range table cannot handle
         part.SetParameter("UNDETECTED", "NOT IN RANGE TABLE");
      }
      else if (!fGeoFilter && (part_to_detect->GetEnergy() < GetMinKECutOff())) {
         part.SetParameter("UNDETECTED", "AT REST");
      }
      else {
         if (IncludeTargetEnergyLoss() && GetTarget() && part.GetZ()) {
            //simulate passage through target material
            auto ebef = part_to_detect->GetE();
            GetTarget()->DetectParticle(part_to_detect);
            auto eLostInTarget = ebef - part_to_detect->GetE();
            part.SetParameter("ENERGY LOSS IN TARGET", eLostInTarget);
            if (part_to_detect->GetE() < GetMinKECutOff())
               part.SetParameter("UNDETECTED", "STOPPED IN TARGET");
         }
 
         if (fGeoFilter || (part_to_detect->GetE() > GetMinKECutOff())) {
 
            auto nvl = PropagateParticle(part_to_detect);
 
            if (nvl.IsEmpty()) {
               if (part.GetZ() == 0) {
                  // tracking
                  part.SetParameter("UNDETECTED", "NEUTRON");
               }
               else {
                  if (part.GetParameters()->HasParameter("DEADZONE")) {
                     // deadzone
                     part.SetParameter("UNDETECTED", "DEAD ZONE");
                  }
                  else {
                     // missed all detectors
                     part.SetParameter("UNDETECTED", "NO HIT");
                  }
               }
            }
            else {
               // check for incomplete stopping of particle => PUNCH THROUGH
               // note that energy losses are not calculated in DEADZONE volume; as soon as a particle enters
               // a DEADZONE volume its propagation stops. therefore particles which lose energy in one or more
               // active material volumes and then hit a DEADZONE may have residual kinetic energy also
               if (part_to_detect->GetE() > GetMinKECutOff()) {
                  part.SetParameter("RESIDUAL ENERGY", part_to_detect->GetE());
                  part.SetParameter("DETECTED", "PUNCH THROUGH");
               }
               else
                  part.SetParameter("DETECTED", "OK");
            }
 
            if (!nvl.IsEmpty()) {
               for (Int_t ii = 0; ii < nvl.GetNpar(); ++ii) {
                  part.SetParameter(nvl.GetNameAt(ii), nvl.GetDoubleValue(ii));
               }
            }
 
            if(part.GetParameters()->HasParameter("MULTIGROUP"))
               multi_group_parts.Add(&part);
 
         }
      }
 
      part_to_detect->SetMomentum(*part_to_detect->GetPInitial());
   }
 
   if(!multi_group_parts.IsEmpty())
   {
      // for each particle with parameter MULTIGROUP we make copies of the particle
      // (1 for each extra group) and add them to the event
      for(auto _p : multi_group_parts)
      {
         auto part = dynamic_cast<KVNucleus*>(_p);
         part->GetParameters()->SetValue("DETECTED","OK");
         auto multigroup = part->GetParameters()->GetIntValue("MULTIGROUP");
         for(int mgrp = 1; mgrp<multigroup; ++mgrp)
         {
            KVNucleus* nuc;
            part->Copy(*(nuc = event->AddNucleus()));
            // now remove/rename parameters in each particle:
            //  - in nuc, GROUP_mgrp -> GROUP, TRAJECTORY_mgrp -> TRAJECTORY, etc.
            nuc->GetParameters()->SetValue("GROUP", nuc->GetParameters()->GetIntValue(Form("GROUP_%d",mgrp)));
            nuc->GetParameters()->SetValue("TRAJECTORY", nuc->GetParameters()->GetStringValue(Form("TRAJECTORY_%d",mgrp)));
            nuc->GetParameters()->SetValue("STOPPING DETECTOR", nuc->GetParameters()->GetStringValue(Form("STOPPING DETECTOR_%d",mgrp)));
            nuc->GetParameters()->RemoveParameter(Form("GROUP_%d",mgrp));
            nuc->GetParameters()->RemoveParameter(Form("TRAJECTORY_%d",mgrp));
            nuc->GetParameters()->RemoveParameter(Form("STOPPING DETECTOR_%d",mgrp));
            part->GetParameters()->RemoveParameter(Form("GROUP_%d",mgrp));
            part->GetParameters()->RemoveParameter(Form("TRAJECTORY_%d",mgrp));
            part->GetParameters()->RemoveParameter(Form("STOPPING DETECTOR_%d",mgrp));
            // parse trajectory into list of detectors
            KVString traj(nuc->GetParameters()->GetStringValue("TRAJECTORY"));
            traj.Begin("/");
            while(!traj.End())
            {
               auto ddet = traj.Next();
               // any parameters in 'part's list which contain the name of this detector are removed
               KVUnownedList to_remove;
               TIter it_par(part->GetParameters()->GetList());
               KVNamedParameter* np;
               while( (np = (KVNamedParameter*)it_par()) )
               {
                  if(TString(np->GetName()).Contains(ddet))
                     to_remove.Add(np);
               }
               if(!to_remove.IsEmpty())
               {
                  for(auto rp : to_remove)
                  {
                     part->GetParameters()->RemoveParameter(rp->GetName());
                  }
               }
            }
            // now do the same for 'nuc's list of parameters
            traj = part->GetParameters()->GetStringValue("TRAJECTORY");
            traj.Begin("/");
            while(!traj.End())
            {
               auto ddet = traj.Next();
               // any parameters in 'nuc's list which contain the name of this detector are removed
               KVUnownedList to_remove;
               TIter it_par(nuc->GetParameters()->GetList());
               KVNamedParameter* np;
               while( (np = (KVNamedParameter*)it_par()) )
               {
                  if(TString(np->GetName()).Contains(ddet))
                     to_remove.Add(np);
               }
               if(!to_remove.IsEmpty())
               {
                  for(auto rp : to_remove)
                  {
                     nuc->GetParameters()->RemoveParameter(rp->GetName());
                  }
               }
            }
         }
      }
   }
 
}
 
 
 
 
 
KVNameValueList KVDetectionSimulator::PropagateParticle(KVNucleus* part)
{
   // Simulate detection of a single particle
   //
   // Propagate particle through the array,
   // calculating its energy losses in all absorbers, and setting the
   // energy loss members of the active detectors on the way.
   //
   // Returns a list containing the name and energy loss of each
   // detector hit in array (list is empty if none i.e. particle
   // in beam pipe or dead zone of the multidetector)
 
   auto nparams = part->GetParameters()->GetNpar();
 
   get_array_navigator()->PropagateParticle(part);
 
   // particle missed all detectors
   if (part->GetParameters()->GetNpar() == nparams ||
         ((part->GetParameters()->GetNpar() - nparams) == 1 && part->GetParameters()->HasParameter("DEADZONE")))
      return KVNameValueList();
 
   // list of energy losses in active layers of detectors
   KVNameValueList NVL;
 
   // trajectories followed in each group crossed by particle
   // (taking into account possibility that particle may traverse more than 1 group)
   std::unordered_map<unsigned int,TString> traj_group_map;
   // find detectors in array hit by particle
   KVDetector* last_detector = nullptr;
   TIter next(part->GetParameters()->GetList());
   KVNamedParameter* param;
   while ((param = (KVNamedParameter*)next())) {
      KVString pname(param->GetName());
      pname.Begin(":");
      KVString pn2 = pname.Next();
      KVString pn3 = pname.Next();
      if (pn2 == "DE") {
         pn3.Begin("/");
         KVString det_name = pn3.Next();
         if (pn3.End() || pn3.Next().BeginsWith("ACTIVE")) {
            // energy loss in active layer of detector
            KVDetector* curDet = fArray->GetDetector(det_name);
            if (curDet) {
               last_detector = curDet;
               Double_t de = param->GetDouble();
               NVL.SetValue(curDet->GetName(), de);
               // add detector to name of trajectory followed by particle
               TString traj;
               if (part->GetParameters()->HasStringParameter("TRAJECTORY")) {
                  traj = part->GetParameters()->GetStringValue("TRAJECTORY");
                  traj.Prepend(Form("%s/", det_name.Data()));
               }
               else {
                  traj = Form("%s/", det_name.Data());
               }
               part->SetParameter("TRAJECTORY", traj);
               // set number of group where detected
               part->SetParameter("GROUP", (int)curDet->GetGroupNumber());
               auto group_num = curDet->GetGroupNumber();
               auto& traj_string = traj_group_map[group_num];
               if(traj_string.IsNull()) // begin new trajectory?
                  traj_string.Form("%s/", det_name.Data());
               else
                  traj_string.Prepend(Form("%s/", det_name.Data()));
            }
         }
      }
   }
   // analyse resulting map of groups & trajectories
   if(traj_group_map.size() > 1)
   {
      // multiple groups hit: need to add copy of particle for each group
      int itraj{0};
      for(auto& tgp : traj_group_map)
      {
         fHitGroups.AddGroup(fArray->GetGroup(tgp.first));
         KVString traj(tgp.second);
         traj.Begin("/"); // first detector in trajectory is stopping detector (for this trajectory)
         if(itraj)
         {
            part->SetParameter(Form("TRAJECTORY_%d",itraj), tgp.second);
            part->SetParameter(Form("GROUP_%d",itraj), (int)tgp.first);
            part->SetParameter(Form("STOPPING DETECTOR_%d",itraj), traj.Next());
         }
         else
         {
            part->SetParameter("TRAJECTORY", tgp.second);
            part->SetParameter("GROUP", (int)tgp.first);
            part->SetParameter("STOPPING DETECTOR", traj.Next());
         }
         ++itraj;
      }
      if(itraj>1) part->SetParameter("MULTIGROUP",itraj);
   }
   else
   {
      // add hit group to list if not already in it
      if (last_detector) {
         fHitGroups.AddGroup(last_detector->GetGroup());
         part->SetParameter("STOPPING DETECTOR", last_detector->GetName());
      }
   }
   return NVL;
}