TTruthTrajectoriesModule.cxx

Go to the documentation of this file.

#include <TTruthTrajectoriesModule.hxx>

// c++
#include <cstdlib>
#include <set>
// Root
#include <TLorentzVector.h>
// oaEvent
#include <HEPUnits.hxx>
#include <TG4Trajectory.hxx>
#include <TND280Log.hxx>
// oaGeomInfo
#include <TGeomInfo.hxx>

ClassImp(ND::TTruthTrajectoriesModule::TTruthTrajectory);

#include "../cvstags/TTruthTrajectoriesModule.cxx"


ND::TTruthTrajectoriesModule::TTruthTrajectoriesModule(const char* name,
                                                       const char* title)
    : fMaxNTrajectories(500),
      fMinLength(0),
      fSaveOnlyP0DTrackerECALTrajectories(kTRUE),
      fSaveParentChain(kTRUE),
      fNTrajLepton(0),
      fNTrajBaryon(0),
      fNTrajMeson(0),
      fNTrajPhoton(0),
      fNTrajOtherCharged(0),
      fNTrajOtherNeutral(0),
      fTrajectories(
          new TClonesArray("ND::TTruthTrajectoriesModule::TTruthTrajectory",
                           fMaxNTrajectories)) {
  SetNameTitle(name, title);
  fIsEnabled = kTRUE;
  fDescription =
      "Module storing the truth information for particle trajectories";
  fCVSTagName = CVSTAG;
  fCVSID = CVSID;
}

ND::TTruthTrajectoriesModule::~TTruthTrajectoriesModule() {
  delete fTrajectories;
}

Bool_t ND::TTruthTrajectoriesModule::Configure(std::string& option) {
  if (option == "saveall" || option == "SaveAll") {
    SetSaveOnlyP0DTrackerECALTrajectories(kFALSE);
    ND280Info("Setting TTruthTrajectoriesModule to Save All Trajectories");
    return kTRUE;
  }
  return kFALSE;
}

void ND::TTruthTrajectoriesModule::InitializeBranches() {
  fOutputTree->Branch("NTraj", &fNTraj, "NTraj/I", fBufferSize)->SetTitle(" Total number of trajectories saved from the event");
  fOutputTree->Branch("NTrajLepton", &fNTrajLepton, "NTrajLepton/I",
                      fBufferSize)->SetTitle("Number of Charged Lepton trajectories saved from the event ");
  fOutputTree->Branch("NTrajBaryon", &fNTrajBaryon, "NTrajBaryon/I",
                      fBufferSize)->SetTitle(" Number of Charged Baryon Trajectories saved from the event ");
  fOutputTree->Branch("NTrajMeson", &fNTrajMeson, "NTrajMeson/I", fBufferSize)->SetTitle("Number of Charged Meson Trajectories saved from the event ");
  fOutputTree->Branch("NTrajPhoton", &fNTrajPhoton, "NTrajPhoton/I",
                      fBufferSize)->SetTitle(" Number of Photon Trajectories saved from the event ");
  fOutputTree->Branch("NTrajOtherCharged", &fNTrajOtherCharged,
                      "NTrajOtherCharged/I", fBufferSize)->SetTitle("Number of Other Charged Trajectories saved from the event ");
  fOutputTree->Branch("NTrajOtherNeutral", &fNTrajOtherNeutral,
                      "NTrajOtherNeutral/I", fBufferSize)->SetTitle(" Number of Other Neutral Trajectories saved from the event");
  fOutputTree->Branch("NTrajOther", &fNTrajOther, "NTrajOther/I", fBufferSize)->SetTitle(" Number of Any Other Trajectories saved from the event");

  fOutputTree->Branch("Trajectories", "TClonesArray", &fTrajectories,
                      fBufferSize, fSplitLevel)->SetTitle(" Clones array of ND::TTruthTrajectoriesModule::TTruthTrajectory  sorted in ascending ID order");
}

Bool_t ND::TTruthTrajectoriesModule::ProcessFirstEvent(ND::TND280Event& event) {
  bool IsDataFile = event.GetContext().IsDetector();
  if (IsDataFile) {
    throw ND::EDataFile();
  }
  return true;
}

void ND::TTruthTrajectoriesModule::FillSaveList(
    ND::THandle<ND::TG4TrajectoryContainer> trajectories) {
  // Create a set of trajectory IDs for trajectories which satisfy the criteria
  // to be saved

  ND::TG4TrajectoryContainer::iterator traj_it = trajectories->begin();
  ND::TG4TrajectoryContainer::iterator traj_end = trajectories->end();
  ND::TG4Trajectory* traj;

  for (; traj_it != traj_end; ++traj_it) {
    traj = &traj_it->second;

    if (!this->SaveTraj(traj)) continue;

    fSaveList.insert(traj->GetTrackId());
  }

  if (fSaveParentChain) {
    // Now need to save the parent chain of every trajectory which was
    // identified in the first set.

    // Create a copy of the original save set, else our iterators will become
    // invalid as we add more elements.
    std::set<int> original_savelist = fSaveList;

    // For each element in the original save set, iterate up its parent chain
    // and make sure the entire chain is in the complete save set. This
    // procedure is not optimised, but it is simple.
    int id;
    std::set<int>::iterator original_savelist_it = original_savelist.begin();
    std::set<int>::iterator original_savelist_end = original_savelist.end();
    for (; original_savelist_it != original_savelist_end;
         ++original_savelist_it) {
      id = *original_savelist_it;

      while (id > 0) {
        fSaveList.insert(id);
        id = trajectories->GetTrajectory(id)->GetParentId();
      }
    }
  }
}

bool ND::TTruthTrajectoriesModule::FillTree(ND::TND280Event& event) {
  if (!event.GetContext().IsMC()) {
    ND280Info(
        "Event Context reports event is non-MC. Assuming entire file is non-MC "
        "and disabling the module.");
    return false;  // disable module
  }

  fNTraj = 0;
  fNTrajLepton = 0;
  fNTrajBaryon = 0;
  fNTrajMeson = 0;
  fNTrajPhoton = 0;
  fNTrajOtherCharged = 0;
  fNTrajOtherNeutral = 0;
  fNTrajOther = 0;
  fTrajectories->Clear();

  if (!event.Has<ND::TG4TrajectoryContainer>("truth/G4Trajectories")) {
    ND280Verbose("No trajectories container found, skipping event.");
    return true;  // this happens occasionally, so don't disable the module
  }

  ND::THandle<ND::TG4TrajectoryContainer> trajectories =
      event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");

  fSaveList.clear();
  FillSaveList(trajectories);

  ND::TG4TrajectoryContainer::iterator traj_it = trajectories->begin();
  ND::TG4TrajectoryContainer::iterator traj_end = trajectories->end();
  ND::TG4Trajectory* traj;
  ND::TTruthTrajectoriesModule::TTruthTrajectory* trajToFill;

  for (; traj_it != traj_end; ++traj_it) {
    traj = &traj_it->second;

    // Check if this trajectory is worth saving
    if (fSaveList.count(traj->GetTrackId()) == 0) continue;

    // Protect from empty trajectories (bug encountered while testing)
    if (traj->GetTrackId() == 0) {
      ND280Error("Found TG4Trajectory with ID = 0");
      continue;
    }

    //
    // Do the saving

    trajToFill = new ((*fTrajectories)[fNTraj]) TTruthTrajectory();
    ++fNTraj;

    trajToFill->ID = traj->GetTrackId();
    trajToFill->PDG = traj->GetPDGEncoding();
    trajToFill->Name = traj->GetParticleName();
    trajToFill->Category = GetCategory(traj);
    trajToFill->ParentID = traj->GetParentId();
    trajToFill->PrimaryID = trajectories->GetPrimaryId(trajToFill->ID);

    this->FillTrajectoryPoints(traj, trajToFill);
    // Trace variables
    this->FillTraces(traj, trajToFill);

    trajToFill->InitMomentum = traj->GetInitialMomentum();
    trajToFill->InitPosition = traj->GetInitialPosition();
    trajToFill->FinalPosition = traj->GetFinalPosition();

    // Charge and Mass

    // ND::G4Trajectory::GetParticle can return 0 if the particle's
    // PDG number is not found in the PDG database. eventAnalysis defines
    // its own extended database but should still check in case a
    // new nucleus is added to the MC, for example.
    // At the moment this is also triggered by isomers - need to handle this
    if (traj->GetParticle() == 0) {
      int pdg = traj->GetPDGEncoding();
      int i = 0;
      int a = pdg / 10 % 1000;
      int z = pdg / 10000 % 1000;
      TGeoElementRN* rn = TGeoElement::GetElementTable()->GetElementRN(a, z, i);
      trajToFill->Mass = rn->MassNo() * 938.272046;  // ROOT TGeoElementRN
                                                     // stores mass as nucleon
                                                     // number therefore convert
                                                     // to MeV/c2
      trajToFill->Charge = static_cast<Int_t>(rn->AtomicNo() * 3);
    } else if (traj->GetParticle()) {
      // In ROOT the TParticlePDG mass is given in GeV, but we
      // want MeV
      trajToFill->Mass = traj->GetParticle()->Mass() * 1000.0;
      trajToFill->Charge = static_cast<Int_t>(traj->GetParticle()->Charge());
    } else {
      ND280Error("PDG number " << traj->GetPDGEncoding()
                               << " does not exist in database. Cannot fill "
                                  "Mass and Charge values.");
    }
  }
  // END LOOP OVER TRAJECTORIES

  // Sorty the clones array so trajectories are in order of ascending ID.
  fTrajectories->Sort();

  return true;
}

bool ND::TTruthTrajectoriesModule::SaveTraj(
    ND::TG4Trajectory* const traj) const {
  // Always save primary trajectories
  if (traj->GetParentId() == 0) {
    return true;
  }

  // Do not save trajectories if they are outside one of the primary detector
  // systems and the flag to just save detector trajectories is on
  if (fSaveOnlyP0DTrackerECALTrajectories) {
    std::vector<TG4TrajectoryPoint>::const_iterator points_it =
        traj->GetTrajectoryPoints().begin();
    std::vector<TG4TrajectoryPoint>::const_iterator points_end =
        traj->GetTrajectoryPoints().end();
    bool save = false;
    std::string volume_name;
    for (; points_it != points_end; ++points_it) {
      volume_name = points_it->GetVolumeName();
      if ((volume_name.find("/Tracker_") != std::string::npos) ||
          (volume_name.find("/P0D_") != std::string::npos) ||
          (volume_name.find("/SFG_") != std::string::npos) ||
          (volume_name.find("TOF_") != std::string::npos) ||
          (volume_name.find("/TopHAT_") != std::string::npos) ||
          (volume_name.find("/BottomHAT_") != std::string::npos) ||
          (volume_name.find("/BrlECal_") != std::string::npos) ||
          (volume_name.find("/DsECal_") != std::string::npos) ||
          (volume_name.find("/P0DECal_") != std::string::npos)) {
        save = true;
        break;
      }
    }

    if (!save) {
      return false;
    }
  }

  // Non-primary trajectories, even those in detectors are only saved if their
  // length is sufficiently long
  TLorentzVector length = traj->GetInitialPosition() - traj->GetFinalPosition();
  if ((length.Vect().Mag() < GetMinimumTrajectoryLengthToSave()) ||
      (traj->GetParentId() == 0)) {
    return false;
  }

  return true;
}

void ND::TTruthTrajectoriesModule::FillTrajectoryPoints(
    ND::TG4Trajectory* const traj,
    ND::TTruthTrajectoriesModule::TTruthTrajectory* trajToFill) {
  const ND::TG4Trajectory::Points& points = traj->GetTrajectoryPoints();

  for (ND::TG4Trajectory::Points::const_iterator pt = points.begin();
       pt != points.end(); ++pt) {
    int proc = 0;
    int subproc = 0;
#if !BEFORE_oaEvent(8, 11, 0)
    switch (pt->GetProcessType()) {
#ifndef SKIP_HADRONIC_TRAJECTORY_POINTS
      case ND::TG4TrajectoryPoint::kProcessHadronic: {
        proc = pt->GetProcessType();
        subproc = pt->GetProcessSubtype();
        break;
      }
#endif
#ifdef INCLUDE_TRANSPORT_TRAJECTORY_POINTS
    case ND::TG4TrajectoryPoint::kProcessTransportation: {
        proc = pt->GetProcessType();
        subproc = pt->GetProcessSubtype();
        break;
      }
#endif
#ifdef INCLUDE_ELECTROMAGNETIC_TRAJECTORY_POINTS
    case ND::TG4TrajectoryPoint::kProcessElectromagnetic: {
        proc = pt->GetProcessType();
        subproc = pt->GetProcessSubtype();
        break;
      }
#endif
#ifdef INCLUDE_DECAY_TRAJECTORY_POINTS
    case ND::TG4TrajectoryPoint::kProcessDecay: {
        proc = pt->GetProcessType();
        subproc = pt->GetProcessSubtype();
        break;
      }
#endif
      default: {
#ifdef FORCE_ALL_TRAJECTORY_POINTS
          proc = pt->GetProcessType();
          subproc = pt->GetProcessSubtype();
#endif
          break;
      }
    }
#else
    // ROOT TObject bit field may be set to indicate interaction type.
    // Bits 14 to 23 are available.  The definitions are.
    //
    //   15: Inelastic hadronic interaction
    //   16: Elastic hadronic interaction
    //   17: Hadronic charge exchange
    //
    //   20: Bremsstrahlung
    //   21: Pair production
    //   22: Compton scattering
    //
    if (pt->TestBit(BIT(15)) != 0) {
      proc = 4;
      subproc = 121;
    } else if (pt->TestBit(BIT(16)) != 0) {
      proc = 4;
      subproc = 111;
    } else if (pt->TestBit(BIT(17)) != 0) {
      proc = 4;
      subproc = 141;
    }
#endif
    if (proc == 0) {
      continue;
    }
    ND::TTruthTrajectoriesModule::TTruthTrajectoryPoint point;
    point.ProcessType = proc + 1000 * subproc;
    point.PositionX = pt->GetPosition().X();
    point.PositionY = pt->GetPosition().Y();
    point.PositionZ = pt->GetPosition().Z();
    point.PositionT = pt->GetPosition().T();
    point.MomentumX = pt->GetMomentum().X();
    point.MomentumY = pt->GetMomentum().Y();
    point.MomentumZ = pt->GetMomentum().Z();
    trajToFill->Points.push_back(point);
  }
}

void ND::TTruthTrajectoriesModule::FillTraces(
    ND::TG4Trajectory* const traj,
    ND::TTruthTrajectoriesModule::TTruthTrajectory* trajToFill) {
  // get the points from the trajectory
  std::vector<ND::TG4TrajectoryPoint> points = traj->GetTrajectoryPoints();
  if (points.size() <= 0) {
    ND280Error("Found TG4Trajectory with no trajectory points");
    return;
  }

  // Set the enumerations as if set to none of them
  eventAnalysisEnums::ESubdetector curr_subdet = eventAnalysisEnums::kNSubdetectors;
  eventAnalysisEnums::ESubdetector prev_subdet = eventAnalysisEnums::kNSubdetectors;

  // used to manipulate values before being passed to TraceInActive
  // std::vector<bool> inactive;
  // used to store the previous point - used for exit position/momenta
  ND::TG4TrajectoryPoint prev_point;
  // these are used to fill the stopping positions (last_ are filled everytime,
  // but only used once)
  TLorentzVector last_pos;
  TVector3 last_mom;

  // Loop over all points in traj
  std::vector<ND::TG4TrajectoryPoint>::iterator pt_iter = points.begin();
  std::vector<ND::TG4TrajectoryPoint>::iterator pt_end = points.end();

  // Variables we'll be using in the loop
  TString curr_path;
  std::string path_str;
  bool is_in_active;
  TLorentzVector pos;
  TVector3 mom;
  TLorentzVector prev_pos;
  TVector3 prev_mom;

  ND280Debug("---Looping over new trajectory:");
  for (; pt_iter != pt_end; ++pt_iter) {
    // get the current geometry path
    curr_path = pt_iter->GetVolumeName();
    path_str = curr_path.Data();

    // convert the path into a subdetector enumeration
    curr_subdet = eventAnalysisEnums::PathToSubdetector(path_str);

    // find out if its in an active volume
    is_in_active = this->GetIsActive(*pt_iter, curr_subdet);

    ND280Debug("     -> is in active = "
               << is_in_active
               << " for det: " << eventAnalysisEnums::DetectorName(curr_subdet)
               << " from path:" << curr_path.Data());

    // Fill vectors if the current subdetector has changed or it is the first
    // point.
    if (curr_subdet != prev_subdet || pt_iter == points.begin()) {
      ND280Debug("     --- above is start of different sub det ---");

      trajToFill->TraceSubdetectors.push_back(static_cast<int>(curr_subdet));
      // this is just to keep the length of inactive correct, will assign
      // correct value later
      // inactive.push_back(false);
      trajToFill->TraceInActive.push_back(false);

      // alliases for current and previous points
      pos = pt_iter->GetPosition();
      mom = pt_iter->GetMomentum();
      prev_pos = prev_point.GetPosition();
      prev_mom = prev_point.GetMomentum();

      // these are used to fill the stopping positions (last_ are filled
      // everytime, but only used once)
      last_pos = pos;
      last_mom = mom;

      // fill entrances with current points
      trajToFill->TraceEntrancePosition.push_back(pos);
      trajToFill->TraceEntranceMomentum.push_back(mom);
      if (pt_iter != points.begin()) {
        // fill exits with previous points
        trajToFill->TraceExitPosition.push_back(prev_pos);
        trajToFill->TraceExitMomentum.push_back(prev_mom);
      }
    }  // end first point or changed subdetector?
    else {
      // If it has not exited the last detector then it has stopped so just get
      // the last point.
      last_pos = pt_iter->GetPosition();
      last_mom = pt_iter->GetMomentum();
    }

    // Record whether have been in an active volume whilst traversing the
    // current detector volume
    if (is_in_active) {
      // int current_volume_index = inactive.size()-1;
      int current_volume_index = trajToFill->TraceInActive.size() - 1;
      ND280Verbose(
          "--Setting current active vector to true for element "
          << current_volume_index << " and sub detector "
          << eventAnalysisEnums::DetectorName(eventAnalysisEnums::ESubdetector(
                 trajToFill->TraceSubdetectors.at(current_volume_index))));

      trajToFill->TraceInActive[current_volume_index] = true;
    }  // end is_in_active?

    // used as the previous points in the next iteration of loop
    prev_subdet = curr_subdet;
    prev_point = *pt_iter;

  }  // end iterator over G4Points

  // Fill final trajectory point position and momentum manually
  trajToFill->TraceExitPosition.push_back(last_pos);
  trajToFill->TraceExitMomentum.push_back(last_mom);

  // check to make sure every vector is the same length
  unsigned int size = trajToFill->TraceSubdetectors.size();
  if ((trajToFill->TraceInActive.size() != size) ||
      (trajToFill->TraceEntrancePosition.size() != size) ||
      (trajToFill->TraceExitPosition.size() != size) ||
      (trajToFill->TraceEntranceMomentum.size() != size) ||
      (trajToFill->TraceExitMomentum.size() != size)) {
    ND280Warn(
        "Sizes of Trace prefixed variables don't match. Don't trust the "
        "results!");
  }
}

// In principle you could write this method without the subdetector
// argument. But because this has to be done from string comparisons at
// the moment, sending the detector give a significant performance
// benefit.
bool ND::TTruthTrajectoriesModule::GetIsActive(
    const ND::TG4TrajectoryPoint& point,
    eventAnalysisEnums::ESubdetector det) const {
  switch (det) {
    case eventAnalysisEnums::kFGD1: {
      double dist = ND::TGeomInfo::Get().FGD().EdgeDistanceFGD1(
          point.GetPosition().X(), point.GetPosition().Y(),
          point.GetPosition().Z());
      if (dist >= 0.0)
        return true;
      else
        return false;
    } break;

    case eventAnalysisEnums::kFGD2: {
      double dist = ND::TGeomInfo::Get().FGD().EdgeDistanceFGD2(
          point.GetPosition().X(), point.GetPosition().Y(),
          point.GetPosition().Z());
      if (dist >= 0.0)
        return true;
      else
        return false;
    } break;

    case eventAnalysisEnums::kTPC1:
    case eventAnalysisEnums::kTPC2:
    case eventAnalysisEnums::kTPC3: {
      TString volume = point.GetVolumeName();
      // I'm not familiar with the TPC geometry. This is the test
      // which was in place before I got here, and I can't obviously
      // see what function to use from oaGeomInfo to get whether
      // or not a point is in the active region of a TPC.
      if (volume.Contains("Half"))
        return true;
      else
        return false;
    } break;

    case eventAnalysisEnums::kDsECal:
    case eventAnalysisEnums::kP0DECalBottom:
    case eventAnalysisEnums::kP0DECalLeft:
    case eventAnalysisEnums::kP0DECalTop:
    case eventAnalysisEnums::kP0DECalRight:
    case eventAnalysisEnums::kBrlECalBottom:
    case eventAnalysisEnums::kBrlECalLeft:
    case eventAnalysisEnums::kBrlECalTop:
    case eventAnalysisEnums::kBrlECalRight: {
      // Temporary hack to deal with Bug 526
      //
      // ND280 does not save sufficient trajectory points that
      // a particle traversing the detector will dfinitley have
      // nodes registered in the Active volume of the detector.
      // Although it is the Active volume the end-user is
      // interested in, we need to temporarily condier the entire
      // ECal volume to be of interest, or nothing will be saved
      // at all.
      //
      // ORIGINAL CODE:
      // TString volume = point.GetVolumeName();
      // if( volume.Contains("Active") ) return true;
      // else return false;
      //
      // TEMPORARY CODE:
      return true;
    } break;

    case eventAnalysisEnums::kP0D: {
      double dist = ND::TGeomInfo::Get().P0D().EdgeDistance(
          point.GetPosition().X(), point.GetPosition().Y(),
          point.GetPosition().Z());
      if (dist >= 0.0)
        return true;
      else
        return false;
    } break;

    default:
      return false;
      break;
  }
}

eventAnalysisEnums::EParticleCategory ND::TTruthTrajectoriesModule::GetCategory(
    const ND::TG4Trajectory* const traj) {
  eventAnalysisEnums::EParticleCategory category =
      eventAnalysisEnums::PDGInfoToCategory(traj->GetParticle());
  switch (category) {
    case eventAnalysisEnums::kChargedLepton:
      ++fNTrajLepton;
      break;
    case eventAnalysisEnums::kChargedBaryon:
      ++fNTrajBaryon;
      break;
    case eventAnalysisEnums::kChargedMeson:
      ++fNTrajMeson;
      break;
    case eventAnalysisEnums::kOtherCharged:
      ++fNTrajOtherCharged;
      break;
    case eventAnalysisEnums::kPhoton:
      ++fNTrajPhoton;
      break;
    case eventAnalysisEnums::kOtherNeutral:
      ++fNTrajOtherNeutral;
      break;
    default:
      ++fNTrajOther;
  }
  return category;
}

//
// Truth Trajectory
// =====================================================================
//

ND::TTruthTrajectoriesModule::TTruthTrajectory::TTruthTrajectory()
    : ID(-1),
      PDG(-1),
      Name("unset"),
      Category(),
      InitMomentum(-999999.9, -999999.9, -999999.9, -999999.9),
      InitPosition(-999999.9, -999999.9, -999999.9, -999999.9),
      FinalPosition(-999999.9, -999999.9, -999999.9, -999999.9),
      ParentID(-1),
      PrimaryID(-1),
      Mass(-999999.9),
      Charge(-999999) {
  TraceSubdetectors.clear();
  TraceInActive.clear();
  TraceEntrancePosition.clear();
  TraceExitPosition.clear();
  TraceEntranceMomentum.clear();
  TraceExitMomentum.clear();
}

Int_t ND::TTruthTrajectoriesModule::TTruthTrajectory::Compare(
    const TObject* obj) const {
  // Sorry about the casting, but I don't want to slow down sorting
  // by checking that everything compared with is of the correct
  // type. In any case, you get what you ask for if you try and
  // use this with anything other than another trajectory.

  if (this->ID <
      static_cast<const ND::TTruthTrajectoriesModule::TTruthTrajectory*>(obj)
          ->ID) {
    return -1;
  } else if (this->ID >
             static_cast<const ND::TTruthTrajectoriesModule::TTruthTrajectory*>(
                 obj)
                 ->ID) {
    return 1;
  } else  // they must be equal
  {
    return 0;
  }
}