KVReconstructedNucleus.cpp

#include "KVReconstructedNucleus.h"
#include "KVIDTelescope.h"
#include "KVGroup.h"
#include "KVMultiDetArray.h"
 
using namespace std;
 
ClassImp(KVReconstructedNucleus);
 
 
 
void KVReconstructedNucleus::init()
{
   //default initialisation
   fReconTraj = nullptr;
   fRealZ = fRealA = 0.;
   fDetNames = "/";
   fIDTelName = "";
   fIDTelescope = nullptr;
   fNSegDet = 0;
   fAnalStatus = 99;
   fTargetEnergyLoss = 0;
   ResetBit(kIsIdentified);
   ResetBit(kIsCalibrated);
   ResetBit(kCoherency);
   ResetBit(kZMeasured);
   ResetBit(kAMeasured);
}
 
 
 
 
 
KVReconstructedNucleus::KVReconstructedNucleus() : fIDResults("KVIdentificationResult", 5)
{
   //default ctor.
   init();
}
 
 
 
 
KVReconstructedNucleus::KVReconstructedNucleus(const KVReconstructedNucleus& obj)
   : KVNucleus(), fIDResults("KVIdentificationResult", 5)
{
   //copy ctor
   init();
#if ROOT_VERSION_CODE >= ROOT_VERSION(3,4,0)
   obj.Copy(*this);
#else
   ((KVReconstructedNucleus&) obj).Copy(*this);
#endif
}
 
 
 
 
 
void KVReconstructedNucleus::Streamer(TBuffer& R__b)
{
   // Stream an object of class KVReconstructedNucleus.
   //Customized streamer.
 
   UInt_t R__s, R__c;
   if (R__b.IsReading()) {
      Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
      if (R__v < 17) {
         // Before v17, fIDResults was a static array: KVIdentificationResult fIDresults[5]
         // We convert this to the new TClonesArray format
         KVNucleus::Streamer(R__b);
         fDetNames.Streamer(R__b);
         fIDTelName.Streamer(R__b);
         if (R__v >= 16) {
            R__b >> fNSegDet;
            R__b >> fAnalStatus;
         }
         R__b >> fRealZ;
         R__b >> fRealA;
         R__b >> fTargetEnergyLoss;
         int R__i;
         KVIdentificationResult id_array[5];
         for (R__i = 0; R__i < 5; R__i++) {
            id_array[R__i].Streamer(R__b);
            id_array[R__i].Copy(*GetIdentificationResult(R__i + 1));
         }
         R__b.CheckByteCount(R__s, R__c, KVReconstructedNucleus::IsA());
      }
      else
         R__b.ReadClassBuffer(KVReconstructedNucleus::Class(), this, R__v, R__s, R__c);
      // if the multidetector object exists, update some informations
      // concerning the detectors etc. hit by this particle
      if (gMultiDetArray) {
         fReconTraj = nullptr;
         RebuildReconTraj();
         fIDTelescope = nullptr;
         if (fIDTelName != "") fIDTelescope = gMultiDetArray->GetIDTelescope(fIDTelName.Data());
         if (fReconTraj) {
            if (R__v < 16) fNSegDet = fReconTraj->GetNumberOfIndependentIdentifications(); // fNSegDet/fAnalStatus non-persistent before v.16
         }
         else {
            TIter next_det(&fDetList);
            KVDetector* det;
            while ((det = (KVDetector*)next_det())) {
               det->AddHit(this);
               if (det->IsDetecting()) { //to be coherent with AddDetector() method
                  //modify detector's counters depending on particle's identification state
                  if (IsIdentified())
                     det->IncrementIdentifiedParticles();
                  else
                     det->IncrementUnidentifiedParticles();
               }
            }
         }
      }
   }
   else {
      R__b.WriteClassBuffer(KVReconstructedNucleus::Class(), this);
   }
}
 
 
 
 
 
void KVReconstructedNucleus::PrintStatusString() const
{
   switch (GetStatus()) {
      case kStatusOK:
         cout <<
              "Particle alone in group, or identification independently of other particles in group is directly possible." << endl;
         break;
 
      case kStatusOKafterSub:
         cout <<
              "Particle reconstructed after identification of others in group and subtraction of their calculated energy losses in common detectors."
              << endl;
         break;
 
      case kStatusOKafterShare:
         cout <<
              "Particle identification estimated after arbitrary sharing of energy lost in common detectors between several reconstructed particles."
              << endl;
         break;
 
      case kStatusStopFirstStage:
         cout <<
              "Particle stopped in first stage of telescope. Estimation of minimum Z."
              << endl;
         break;
 
      case kStatusPileupDE:
         cout <<
              "Undetectable pile-up in first member of identifying telesscope (apparent status=OK). Would lead to incorrect identification by DE-E method (Z and/or A overestimated)."
              << endl;
         break;
 
      case kStatusPileupGhost:
         cout <<
              "Undetectable ghost particle in filtered simulation. Another particle passed through all of the same detectors (pile-up)."
              << endl;
         break;
 
 
      default:
         cout << GetStatus() << endl;
         break;
   }
}
 
 
 
 
Bool_t KVReconstructedNucleus::InArray(const TString& name) const
{
   // Returns kTRUE if particle was detected in array with given name
   return GetArrayName() == name;
}
 
 
 
 
TString KVReconstructedNucleus::GetArrayName() const
{
   // Returns name of array particle was detected in (if known)
   if (GetParameters()->HasStringParameter("ARRAY")) return GetParameters()->GetStringValue("ARRAY");
   if (GetStoppingDetector()) return GetStoppingDetector()->GetGroup()->GetParentStructure<KVMultiDetArray>()->GetName();
   return "";
}
 
 
 
 
void KVReconstructedNucleus::Print(Option_t*) const
{
 
   int ndets = GetNumDet();
   if (ndets) {
 
      for (int i = ndets - 1; i >= 0; i--) {
         KVDetector* det = GetDetector(i);
         if (det) det->Print("data");
      }
      for (int i = 1; i <= GetNumberOfIdentificationResults(); i++) {
         KVIdentificationResult* idr = const_cast<KVReconstructedNucleus*>(this)->GetIdentificationResult(i);
         if (idr->IDattempted) idr->Print();
      }
   }
   if (GetStoppingDetector()) cout << "STOPPED IN : " <<
                                      GetStoppingDetector()->GetName() << endl;
   if (IsIdentified()) {
      if (GetIdentifyingTelescope()) cout << "IDENTIFIED IN : " <<
                                             GetIdentifyingTelescope()->GetName() << endl;
      cout << " =======> ";
      if(GetZ())
      {
         cout << " Z=" << GetZ() << " A=" << GetA();
         if (IsAMeasured()) cout << " Areal=" << GetRealA();
         else cout << " Zreal=" << GetRealZ();
      }
      else if(GetA()==1)
      {
         cout << "neutron";
      }
      else
      {
         cout << "gamma";
      }
   }
   else {
      cout << "(unidentified)" << endl;
   }
   if (IsCalibrated()) {
      cout << " Total Energy = " << GetEnergy() << " MeV,  Theta=" << GetTheta() << " Phi=" << GetPhi() << endl;
      cout << "    Target energy loss correction :  " << GetTargetEnergyLoss() << " MeV" << endl;
   }
   else {
      cout << " (uncalibrated)" << endl;
   }
   cout << "RECONSTRUCTION STATUS : " << endl;
   PrintStatusString();
   if (fReconTraj) fReconTraj->ls();
   if (GetParameters()->GetNpar()) GetParameters()->Print();
}
 
 
 
#if ROOT_VERSION_CODE >= ROOT_VERSION(3,4,0)
 
 
void KVReconstructedNucleus::Copy(TObject& obj) const
#else
void KVReconstructedNucleus::Copy(TObject& obj)
#endif
{
   //
   //Copy this to obj
   //
   KVNucleus::Copy(obj);
   KVReconstructedNucleus& robj = (KVReconstructedNucleus&)obj;
   robj.SetIdentifyingTelescope(GetIdentifyingTelescope());
   robj.fDetNames = fDetNames;
   robj.SetRealZ(GetRealZ());
   robj.SetRealA(GetRealA());
   robj.SetTargetEnergyLoss(GetTargetEnergyLoss());
   robj.fReconTraj = fReconTraj;
   robj.fIDTelescope = fIDTelescope;
   robj.fNSegDet = fNSegDet;
   robj.fAnalStatus = fAnalStatus;
   // copy id results
   Int_t nidres = GetNumberOfIdentificationResults();
   for (int i = 1; i <= nidres; ++i) {
      GetIdentificationResult(i)->Copy(*robj.GetIdentificationResult(i));
   }
   robj.SetBit(kIsIdentified, TestBit(kIsIdentified));
   robj.SetBit(kIsCalibrated, TestBit(kIsCalibrated));
   robj.SetBit(kCoherency, TestBit(kCoherency));
   robj.SetBit(kZMeasured, TestBit(kZMeasured));
   robj.SetBit(kAMeasured, TestBit(kAMeasured));
}
 
 
 
 
void KVReconstructedNucleus::CopyAndMoveReferences(const KVReconstructedNucleus* other)
{
   // Copy all characteristics of 'other' and also change all references to
   // 'other' to references to 'this' (i.e. in detectors hit by particle).
   // 'other' will not be fully valid after this operation (shouldn't be used further)
 
   other->Copy(*this);
   KVGeoDetectorNode* node;
   const KVReconNucTrajectory* traj = other->GetReconstructionTrajectory();
   traj->IterateFrom();
   while ((node = traj->GetNextNode())) {
      KVDetector* d = node->GetDetector();
      // if 'other' was in the detector's hitlist, we add 'this' in its place
      if(d->RemoveHit(const_cast<KVReconstructedNucleus*>(other)))
         d->AddHit(this);
   }
}
 
 
 
 
void KVReconstructedNucleus::Clear(Option_t* opt)
{
   // Reset nucleus. Calls KVNucleus::Clear.
   // if opt!="N": Calls KVGroup::Reset for the group where it was reconstructed.
 
   KVNucleus::Clear(opt);
   if (GetGroup() && strncmp(opt, "N", 1))
      GetGroup()->Reset();
   fDetList.Clear();
   fIDResults.Clear("C");
   init();
}
 
 
 
 
void KVReconstructedNucleus::ModifyReconstructionTrajectory(const KVReconNucTrajectory* t)
{
   // this method modifies the reconstructed trajectory
   // probably called during the identification coherency check
   fReconTraj = t;
   fDetNames = t->GetPathString();//to store/retrieve reconstruction trajectory on file
}
 
 
 
 
void KVReconstructedNucleus::SetReconstructionTrajectory(const KVReconNucTrajectory* t)
{
   // Method called in initial reconstruction of particle.
 
   fReconTraj = t;
   fNSegDet = t->GetNumberOfIndependentIdentifications();
   t->AddUnidentifiedParticle(0);//add 1 unidentified particle to each detector on trajectory
   fDetNames = t->GetPathString();//to store/retrieve reconstruction trajectory on file
}
 
 
 
 
void KVReconstructedNucleus::GetAnglesFromReconstructionTrajectory(Option_t* opt)
{
   // Calculate angles theta and phi for nucleus based on the detectors on its reconstruction trajectory.
   // The momentum is set using these angles, its mass and its kinetic energy.
   //
   // The detector with the smallest solid angle along the trajectory is the one which defines the
   // angles for the reconstructed particle.
   //
   // The (optional) option string can be "random" or "mean":
   //
   // If "random" (default) the angles are drawn at random between the over the surface of the detector.
   //
   // If "mean" the (theta,phi) position of the centre of the detector is used to fix the nucleus' direction.
 
 
   if (GetEnergy() <= 0.0)//don't try if particle has no correctly defined energy
      return;
   if (!GetReconstructionTrajectory())
      return;
 
   KVDetector* angle_det = nullptr;
   double small_solid = 1.e+09;
   GetReconstructionTrajectory()->IterateFrom();
   KVGeoDetectorNode* n;
   while ((n = GetReconstructionTrajectory()->GetNextNode())) {
      if (n->GetDetector()->GetSolidAngle() < small_solid) {
         angle_det = n->GetDetector();
         small_solid = n->GetDetector()->GetSolidAngle();
      }
   }
   if (!strcmp(opt, "random")) {
      //random angles
      TVector3 dir = angle_det->GetRandomDirection("random");
      SetMomentum(GetEnergy(), dir);
   }
   else {
      //middle of telescope
      TVector3 dir = angle_det->GetDirection();
      SetMomentum(GetEnergy(), dir);
   }
}
 
 
 
 
void KVReconstructedNucleus::RebuildReconTraj()
{
   // Called by Streamer when reading in data
   // The fDetNames string is used to associate the particle with its reconstruction trajectory
 
   if (gMultiDetArray) {
      // for old data, detnames was written as "/DET_1/DET_2/..."
      // trajectory paths are written as "DET_1/DET_2/..."
      if (fDetNames[0] == '/') fDetNames.Remove(0, 1);
      fDetNames.Begin("/");
      KVString n_stop_det = fDetNames.Next();
      KVDetector* stop_det = gMultiDetArray->GetDetector(n_stop_det);
      if (stop_det) {
         KVGroup* gr = stop_det->GetGroup();
         gr->AddHit(this);//for backwards compatibility, group coherency analysis performed in KVReconstructedEvent::Streamer (old data)
         fReconTraj = (KVReconNucTrajectory*)gr->FindReconTraj(fDetNames);
         if (fReconTraj) {
            fReconTraj->IterateFrom();
            KVGeoDetectorNode* n;
            while ((n = fReconTraj->GetNextNode())) {
               n->GetDetector()->AddHit(this);
               if (IsIdentified()) n->GetDetector()->IncrementIdentifiedParticles();
               else n->GetDetector()->IncrementUnidentifiedParticles();
            }
         }
      }
   }
}
 
 
 
 
void KVReconstructedNucleus::ReplaceReconTraj(const TString& traj_name)
{
   // Change the particle's reconstruction trajectory to a different one starting from the
   // same stopping detector (and therefore in the same group).
   //
   // trajectory paths are written as "DET_1/DET_2/..."
 
   KVReconNucTrajectory* new_traj = (KVReconNucTrajectory*)GetGroup()->FindReconTraj(traj_name);
   if (!new_traj) {
      Error("ReplaceReconTraj", "Trajectory %s not found - meant to replace %s", traj_name.Data(), fReconTraj->GetTitle());
      return;
   }
   // first iterate over existing trajectory and reset detectors
   fReconTraj->IterateFrom();
   KVGeoDetectorNode* n;
   while ((n = fReconTraj->GetNextNode())) {
      n->GetDetector()->RemoveHit(this);
      if (IsIdentified()) n->GetDetector()->IncrementIdentifiedParticles(-1);
      else n->GetDetector()->IncrementUnidentifiedParticles(-1);
   }
   // set new trajectory
   fReconTraj = new_traj;
   fReconTraj->IterateFrom();
   while ((n = fReconTraj->GetNextNode())) {
      n->GetDetector()->AddHit(this);
      if (IsIdentified()) n->GetDetector()->IncrementIdentifiedParticles();
      else n->GetDetector()->IncrementUnidentifiedParticles();
   }
}
 
 
 
 
void KVReconstructedNucleus::SetIdentification(KVIdentificationResult* idr, KVIDTelescope* idt)
{
   // Set identification of nucleus from informations in identification result object
   // The mass (A) information in KVIdentificationResult is only used if the mass
   // was measured as part of the identification. Otherwise the nucleus' mass formula
   // will be used to calculate A from the measured Z.
   //
   // The identifying telescope is set to idt.
 
   SetIdentifyingTelescope(idt);
   SetIDCode(idr->IDcode);
   SetZMeasured(idr->Zident);
   SetAMeasured(idr->Aident);
   SetZ(idr->Z);
   if (idr->Aident) {
      SetA(idr->A);
      SetRealA(idr->PID);
   }
   else {
      SetRealZ(idr->PID);
   }
}
 
 
 
 
void KVReconstructedNucleus::ls(Option_t*) const
{
   printf(" A:%6s", GetParameters()->GetStringValue("ARRAY"));
   if (GetStoppingDetector()) printf("  D:%10s", GetStoppingDetector()->GetName());
   printf(" IDCODE=%2d", GetIDCode());
   if (IsIdentified()) {
      if (GetIdentifyingTelescope()) printf(" ID:%15s", GetIdentifyingTelescope()->GetName());
      if (IsZMeasured()) printf(" Z=%2d", GetZ());
      else printf("     ");
      if (IsAMeasured()) printf(" A=%3d  : ", GetA());
      else printf("        : ");
      if (IsCalibrated()) printf(" E=%g MeV", GetEnergy());
      if (GetParameters()->IsValue("Coherent", false)) printf("/not coherent/");
      if (GetParameters()->IsValue("Pileup", true)) printf("/pileup/");
   }
   printf("\n");
}